Modelos estadísticos del VIH usando software libre

Juan Pablo Acuña González1, Dr. Octaviano Juárez Romero2, Dr. Santiago Marquina Benítez3, Dr. Juan Villagómez Méndez4, Dr. Miguel Apolonio Herrera Miranda4 and Dr. Virgilio Cruz Guzmán4

1Maestría en Métodos Estadísticos Aplicados. Universidad Autónoma de Guerrero
2Director. Maestría en Métodos Estadísticos Aplicados. Universidad Autónoma de Guerrero
3Co-Director. Maestría en Métodos Estadísticos Aplicados. Universidad Autónoma de Guerrero
4Asesor. Maestría en Métodos Estadísticos Aplicados. Universidad Autónoma de Guerrero

27 febrero 2023

Abstract

Este trabajo presenta una implementación estadística en shiny (Chang et al. 2022) para modelar sistemas en biología a través de modelos deterministas de ecuaciones diferenciales ordinarias (ODEs) así como incluir la modelación estocástica a través de ecuaciones diferenciales estocásticas (SDEs) que permitan investigar la estructura dinámica de las reacciones bioquímicas que subyacen al virus de inmunodeficiencia humana (VIH) tales como el ciclo de replicación del virus o la respuesta inmune durante la etapa asintomática de la infección por VIH. En las últimas tres décadas, el lenguaje de programación $\texttt{R}$ (R Core Team 2022) se ha convertido en la herramienta estándar en estadística al usar un lenguaje equilibrado de alto nivel junto a uno de los mejores soportes para simulación estocástica y análisis estadístico. Las grandes ventajas de $\texttt{R}$ frente a otras tecnologías son su naturaleza de software de código abierto completamente gratuito y su uso cada vez mayor en bioinformática y varias aplicaciones de biología de sistemas (Giorgi, Ceraolo, and Mercatelli 2022). Muchas librerias se han escrito usando el lenguaje de programación $\texttt{R}$ para poder simular y hacer inferencia sobre los parámetros de estos modelos. Sin embargo, la pronunciada curva de aprendizaje requerida para aplicar metodologías estadísticas en $\texttt{R}$ de manera efectiva hace que estas librerías sean potencialmente difíciles de usar. Esta aplicación permite la generación de trayectorias deterministas y estocásticas para un modelo de ecuaciones diferenciales con los parámetros iniciales especificados mediante la implementación de dos algoritmos principales: 1) El método de Euler para representar la modelación clásica determinista y, 2) El algoritmo de Gillespie (Gillespie 1977) para evaluar las trayectorias por simulación estocástica. Finalmente se ha creado una librería, $\texttt{randomverse}$, para poder trabajar con la app localmente y brindar documentación para este tipo de implementaciones.

Package

randomverse 0.1.0

1 Introducción

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 de la simulación e inferencia. $\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 computacional, 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 ciclo de vida del VIH en una célula T CD4 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 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 Antecedentes

La modelación estocástica de los fenómenos naturales es una herramienta indispensable en la comprensión de enfermedades complejas al permitir incorporar factores como los mecanismos de herencia y mutación, las causas de la senescencia o la aparición de infecciones virales fatales. La funcionalidad de estos procesos se basa en reacciones químicas en los organismos y las estructuras de las células vivas.

Las áreas superpuestas de biología molecular, bioquímica y genética se han estudiado a partir de un enorme conjunto de técnicas de modelación. En este trabajo se enfoca en abordar aplicaciones que utilizan la metodología de las reacciones bioquímicas para la representación de procesos biológicos.

Esto produce una representación de los sistemas biológicos como redes de reacciones bioquímicas acopladas que pueden simularse in silico mediante diferentes algoritmos dependiendo de las suposiciones hechas sobre la cinética subyacente.

Las reacciones bioquímicas generalmente se especifican mediante ecuaciones de reacción de la forma como en Ecuación (1).

\[\begin{equation} \begin{aligned} a_1 A_1 + . . . + a_k A_k → b_1 B_1 + . . . + b_l B_l, \end{aligned} \tag{1} \end{equation}\]

donde $k$ es el número de reactivos y $l$ es el número de productos. $A_i$ representa la $i$ésima molécula de reactivo y $B_j$ es la $j$ésima molécula de producto. $a_i$ es el número de moléculas de $A_i$ consumidas en un solo paso de reacción, y $b_j$ es el número de moléculas de $B_j$ producidas en un solo paso de reacción. Los coeficientes $a_i$ y $b_j$ se conocen como estequiometrías. Se supone que las estequiometrías son números naturales con el máximo común divisor igual a uno.

Los fenómenos biológicos se representan como modelos especificando una red bioquímica con una lista de reacciones correspondientes al sistema de interés. Además, en esta representación se necesitan la velocidad de cada reacción y las cantidades iniciales de cada especie que reacciona.

Cada reacción tiene asociada una ley de velocidad que cuantifica la propensión de la reacción a ocurrir y que depende de las cantidades actuales de reactivos disponibles. Con todos los elementos anteriores, las reacciones, las leyes de velocidad y las cantidades iniciales, se especifica el modelo y se puede simular dinámicamente en una computadora.

Una vez que se establece la representación analítica, se construye una representación gráfica del modelo para comprender las redes de reacción como diagramas de ruta. El modelo gráfico para un conjunto de reacciones bioquímicas acopladas utiliza un gráfico donde los nodos se dividen en un conjunto que representa las especies (nodos elípticos) y otro conjunto que representa las reacciones (nodos rectangulares). Los arcos de un nodo de especies a un nodo de reacción indican que esas especies son reactivos para esas reacciones y los arcos de un nodo de reacciones a nodos de especies indican que esas especies son productos de esas reacciones. Los pesos asociados con los arcos representan las estequiometrías asociadas con los reactivos y productos (se supone que los arcos no numerados tienen un peso de uno). Finalmente, hay un número entero de moléculas de cada uno de los nodos que representan una especie de un lugar (P) a una transición (T) que indica la abundancia de esa especie asociada con él en un momento dado. Un conjunto de números de moléculas en cualquier momento dado se conoce como la marca actual de la red (estado del proceso).

Para medir el cambio en un sistema de este tipo, se comparan en un formato tabular una matriz Pre que contiene los pesos de los arcos que van de los lugares a las transiciones y una matriz Post que contiene los pesos de los arcos de las transiciones a los lugares. El número de moléculas asociadas con cada especie disminuirá de acuerdo con los números de Pre y aumentará de acuerdo con los números de Post cada vez que se produzca una reacción. Evaluando la diferencia entre las matrices Post y Pre se obtiene el cambio de estado del sistema.

La representación tabular del modelo de general se muestra en la Tabla 1.

Table 1: Tabular representation of biological systems.
	Reactants (Pre)				Products (Post)
Reactions	$x_1$	$x_2$	$...$	$x_u$	$x_1$	$x_2$	…	$x_u$
$r_1$	$a_{11}$	$a_{21}$	…	$a_{u1}$	$b_{11}$	$b_{21}$	…	$b_{u1}$
$r_2$	$a_{12}$	$a_{22}$	….	$a_{u2}$	$b_{12}$	$b_{22}$	…	$b_{u2}$
…	…	…	…	…	…	…	…	…
$r_v$	$a_{1v}$	$a_{2v}$	…	$a_{uv}$	$b_{1v}$	$b_{2v}$	…	$b_{uv}$

Para dar una representación matricial de este efecto neto, la ecuación (2) introduce la matriz $A$ como la diferencia entre Post y Pre, que se denomina matriz de reacción. La transpuesta de esta matriz se denota $S$ y se llama matriz estequiométrica.

\[\begin{equation} \overline{M} = M + Sr. \tag{2} \end{equation}\]

2.1 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 (3).

\[\begin{equation} \frac{dX}{dt} = f(X), \tag{3} \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 (4)

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

con un pequeño incremento $∆t$ Ecuación (5) es obtenida

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

y reacomodando, Ecuación (6) 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{6} \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:

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)$.
Un $\Delta t$ de tiempo es fijada
El tiempo al siguiente evento es calculado como $t : = t + \Delta t$.
EL siguiente evento es calculado como $X = X + f(X(t)) \Delta t$.
Los valores actualizados para $X$ y $t$ se apilan en forma de lista.
Si $t < T_{max}$, el algoritmo regrea al paso $2$.
La salida, tanto de $X$ como de $t$ se encuadra como un arreglo de datos una vez que el bucle para.

2.2 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 (7) actualiza el estado del sistema.

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

Aquí,

\[\begin{equation} S = (Post − Pre)^T, \tag{8} \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 Ecuación (9).

\[\begin{equation} h_0 (x, c) ≡ \sum_{i=1}^v h_i (x, c_i ). \tag{9} \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:

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$ .
Una función de riesgo $h_i (x, c_i )$ es calculada para cada $i = 1, 2, . . . , v$ basada en el estado actual $x$.
La función de riesgo combinada $h_0 (x, c) ≡ \sum_{i=1}^v h_i (x, c_i )$ es obtenida.
Se usa una muestra de una cantidad aleatoria $Exp(h_0 (x, c))$ para calcular el tiempo alsiguiente evento $t'$.
El timepo se actualiza como $t : = t + t'$ .
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.
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$.
Los valores actualizados de $x$ y $t$ se apilan en la forma de una lista.
Si $t < T_{max}$, el algoritmo regresa al paso $2$.
La salida de $x$ and $t$ se enmarca en un erreglo de datos al terminar el bucle.

3 Métodología

3.1 Ciclo de vida del VIH-1 en una célula T CD4

El modelo en (Sazonov et al. 2021), especifica las reacciones bioquímicas que describen el ciclo de vida dentro de una célula T CD4. Consiste de las siguientes etapas.

Entrada viral:
- Unión del virión a los receptores CD4 (la glicoproteína viral gp120 se une a los receptores CD4 en la superficie de las células T);
- Unión al correceptor CCR5 o CXCR4;
- Fusión de la membrana del virus y la membrana celular, es decir, la nucleocápside no está recubierta y el ARN viral se inyecta en la célula.

La dinamica anterior se describe mediante las siguientes ecuaciones:

\[\begin{equation} \begin{aligned} \frac{dx_1}{dt} &= -(k_{bound} + d)x_1 \\ \frac{dx_2}{dt} &= k_{bound}x_1 - (k_{fuse}+d_{bound})x_2. \end{aligned} \tag{10} \end{equation}\]

dónde: $x_1 = [V_{free}]$ es el número de viriones libres fuera de la célula; $x_2 = [V _{bound}]$ es el número de viriones unidos a CD4 y al co-receptor. Los parámetros del modelo $k_{bound} = 3.1$ $h^{- 1}$ ; $k_{fuse} = 0.7$ $h^{- 1}$ ; $d = 0.38$ $h^{-1}$ ;$d _{bound} = 0.008$ $h^{-1}$ ; representan la tasa de unión del virión a la membrana de la célula T CD4+, la tasa de fusión del virión con la célula, la tasa de eliminación de los viriones maduros libres y la tasa de degradación de los viriones unidos, respectivamente.

Transcripción inversa
- Síntesis de ADN de cadena negativa a partir de ARN viral;
- Síntesis de ADN de cadena positiva;
- Formación de ADN de doble cadena.

Esta etapa está modelada por las siguientes ecuaciones:

\[\begin{equation} \begin{aligned} \frac{dx_3}{dt} &=k_{fuse}x_2 -(k_{RT} + d_{RNA_{cor}})x_3 \\ \frac{dx_4}{dt} &= k_{RT}x_3 - (k_{DNA_{t}}+d_{bound})x_4. \end{aligned} \tag{11} \end{equation}\]

dónde: $x_3 = [RNA_{cor}]$ es el número de moléculas de RNA genómico en el citoplasma; $x_4 = [DNA_{cor}]$ es el número de moléculas de ADN provirales sintetizadas por transcripción inversa. Los parámetros del modelo $k_{RT} = 0.43$ $h^{- 1}$ ; $k_{ADN_{t}} = 0.12$ $h^{-1}$ ; $d_{ARN_{cor}} = 0.21$ $h^{-1}$ ; $d_{ADN_{cor}} = 0.03$ $h^{-1}$ ; representan la tasa de transcripción inversa, la tasa de transporte del ADN desde el citoplasma al núcleo, la tasa de degradación del ARN en el citoplasma y la tasa de degradación del ADN en el citoplasma, respectivamente.

Integración. El cambio en la cantidad de ADN proviral en el núcleo y la cantidad de ADN integrado se modelan con las siguientes ecuaciones:

\[\begin{equation} \begin{aligned} \frac{dx_5}{dt} &= (k_{DNA_{t}}x_4 -(k_{int} + d_{DNA_{nuc}})x_5 \\ \frac{dx_6}{dt} &= k_{int}x_5 - d_{DNA_{int}}x_6. \end{aligned} \tag{12} \end{equation}\]

dónde: $x_5 = [DNA_{nuc}]$ es el número de moléculas de ADN en el núcleo; $x_6 = [DNA int]$ es el número de ADN integrado. Los parámetros del modelo $k_{int} = 0.14$ $h^{- 1}$ ; $d_{ADN_{nuc}} = 0.01$ $h^{-1}$ ; $d_{ADN_{int}} = 0.00002$ $h^{- 1 }$; representan la tasa de integración, la tasa de degradación del ADN en el núcleo y la tasa de degradación del ADN integrado en el cromosoma, respectivamente.

Transcripción. Esta etapa se describe mediante las siguientes ecuaciones:

\[\begin{equation} \begin{aligned} \frac{dx_7}{dt} &= f_{TR}x_6 - (k_{ssRNA_g} g_{Rev} + k_{eRNA_g} f_{Rev} + d_{RNA_{g}})x_7 \\ \frac{dx_8}{dt} &= k_{ssRNA_g} g_{Rev}x_7 -(k_{dsRNA_{ss}} g_{Rev} + k_{eRNA_{ss}} f_{Rev} + d_{RNA_{ss}})x_8 \\ \frac{dx_9}{dt} &= k_{dsRNA_{ss}} g_{Rev}x_8 - (k_{eRNA_{ds}} + d_{RNA_{ds}})x_9 \\ \frac{dx_{10}}{dt} &= k_{eRNA_g} f_{Rev}x_7 - (k_{tp,RNA} + d_{RNA_g})x_{10}\\ \frac{dx_{11}}{dt} &= k_{eRNA_{ss}} f_{Rev}x_8 - d_{RNA_{ss}}x_{11}\\ \frac{dx_{12}}{dt} &= k_{eRNA_{ds}}x_9 - d_{RNA_{ds}}x_{12}. \end{aligned} \tag{13} \end{equation}\]

dónde: $x_7 = [mRNA_g]$ es el número de moléculas de ARNm de VIH en el núcleo: g para genómico o de longitud completa; $x_8 = [mRNA_{ss}]$ es el número de moléculas de ARNm de VIH empalmadas individualmente (ss) en el núcleo; $x_9 = [mRNA_{ds}]$ es el número de moléculas de ARNm de VIH doblemente empalmados (ds) en el núcleo; $x_{10} = [mRNAc_g]$ es el número de moléculas de ARNm de VIH en el citoplasma: g para genómico o de longitud completa; $x_{11} = [mRNAc_{ss}]$ es el número de moléculas de ARNm de VIH empalmadas individualmente (ss) en el citoplasma; $x_{12} = [mRNAc_{ds}]$ es el número de moléculas de ARNm de VIH doblemente empalmados (ds) en el citoplasma. $f_{TR} = f_{TR} ( x_{16} ) = TR_{cell} + f_{Tat} ( x 16 ) TR_{Tat} ; f_{Tat} = f_{Tat} ( x 16 ) = x 16 / ( θ_{Tat} + x 16 ) ;$ $f_{Rev} = f_{Rev} ( x 17 ) = x 17 / ( θ_{Rev} + x_{17} )$ ; $g_{Rev}= g_{Rev} ( X 17 ) = 1 - β f_{Rev}$. Los parámetros del modelo: $TR_{cell} = 15$ $h^{-1}$ ; $TR_{Tat} = 1500$ $h^{1}$ ; $θ_{Tat} = 10^3$ mol.; $\theta_{Rev} = 7.7 × 10^4$ mol.; $β = 0.9$; $k_{ssRNA_{g}} = 2.4$ $h^{-1}$ ; $k_{eRNA_g} = 2.3$ $h^{-1}$ ; $k_{RNA_{ss}} = 2.3$ $h^{-1}$ ; $k_{dsRNA_{ss}} = 2.4$ $h^{- 1}$ ; $k_{eRNA_{ds}} = 4.6$ $h^{- 1}$ ; $k_{tp,RNA} = 2.8$ $h^{- 1 }$; $d_{RNA_g} = d_{RNA_{ss}}= d_{RNA_{c_{ds}}} = 0.12$ $h^{- 1 }$ ; representan la tasa de transcripción basal intrínseca de la célula, el nivel de transcripción inducido por la transactivación de Tat, el efecto inhibidor de Rev sobre las tasas de empalme, lo que implica su reducción de 1/ (1 - β) veces en el nivel de saturación de Rev, la tasa de empalme para el ARN del virus de longitud completa, la tasa de $[mRNA_{g}]$ exportado desde el núcleo, la tasa de $[mRNA_{ss}]$ exportado desde el núcleo, la tasa de empalme para el ARN del virus empalmado individualmente, la tasa de $[mRNA_{ds}]$ exportan desde el núcleo, la tasa de transporte de $[mRNA_g]$ a la membrana celular y las tasas de degradación de $[mRNA_i ]$ , $i ∈ \{ g, ss, ds \}$ , respectivamente.

Traducción. Se considera la cinética de las proteínas Gag-Pol, Gag, gp160, Tat y Rev. Se utiliza el siguiente conjunto de ecuaciones diferenciales:

\[\begin{equation} \begin{aligned} \frac{dx_{13}}{dt} &= k_{trans} f_{g, Gag-Pol} x_{10} - ( k_{tp, Gag-Pol} + d_{p, Gag-Pol} ) x_{13}\\ \frac{dx_{14}}{dt} &= k_{trans} f_{g,Gag} x_{10} − ( k_{tp,Gag} + d_{p,Gag} ) x_{14} \\ \frac{dx_{15}}{dt} &= k_{trans} f_{ss,gp160} x_{11} − ( k_{tp,gp160} + d_{p,gp160} ) x_{15}\\ \frac{dx_{16}}{dt} &= k_{trans} f_{ds,Tat} x_{12} − d_{p,Tat} x_{16}\\ \frac{dx_{17}}{dt} &= k_{trans} f_{ds,Rev} x_{12} − d_{p,Rev} x_{17} \end{aligned} \tag{14} \end{equation}\]

dónde: $x_{13} = [ P_{Gag-Pol} ]$ es el número de moléculas de proteína: Gag-Pol; $x_{14} = [ P_{Gag} ]$ es el número de moléculas de proteína: Gag; $x_{15} = [ P_{gp160} ]$ es el número de moléculas de proteína: gp160; $x_{16} = [ P_{Tat} ]$ es el número de moléculas de proteína: Tat; $x_{17} = [ P_{Rev} ]$ es el número de moléculas de proteína: Rev. Los parámetros del modelo $k_{trans} = 524$ $h^{-1}$ ; $f_{g, Gag-Pol} = 0.05$; $f_{g, Gag} = 0.95$; $f_{ss,gp160} = 0.64$; $f_{ds,Tat} = 0.025$; $f_{ds,Rev }= 0.2$ ; $k_{tp, Gag-Pol} = 2.8$ $h^{-1}$ ; $k_{tp, Gag} = 2.8$ $h^{- 1}$ ; $k_{tp,gp160} = 2.8$ $h^{- 1 }$; $d_{p, Gag-Pol} = 0.09$ $h^{- 1}$ ; $d_{p,Gag} = 0.09$ $h^{-1}$ ; $d_{p,gp160} = 0.02$ $h^{-1}$ ; $d_{p,Tat} = 0.04 $ $h^{- 1}$ ; $d_{p,Rev} = 0.07$ $h^{ 1}$ ; representan la tasa de traducción de mRNA a proteínas, $f _{ij}$ representan la fracción de $[mRNA_i]$ que codifica $[ P_j ]$ , $i ∈ \{ g, ss, ds \}$ , $j ∈ \{ Gag-Pol, Gag, gp160, Tat, Rev \}$ . Siguiéndolos, los parámetros definen las tasas de transporte de proteína $[P_j]$ a membrana, $j ∈ \{ Gag-Pol, Gag, gp160 \}$ , y las tasas de degradación de las proteínas Gag-Pol, Gag, gp160, Tat y Rev, respectivamente.

Montaje, Brotación y Maduración. Las últimas etapas del ciclo de replicación del VIH-1 consisten en el transporte de las proteínas a la membrana plasmática celular. La brotación de los viriones sintetizados termina con su maduración. Las respectivas ecuaciones del modelo determinista dicen:

\[\begin{equation} \begin{aligned} \frac{dx_{13}}{dt} &= k_{trans} f_{g, Gag-Pol} x_{10} - ( k_{tp, Gag-Pol} + d_{p, Gag-Pol} ) x_{13}\\ \frac{dx_{18}}{dt} &= k_{tp,RNA} x_{10} − k_{comb} N_{RNA} f_c - d_{RNA_g} x_{18} \\ \frac{dx_{19}}{dt} &= k_{tp, Gag-Pol} x_{13} − k_{comb} N_{Gag-Pol} f_c - d_{mem, Gag-Pol} x_{19} \\ \frac{dx_{20}}{dt} &= k_{tp,Gag} x_{14} - k_{comb} N_{Gag} f_c - d_{mem,Gag} x_{20}\\ \frac{dx_{21}}{dt} &= k_{tp,gp160} x_{15} - k_{comb} N_{gp160} f_c - d_{mem,gp160} x_{21} \end{aligned} \tag{15} \end{equation}\]

donde: $x_{18} = [ RNA_{mem} ] , x_{19} = [ P_{mem, Gag-Pol} ] , x_{20} = [ P_{mem,Gag} ] , x_{21} = [ P_{mem,gp160} ]$ son, respectivamente, el número de ARN viral y las moléculas de proteína viral Gag-Pol, Gag y gp160 en la membrana. Los parámetros del modelo $k_{tp,ARN} = 2.8$ $h^{- 1}$ ; $k_{tp, Gag-Pol} = 2.8$ $h^{- 1}$; $k_{tp, Gag} = 2.8 h^{- 1}$ ; $k_{tp,gp160} = 2.8$ $h^{-1$ ; $k_{comb} = 8.0$ $h^{- 1}$ ;$N_{ARN} = 2$; $N_{Gag-Pol} = 250$; $N_{Gag} = 5000$; $N_{gp160} = 24$; representan las tasas de transporte de ARN y proteínas $[ P_j ]$ a la membrana, $j ∈ \{ Gag-Pol, Gag, gp160 \}$ , la tasa de incorporación de moléculas en complejos pre-virión, el número de transcritos de ARN viral en un nuevo virión, el número de moléculas de Gag-Pol en un nuevo virión, el número de moléculas de Gag en un nuevo virión y el número de moléculas de gp160 en un nuevo virión, respectivamente.

El parámetro del modelo $K_{V_{rel}}$ que ingresa a la función de ensamblaje del complejo del virión,

\[f_c=f_c(x_{18},...,x_{21})=x_{18}\cdot \frac{x_{19}}{x_{19}k_{V_{rel}}N_{Gag-Pol}}\cdot \frac{x_{20}}{x_{20}k_{V_{rel}}N_{Gag}}\cdot \frac{x_{21}}{x_{21}k_{V_{rel}}N_{gp160}}\] $K_{V_{rel}} = 10^3$ se refiere a la escala característica de la progenie viral por ciclo de replicación. Los parámetros del modelo $d_{ARN_g} = 0.12$ $h^{- 1}$ ; $d_{mem, Gag-Pol} = 0.004$ $h ^{- 1}$ ; $d_{mem, Gag} = 0.004$ $h^{-1}$ ; $d_{mem,gp160} = 0.014$ $h^{-1}$ ; representan las tasas de degradación de $[RNA_{mem}]$, la proteína Gag-Pol anclada a la membrana, la proteína Gag anclada a la membrana y la gp160 asociada a la membrana, respectivamente.

Para este modelo de replicación del VIH, la parametrización de tipo Michaelis-Menten se usa para limitar la tasa de ensamblaje de los viriones de la progenie por el componente proteico menos abundante de los tres considerados $[P_{mem, Gag-Pol}], [P_{mem, Gag} ] , [ P_{mem,gp160} ]$ . Los parámetros $K_{V_{rel}}$, $N_{Gag-Pol}$, $N_{Gag}$, $N_{gp160}$ especifican la cantidad de referencia de la progenie viral por célula y el número de moléculas de proteína Gag-Pol, Gag y gp160 requeridas para cada virión, respectivamente.

En la fase tardía de la replicación del virus, el ARN viral y las proteínas virales en la membrana se asocian con el complejo pre-virión y luego se combinan para generar un nuevo virión.

\[\begin{equation} \begin{aligned} \frac{dx_{22}}{dt} &= k_{comb} f_c - ( k_{bud} + d_{comb} ) x_{22}\\ \frac{dx_{23}}{dt} &= k_{bud} x_{22} - ( k_{mat} + d_{bud} ) x_{23}\\ \frac{dx_{24}}{dt} &= k_{mat}x_{23} - dx_{24} \end{aligned} \tag{16} \end{equation}\]

dónde: $x_{22} = [V_{pre-virión}]$ es el número de viriones en la membrana; $x_{23} = [V_{bud}]$ es el número de virus libres después de brotar de la célula; $x_{24} = [ V_{mat} ]$ es el número de viriones maduros fuera de la célula. Los parámetros del modelo $k_{comb} = 8$ $h^{-1}$ ; $k_{bud} = 2.0$ $h^{-1}$ ; $k_{mat} = 2.4$ $h^{- 1}$ ; $d _{comb} = 0.2$ $h^{-1}$ ; $d _{bud} = 0.38$ $h^{- 1}$ ; $d = 0.24$ $h^{- 1}$ ; representan la tasa de incorporación de moléculas en los complejos pre-virión, la tasa de gemación de nuevos viriones, la tasa de maduración, la tasa de degradación del complejo pre-virión ensamblado, la tasa de degradación de las partículas similares a virus inmaduras y la tasa de eliminación de viriones maduros libres, respectivamente.

3.2 Implementación en shiny

Las variables de entrada se especifican mediante parámetros de modelo. Estas incluyen el cuerpo de la función, los argumentos de la función y sus valores, las matrices de reacción previa y posterior del sistema, las dimensiones del modelo que especifican el número de reacciones y especies, así como los nombres de especies y del modelo. Con todos los parámetros especificados, se establece un horizonte temporal y las realizaciones correspondientes para un modelo se devuelven como salida, ya sea del sistema general o de una especie en particular.

3.2.1 Modelo determinista

Para el modelo de replicación se ingresan los valores de la Tabla como parámetros en la app.

Table 2: Shiny components of HIV-1 life-cycle in a T CD4 cell. Deterministic model.
Component	Input
Function Body	\(\texttt{-(kbound+d)x1,kboundx1-(kfuse+dbound)x2, kfusex2-(kRT+dRNAcor)x3,kRTx3-(kDNAt+dDNAcor)x4, kDNAtx4-(kint+dDNAnuc)x5,kintx5-dDNAintx6, (TRcell+(x16/(thetaTat+x16))TRTat)x6-(kssRNAg(1-betax17/(x17+thetaRev))+keRNAgx17/(x17+thetaRev)+dRNAg)x7,kssRNAg(1-betax17/(x17+thetaRev))x7-(kdsRNAss(1-betax17/(x17+thetaRev))+keRNAssx17/(x17+thetaRev)+dRNAss)x8,kdsRNAss(1-betax17/(x17+thetaRev))x8-(keRNAds+dRNAcds)x9,(keRNAgx17/(x17+thetaRev))x7-(ktpRNA+dRNAg)x10,(keRNAssx17/(x17+thetaRev))x8-dRNAssx11,keRNAdsx9-dRNAcdsx12,ktransfgGagPolx10-(ktpGagPol+dpGagPol)x13, ktransfgGagx10-(ktpGag+dpGag)x14,ktransfssgp160x11-(ktpgp160+dpgp160)x15,ktransfdsTatx12-dpTatx16,ktransfdsRevx12-dpRevx17, ktpRNAx10-kcombNRNAx18x19/(x19+KVrelNGagPol)x20/(x20+KVrelNGag)x21/(x21+KVrelNgp160)-dRNAgx18,ktpGagPolx13-kcombNGagPolx18x19/(x19+KVrelNGagPol)x20/(x20+KVrelNGag)x21/(x21+KVrelNgp160)-dmemGagPolx19,ktpGagx14-kcombNGagx18x19/(x19+KVrelNGagPol)x20/(x20+KVrelNGag)x21/(x21+KVrelNgp160)-dmemGagx20,ktpgp160x15-kcombNgp160x18x19/(x19+KVrelNGagPol)x20/(x20+KVrelNGag)x21/(x21+KVrelNgp160)-dmemgp160x21, kcombx18x19/(x19+KVrelNGagPol)x20/(x20+KVrelNGag)x21/(x21+KVrelNgp160)-(kbud+dcomb)x22,kbudx22-(kmat+dbud)x23,kmatx23-d*x24}\)
Function Arguments	$\texttt{x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,kbound,kfuse,d,dbound,kRT,dRNAcor,kDNAt,dDNAcor,kint,dDNAnuc,dDNAint,TRcell,TRTat,thetaTat,thetaRev,beta,kssRNAg,keRNAg,keRNAss,kdsRNAss,keRNAds,ktpRNA,dRNAg,dRNAss,dRNAcds,ktrans,fgGagPol,fgGag,fssgp160,fdsTat,fdsRev,ktpGagPol,ktpGag,ktpgp160,dpGagPol,dpGag,dpgp160,dpTat,dpRev,kcomb,NRNA,NGagPol,NGag,Ngp160,KVrel,dmemGagPol,dmemGag,dmemgp160,kbud,kmat,dcomb,dbud}$
Initial Markings with Hazards	$\texttt{4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.1,.7,.38,.0008,.43,.12,.21,.03,.14,.001,.00002,15,1500,10000,77000,0.9,2.4,2.3,2.3,2.4,4.6,2.8,0.12,0.12,0.12,524,0.05,0.95,0.64,0.025,0.2,2.8,2.8,2.8,0.09,0.09,0.02,0.04,0.07,8,2,250,5000,24,10000,0.004,0.004,0.014,2,2.4,.52,0.38}$
Pre matrix	\(\texttt{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}\)
Post matrix	\(\texttt{-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, -1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,-1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,-1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,-1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1}\)
Reactions	52
Species	24
Species names	Vfree,Vbound,RNAcor,DNAcor,DNAnuc,DNAint,mRNAg,mRNAss,mRNAds,mRNAcg,mRNAcss,mRNAcds,PGag-Pol,PGag,Pgp160,PTat,PRev,RNAmem,PmemGagPol,PmemGag,Pmemgp160,Vprevirion,Vbud,Vmat
Model name	HIV-1 life cycle

La solución del sistema determinista (10) - (16) con condición inicial $[ V_{free} ]( 0 ) = 4$ se presenta en la Figura 1.

Figure 1: Solución numérica del modelo determinista de la replicación del VIH-1

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] kableExtra_1.3.4 BiocStyle_2.24.0
## 
## loaded via a namespace (and not attached):
##  [1] bslib_0.4.0         compiler_4.2.2      BiocManager_1.30.18
##  [4] jquerylib_0.1.4     highr_0.9           tools_4.2.2        
##  [7] digest_0.6.29       jsonlite_1.8.0      evaluate_0.16      
## [10] lifecycle_1.0.3     viridisLite_0.4.1   rlang_1.0.6        
## [13] cli_3.4.0           rstudioapi_0.14     yaml_2.3.5         
## [16] xfun_0.33           fastmap_1.1.0       httr_1.4.4         
## [19] stringr_1.4.0       knitr_1.40          xml2_1.3.3         
## [22] sass_0.4.2          systemfonts_1.0.4   webshot_0.5.2      
## [25] svglite_2.1.0       glue_1.6.2          R6_2.5.1           
## [28] rmarkdown_2.16      bookdown_0.29       magrittr_2.0.3     
## [31] scales_1.2.1        htmltools_0.5.4     rvest_1.0.0        
## [34] colorspace_2.0-3    stringi_1.7.8       munsell_0.5.0      
## [37] cachem_1.0.6

References

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.

Giorgi, Federico M., Carmine Ceraolo, and Daniele Mercatelli. 2022. “The r Language: An Engine for Bioinformatics and Data Science.” Life 12 (5). https://doi.org/10.3390/life12050648.

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.

R Core Team. 2022. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Sazonov, Igor, Dmitry Grebennikov, Andreas Meyerhans, and Gennady Bocharov. 2021. “Markov Chain-Based Stochastic Modelling of HIV-1 Life Cycle in a CD4 t Cell.” Mathematics 9 (17): 2025. https://doi.org/10.3390/math9172025.

Wilkinson, Darren. 2018. Smfsb: Stochastic Modelling for Systems Biology. https://CRAN.R-project.org/package=smfsb.

	Reactants (Pre)				Products (Post)
Reactions	\(x_1\)	\(x_2\)	\(...\)	\(x_u\)	\(x_1\)	\(x_2\)	…	\(x_u\)
\(r_1\)	\(a_{11}\)	\(a_{21}\)	…	\(a_{u1}\)	\(b_{11}\)	\(b_{21}\)	…	\(b_{u1}\)
\(r_2\)	\(a_{12}\)	\(a_{22}\)	….	\(a_{u2}\)	\(b_{12}\)	\(b_{22}\)	…	\(b_{u2}\)
…	…	…	…	…	…	…	…	…
\(r_v\)	\(a_{1v}\)	\(a_{2v}\)	…	\(a_{uv}\)	\(b_{1v}\)	\(b_{2v}\)	…	\(b_{uv}\)