Una aplicación web para la educación a distancia en biomatemáticas al nivel media superior.

Juan Pablo Acuña González1* and Dra. María del Pilar Alonso Reyes.2

1Maestría en Docencia para la Educación Media Superior. Facultad de Estudios Superiores Acatlán, Universidad Nacional Autónoma de México.
2Tutora Principal. Facultad de Ciencias, Universidad Nacional Autónoma de México.

*jp_6@ciencias.unam.mx

Abstract

Se sitúa el contexto actual en México de la Educación a Distancia (EaD) y de la trascendencia de las Tecnologías de Información y Comunicación (TIC’s) en la Educación Media Superior (EMS) a partir de la Reforma Integral de la Educación Media Superior (RIEMS), ya que sintetiza los cambios y las fuerzas principales que rigen la educación a nivel nacional. Con esto se introduce una apĺicación en shiny (Chang et al. 2022) para atacar problemas en biomatemáticas mediados por tecnología a nivel EMS usando tanto una modelación determinista como estocástica.

1 Educación a distancia en México

Como principal razón o justificación de la Educación a Distancia (EaD) se encuentra el uso de herramientas de software libre que permiten el desarrollo de tecnología para resolver problemas dentro de la sociedad mexicana. La fuerte tendencia hacia las actividades asincrónicas y la generación de contenido educativo permiten una mayor cobertura que la educación tradicional. Particularmente, se vuelve fundamental la educación a distancia para integrar currículas y planes de estudio neurodiversos. Como se menciona en García Palacios (2015), la educación tradicional e institucionalizada tiende a: “… deja fuera a quienes no pueden o no quieren sujetarse a tiempos, lugares y modos de aprender determinados”. Y en cambio, también en García Palacios (2015) se comenta que: “…la principal intención que anima a la educación a distancia es que la distancia no exista”.Así, la EaD propicia nuevas condiciones académicas y se aleja de las prácticas burocráticas escolares, lo que brinda la oportunidad de trabajar con perfiles de estudiantes más diversos que permiten construir una perspectiva más grande de la educación.

La riqueza y trascendencia histórica que ha marcado la EaD en México a a través del Centro de Actualización del Magisterio en la Ciudad de México (CAMCM), se puede mostrar en la siguiente línea temporal (CAMCM 2015).

• 1921: Se inicia el proceso de expansión de la matrícula en educación básica y normal.
• 1945: Creación del Instituto Federal de Capacitación del Magisterio (IFCM) que porcionó la formación profesional a más de cien mil profesores en servicio en escuelas primarias en todo el país, publicó más de doscientos títulos que forman parte de la Biblioteca Pedagógica del Magisterio y distribuyó diversos materiales impresos para las etapas de educación a distancia.
• 1971: Nace la Dirección General de Mejoramiento Profesional del Magisterio (DGMPM), como proyecto de formación continua del magisterio.
• 1972: En la Universidad Nacional Autónoma de México (UNAM) (Coordinación de Universidad Abierta y Educación a Distancia, 2011-2013) nace el Sistema Universidad Abierta en 1972, ahora llamado Sis- tema de Universidad Abierta y Educación a Distancia.
• 1974: Se funda el Sistema Abierto de Enseñanza del Instituto Politécnico Nacional.
• 1979: La UPN (Universidad Pedagógica Nacional) crea el Sistema de Educación a Distancia
• 1978: La institución da paso a la denominación de Dirección General de Capacitación y Mejoramiento Profesional del Magisterio. Detectó que quedaban algunos profesores de educación primaria sin formación profesional y proporcionó las licenciaturas en Educación Preescolar y en Educación Primaria. Se ofrece un programa de Nivelación Psicopedagógica a profesores de educación secundaria y la Licenciatura en Docencia Tecnológica para profesores que se desempeñan en la asignatura de Educación Tecnológica. Además, se desarrollaron las licenciaturas siguientes: Educación Secundaria por Televisión y Educación Normal en diferentes especialidades; así como la formación de profesor en: Educación Preescolar (Bilingüe Bicultural), Educación Preescolar rural, Educación primaria (Instructor Comunitario).
• 1980: Se establece el Sistema de Enseñanza a Distancia de la Universidad Veracruzana y el Telebachillerato de Veracruz.
• 1989: Se fusionan las direcciones de Educación Normal y de Capacitación y Mejoramiento Profesional del Magisterio en la institución denominada Dirección General de Educación Normal y Actualización del Magisterio en el Distrito Federal (DGENAMDF), de la que depende orgánicamente el Centro de Actualización del Magisterio en el Distrito Federal (CAMDF).
• 1997: El Instituto Tecnológico y de Estudio Superiores de Monterrey (ITESM) funda la Universidad Virtual.
• 2004: La UDG (Universidad de Guadalajara) dictaminó el Sistema de Universidad Virtual.
• 2007: Surge la Universidad Virtual de Guanajuato.
• 2008: Los Comités Interinstitucionales para la Evaluación de la Educación Superior (CIEES) evalúan las primeras licenciaturas a distancia.
• 2011: Surge la Universidad Virtual de Michoacán.
• 2012: Universidad Digital del Estado de México y también la Universidad Abierta y a Distancia de México.
• 2014: Fueron aceptados por el Consejo Nacional de Ciencia y Tecnología (Conacyt) los posgrados en modalidad a distancia.
• 2018: Cambia el nombre de atribución de Centro de Actualización en el Distrito Federal a Centro de Actualización en la Ciudad de México. Realiza las funciones de docencia, investigación y extensión, y difusión, y tiene bajo su responsabilidad los proyectos de actualización, capacitación y superación profesional del magisterio de educación básica.

Una de las principales ideas de la política educativa desde principio de siglo es el debate entre docentes acerca del impacto de la globalización en buscar la equidad y la inclusión para cada estudiante. El binomio entre globalización y TIC’s permite un acceso a la información de manera veloz y abrumadora por lo que el rol de la docencia es poder ayudar a discriminar esta información. Por otro lado, el factor social que funge la globalización en la educación es el de preparar para la vida en democracia y ejercer la ciudadanía entendiendo diferentes realidades a nivel global. Esta es una de las pautas más sobresalientes de las reformas educativas en los últimos años, ya que la educación mediada por tecnología se he vuelto el estándar para poder evaluar las mejoras y las deficiencias en el sistema y en los procesos de aprendizaje.

Con la llegada de la década de los veinte la EMS se transforma y rompe con el esquema de reproducción cultural para convertirse en un espacio de generación de autonomía. Se retoma el trabajo desarrollado desde la década de los noventa por la Educación Media Superior a Distancia (EMSAD) para poder hacer llegar los servicios educativos a la población en confinamiento derivado de la infección por coronavirus. Se han volteado las miradas hacia las ventajas de desarrollar un enfoque educativo híbrido desde la perspectiva y los recursos de la EMSAD, ya que este puede crear mejores estrategias para generar estudiantes independientes a través de la manipulación de paquetes pedagógicos y generación de contenido, tanto en forma sincrónica como en la modalidad asincrónica para ajustarse lo más posible a la población.

Al contar con una mayor implantación de sistemas híbridos se produce una menor necesidad de gastos en recursos humanos por lo que se pueden ampliar más eficientemente las redes educativas. También se disminuyen los costos de infraestructura pues se pueden emplear espacios que ya están funcionando con otras actividades (escuelas, oficinas gubernamentales, etc.). La generación de habilidades de auto-gestión vuelven a los sistemas híbridos más sostenibles ya que la eficiencia terminal mejora pues cada estudiante se responsabiliza de su propio aprendizaje (Maciel and Ángeles 1970). El reto es superar la curva entre el desinterés al interés por estudiar en cada estudiante, por lo que se vuelve una estrategia que va de la mano, el contar con una educación inclusiva que atienda a las necesidades y deseos de la juventud actual y sea capaz de preparar para las nuevas esquemáticas de consumo y producción de la sociedad del tercer milenio.

Una de las acciones interesantes de la RIEMS en 2021 para atacar estos nuevos paradigmas educativos es la propuesta de la opción educativa auto-planeada pues ofrece un mayor alcance para obtener el certificado de bachillerato para personas que por razones diversas no pudieron seguir en la modalidad escolarizada. Es una propuesta en la modalidad mixta con un tiempo mínimo de dos años y una mediación docente del 30% del plan de estudios, con lo que las habilidades auto gestivas se van formando con acompañamiento. Además es una opción bivalente que permite la posibilidad de continuar con estudios de nivel superior así como capacitar para el trabajo (“Sitio Web Oficial Del CBTIS No. 179-a,” n.d.).

Otro aspecto muy importante en la RIEMS 2021 es el Programa Atención de Planteles Federales de Educación Media Superior con Estudiantes con Discapacidad (PAPFEMS) que se enfoca en dos objetivos generales: 1) fortalecer la atención y educación que las personas con discapacidad merecen al acceder a la Educación Media Superior y 2) mejorar los servicios educativos que los planteles federales de EMS ofrecen a este sector vulnerable de la población. Los lineamientos de la RIEMS se enfocan en crear estrategias para identificar y reducir las barreras que limitan el acceso universal a la población estudiantil a través de la eliminación de prácticas de discriminación basándose en la valoración de la diversidad y la adaptación de las sistemáticas educativas a las habilidades y estilos de aprendizaje de cada estudiante. Estas acciones tienen seguimiento por parte de la llamada Contraloría Social, que principalmente ejerce acciones para aplicar apoyos económicos proporcionados por PAPFEMS y el contar con una herramienta informática donde la población beneficiada del programa pueda realizar actividades (“Sitio Web Oficial Del CBTIS No. 179-b,” n.d.).

La aplicación en shiny se ha documentado como la librería \(\texttt{randomverse}\) del software libre \(\texttt{R}\) (R Core Team 2022) y se ha alojado en github (https://github.com/jp-rgb/randomverse) para ofrecer una auto planeación a través de la documentación y cumplir con ser un paquete pedagógico disponible a la población en general ya que se puede emplear de manera local empleando la librería o de manera remota mediante el acceso a la url: https://randomverse.shinyapps.io/RV_int/. Así, el uso de la aplicación puede servir como apoyo ante una audiencia tradicional o como una herramienta para estudiar independientemente.

En las últimas tres décadas \(\texttt{R}\) se ha convertido en la herramienta estándar en estadística utilizando un equilibrio de un lenguaje de alto nivel con uno de los mejores soportes para simulación estocástica y análisis estadístico (Giorgi, Ceraolo, and Mercatelli 2022). La ventaja de \(\texttt{R}\) frente a otras tecnologías son su naturaleza de software de código abierto completamente gratuito y su fuerte uso en bioinformática y biología computacional, lo que facilita su integración en asignaturas de EMS.

La modelación de sistemas dinámicos en biomatemáticas a nivel EMS se ha desarrollado partiendo de metodologías deterministas clásicas como el método de Euler, pues sólo requiere de conocimientos previos de cálculo diferencial e integral (CUAED-UNAM 2020). El método de Euler se integra a la aplicación a través de la pestaña \(\texttt{Deterministic}\).

Además se incluye la pestaña \(\texttt{Stochastic}\) para implementar el algoritmo de simulación estocástica (Gillespie 1977). La pestaña \(\texttt{Summary}\) brinda las estadísticas resumen de la simulación y sus salidas gráficas son la base para ilustrar conceptos de experimentos aleatorios y convergencia así como la identificación de dsitribuciones de probabilidad.

2 Introducción a la modelación en biomatemáticas

En aplicaciones en biología se utilizan términos algebraicos para representar especies de interés como lo muestra la Tabla 1.

Table 1: Representación algebraica de especies biológicas

Especie	término
Gen	\(g\)
Ácido ribonucleico (ARN)	\(r\)
Proteína	\(P\)

La modelación matemática de aplicaciones en biología es importante pues la forma en que interactúan las especies de interés en cierto fenómeno están reguladas por una serie de reacciones bioquímicas las cuales tienen una representación algebraica o analítica.

Así, para expresar el proceso por el cual la información de un gen se utiliza para generar una copia de ARN se plantea la reacción expresada en la Ecuación (1).

\[\begin{equation} \begin{aligned} g \rightarrow g + r \,\ (\text{transcripción}) \end{aligned} \tag{1} \end{equation}\]

Análogamente, en la Ecuación (2) se expresa el proceso de convertir la información contenida en el ARN mensajero (ARNm) en proteínas.

\[\begin{equation} \begin{aligned} r \rightarrow r + P \,\ (\text{traducción}) \end{aligned} \tag{2} \end{equation}\]

Se puede hacer uso de las matrices para representar fenómenos en biología al asociar las reacciones bioquímicas de un proceso con las filas de la matriz, y asociar las especies involucradas en las reacciones con las columnas de la matriz.

Si por ejemplo se considera que hay un evento de trasncripción y uno de traducción, el número de especies presentes en la reacción se pueden analizar a través de dos matrices, una previa al evento (Tabla 2) y otra posterior al evento (Tabla 3).

Table 2: Representación matricial previa al evento (\(Pre\))

	Especies
Reacción	\(g\)	\(r\)	\(P\)
Tanscripción	1	0	0
Traducción	0	1	0

Table 3: Representación matricial posterior al evento (\(Post\))

	Especies
Reacción	\(g\)	\(r\)	\(P\)
Tanscripción	1	1	0
Traducción	0	1	1

Para iniciar un algoritmo que permita actualizar un sistema en biología con respecto al tiempo, se considera la siguiente secuencia de pasos:

• Algoritmo 1 (de Actualización)

1. Se toma el estado inicial del sistema como un vector, \(M_0\), que registra las cantidades iniciales de cada especie biológioca del sistema en \(t=0\).
2. Se calcula la matriz de efecto neto del sistema \(A=Post-Pre\)
3. Se contruye la matriz de estequiometrías, \(S=A^T\), como la transpuesta de \(A\)
4. Se toma el vector \(r\) que indica las reacciones que ocurren a tiempo \(i\)
5. Se calcula el valor del nuevo estado del sistema como \(\tilde{M}=M+Sr\)

En el contexto del ejemplo de transcripción y traducción de las Ecuaciones (1) y (2), se puede calcular el valor del sistema a partir del estado inicial como

\[\begin{equation} \begin{aligned} \tilde{M} &= M+Sr \\ &= M+A^Tr \\ &= \begin{pmatrix} 1 \\ 2 \\ 10 \end{pmatrix} + \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{pmatrix}^T \begin{pmatrix} 1 \\ 1 \\ \end{pmatrix} \\ &= \begin{pmatrix} 1 \\ 2 \\ 10 \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ 1 & 0 \\ 0 & 1 \\ \end{pmatrix} \begin{pmatrix} 1 \\ 1 \\ \end{pmatrix} \\ &= \begin{pmatrix} 1 \\ 2 \\ 10 \end{pmatrix} + \begin{pmatrix} 0 \\ 1 \\ 1 \\ \end{pmatrix} \\ &= \begin{pmatrix} 1 \\ 3 \\ 11 \\ \end{pmatrix} \end{aligned} \tag{3} \end{equation}\]

2.1 Enfoque determinista

La cinética determinista de los sistemas biológicos se rige por leyes de reacción representadas por sistemas de ecuaciones diferenciales. Cuando las ecuaciones diferenciales son demasiado complicadas para resolverlas analíticamente, la solución se puede examinar numéricamente integrando un sistema de Ecuaciones Diferenciales Ordinarias (EDO). El enfoque más simple y práctico emplea el método de Euler de primer orden.

La idea básica de integrar numéricamente un sistema de ecuaciones diferenciales que representan un sistema en biología comienza escribiendo el sistema como un vector de EDO como en la Ecuación (4).

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

donde \(X\) es un vector \(p\)-dimensional y \(f (·) : R^p → R^p\) es una función arbitraria (no-lineal) \(p\)-dimensional de \(X\). En las aplicaciones biológicas esta función está dada por, \(f(X) = Sr(X)\).

Empleando la definición de derivada, definida en la Ecuación (5)

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

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

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

y reacomodando, la Ecuación (7) da una froma simple de calcular \(X(t + ∆t)\) a partir de \(X(t)\). Iniciando con la condición inicial \(X(0)\), entonces \(X(∆t), X(2∆t), X(3∆t), . . .\) pueden calcularse y así obtener la dinámica completa del sistema.

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

Esta rutina se calcula en la pestaña \(\texttt{Deterministic}\) en función del horizonte temporal (con las especificaciones de parámetros previas) y se resume en los siguientes pasos:

• Algoritmo 2 (método de Euler)

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.

2.2 Simulación estocástica

Para llevar a cabo el estudio de simulación se toma un sistema con \(u\) especies \(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\) es \(X_t = (X_{1t} , X_{2t} , . . . , X_{ut})^T\). El número de 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 (8) actualiza el estado del sistema.

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

Aquí,

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

\[\begin{equation} h_0 (x, c) ≡ \sum_{i=1}^v h_i (x, c_i ). \tag{10} \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 secuencia de pasos para llevar a cabo esta rutina está implementada bajo la pestaña \(\texttt{Stochastic}\) y se resume como sigue:

• Algoritmo 3 (de Gillespie)

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 tiempo 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 índice 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\) y \(t\) se enmarca en un erreglo de datos al terminar el bucle.

3 Metodología

3.1 Diseño de la aplicación en shiny

La principal meta que se buscó al diseñar la aplicación fue contar con la mayor reactividad posible para ofrecer a la parte usuaria una interfaz potente para representar modelos básicos en biomatemáticas con suficiente libertad. La interfaz cuenta con un diseño despejado común dentro de estas aplicaciones. La disposición de los elementos de la aplicación se divide en: la sección que contiene un título dispuesto horizontalmente en la parte superior de la aplicación (Región A), a la izquierda un panel de barra lateral que contiene principalmente controladores que solicitan la entrada de parámetros (Región B), y a la derecha se encuentra el panel principal para mostrar visualizaciones y alternar entre tareas específicas mediante las pestañas (Región C). Este diseño se ilustra a continuación en la Figura 1.

Figure 1: Diseño de la aplicación shiny
Región A: encabezado. Región B: entrada de parámetros. Región C: salida gráfica.

Este diseño es relativamente simple de implementar, al mismo tiempo que proporciona una interfaz usuaria limpia e intuitiva.

Figure 2: Vista global de la aplicación
Salida gráfica para los parámetros del modelo y análisis de simulación especificados.

A continuación se describen los componenetes de la barra lateral y las opciones para la salida gráfica.

3.1.1 Parámetros del modelo

En este primer sub apartado se especifican los parámeteros del modelo que se desea estudiar. Las especificaciones a ingresar se pueden ver en la parte superior de la captura de pantalla de la Figura 3.

• Cuerpo de la función: Se refiere a las expresiones dadas por un sistema de EDO o al conjunto de transiciones que especifican una cadena de markov (CM), separando cada término con una coma y sin dejar espacios en blanco.
• Argumentos de la función: Se refiere a la especificación de los símbolos a usar para expresar cada variable y cada parámetro del modelo. Los símbolos especificados en este campo son los que se deben usar al especificar el cuerpo de la función. Cada símbolo se va separando con una coma y sin dejar espacios en blanco, iniciando con los símbolos para expresar a las variables del modelo y siguiendo con la especificación de los símbolos para especificar a los parámetros del modelo.
• Valores iniciales: En este campo se especifican los valores numéricos que cada símbolo en el argumento de la fucnión del modelo toma al inicio del sistema. Se ingresan valores separados por coma sin dejar espacios en blanco.
• Matriz Pre: En este campo se ingresa los valores de la transpuesta a la matriz previa al evento, siguiendo el orden de aparición de valores por columna. Los valores numéricos se van ingresando separados por coma sin dejar espacios en blanco.
• Matriz Post: En este campo se ingresa los valores de la transpuesta a la matriz posterior al evento, siguiendo el orden de aparición de valores por columna. Los valores numéricos se van ingresando separados por coma sin dejar espacios en blanco.
• Número de reacciones: En este campo se selecciona la cantidad de reacciones bioquímicas de las que se compone el sistema.
• Número de especies: En este campo se selecciona la cantidad de especies biológicas de las que se compone el sistema.
• Nombres de las especies: Como su etiqueta lo indica, en este campo se ingresan los nombres asignados a las especies biológicas dentros del sistema. Se ingresan los nombres separándolos por comas y sin dejar espacios en blanco.
• Nombre del modelo: Aquí se ingresa el nombre asignado al modelo que se va a estudiar. Se pueden incluir espacios en blanco.

Figure 3: Panel superior de la barra lateral
Introducción de los parámetros del modelo a simular.

3.1.2 Análisis de simulación

El segundo sub apartado en la barra lateral permite trabajar con el análisis de trrayectorias de simulación para un modelo especificado por una CM. Las especificaciones a ingresar se pueden ver en la parte inferior de la captura de pantalla de la Figura 4.

• Trayectorias: Aquí se ingresa la cantidad de trayectorias a simular.
• Seleccionar especie: Esta línea permite seleccionar una especie biológica dentro del sistema para analizar la evolución y comportamiento de sus trayectorias. El número de simulaciones es el especificado en la cantidad de trayectorias.
• Histograma: Esta opción permite visualizar la distribución para una especie seleccionada y una cantidad de trayectorias especificada mediante un histograma.
• Horizonte temporal: En este campo se especifica el horizonte de tiempo para simular el sistema (en unidades del modelo especificado).
• Facet: Esta opción permite visualizar el sistema a detalle por especie biológica.
• Graficar: Este botón refresca todas las especificaciones ingresadas, tanto de parámetros como de la simulación, y produce la salida gráfica en el panel principal.

Figure 4: Panel inferior de la barra lateral
Especificaciones de la simulación.

3.1.3 Modelo determinista

Esta es la opción gráfica al seleccionar la pestaña \(\texttt{Deterministic}\). Produce visualizaciones del comportamiento respecto al tiempo por especie biológica de algún sistema especificado en los parámetros del modelo de la barra lateral. Una vez ingresados las especificaciones, se indica el horizonte de tiempo y se presiona el botón para refrescar y graficar el modelo. Esta pestaña produce las trayectorias deterministas a la solución numérica de un sistema de EDO a través del método de Euler. La Figura 5 muestra la salida gráfica bajo el enfoque determinista y la Figura 6 muestra la visualización por especie.

Figure 5: Salida gráfica del moodelo determinista

Figure 6: Salida gráfica del moodelo determinista
Visualización por especie biológica.

3.1.4 Modelo estocástico

Bajo la pestaña \(\texttt{Stochastic}\) se producen visualizaciones sujetas a la incorporación de “ruido” en los modelos para representar mejor la variabilidad en bioología. Para esta representación, se mantienen los parámetros ya especificados en la barra lateral, pero se debe cambiar el cuerpo de la función al conjunto de transiciones que especifican un modelo de CM. Al refrescar y graficar se producen las trayectorias de estocásticas a la solución exacta o completa por simulación estocástica mediante el algoritmo de Gillespie. Las Figuras 7 y 8 muestran la salida gráfica bajo el enfoque estocástico.

Figure 7: Salida gráfica del moodelo estocástico

Figure 8: Salida gráfica del moodelo estocástico
Visualización por especie biológica.

3.1.5 Resumen de simulación

La tercer pestaña, \(\texttt{Summary}\), brinda las visualizaciones de resumen para las especificaciones de la generación de trayectorias. Se cuenta con la opción de mostrar el histograma de la distribución por cada especie biológica. La Figura 9 muestra la visualización de las trayectorias y el histograma dadas las especificaciones de la barra

Figure 9: Salida gráfica de la simulación
Gráficas de las trayectorias y del histograma de las realizaciones son producidas para un especie biológicaespecificada.

3.2 Red genética auto-regulatoria

Un punto de partida para describir fenómenos biológicos mediante el uso de ecuaciones de reacción y su teoría matemática asociada es describir las leyes naturales que subyacen en el proceso esencial de transferencia de información entre moléculas portadoras de información.

Las células del cuerpo humano se especializan en diferentes tareas, y todos estos diferentes tipos de células provienen de un embrión unicelular. Con cada división de esa célula, las instrucciones se transmiten a nuevas células.

Estas instrucciones se pueden codificar en una cadena que representa un polímero hecho de nucleótidos. Los cuatro nucleótidos de las moléculas de ADN, adenina (A), guanina (G), citosina (C) y timina (T), dispuestos en secuencias específicas, almacenan la información para el funcionamiento de los organismos vivos.

La ruta básica de la replicación del ADN consiste en que, si un gen se activa, se transcribe en ARN mensajero (ARNm), y si este gen codifica proteínas, el ARNm se traduce en una proteína. La transcripción se produce en el núcleo de los eucariotas y la traducción se produce en el citoplasma.

El gen es un concepto central para la codificación de la información de las proteínas, ya que todas las células vivas dependen de ellos porque son esenciales para las tareas estructurales y funcionales de los organismos.

A través de un proceso estocástico, los efectos predominantes a escala de redes genéticas y bioquímicas pueden modelarse para representar las características ruidosas de los sistemas biológicos. Los procesos involucrados en la expresión y regulación génica pueden simplificarse generando redes de autorregulación. Estas redes describen los mecanismos reguladores que determinan la cantidad requerida de proteínas producidas por cada célula a través de un conjunto de reacciones bioquímicas acopladas. Por lo tanto, la transcripción y la traducción se pueden ensamblar para construir redes genéticas. Los componentes de la Ecuación (11) resumen una red genética auto-regulatoria (AR) muy simple.

\[\begin{equation} \begin{aligned} g + P_2 \leftarrow \rightarrow g · P_2 \,\ (\text{represión}) \\ g \rightarrow g + r \,\ (\text{transcripción}) \\ r \rightarrow r + P \,\ (\text{traducción}) \\ 2P \leftarrow\rightarrow P_2 \,\ (\text{dimerisación}) \\ r \rightarrow ∅ \,\ (\text{degradación} \,\ \text{mRNA} ) \\ P \rightarrow ∅ \,\ (\text{degradación} \,\ \text{proteína} ). \end{aligned} \tag{11} \end{equation}\]

En estas ecuaciones, \(P\) representa una proteína, \(P_2\) el compuesto de dos de estas proteínas ( dímero), \(g\) un gen y \(r\) una transcripción de \(g\). El conjunto vacío, \(∅\), indica la degradación del producto de una reacción y el punto, \(\cdot\), representa el compuesto de dos reactivos. La aplicación antes mencionada es intrínsecamente estocástica.

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.

Por ejemplo, en la red AR hay ocho reacciones, ya que dos de ellas pueden ocurrir en ambas direcciones, proceso conocido como reversibilidad. 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. Además, el sistema de autorregulación se define con una cantidad inicial para cada una de las cinco especies biológicas involucradas (\(g·P_2\), \(g\), \(r\), \(P\) y \(P_2\)). 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 in silico.

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 (L) 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 el marcado actual de la red.

Por lo tanto el modelo gráfico de la red genética AR tiene la siguiente representación:

Figure 10: Reacciones de la red genética auto-regulatoria

La red de reacciones en la Figura 10 implican un bucle. Muestra la evolución de reacciones en particular cuando estas ocurren, así en cada etapa del proceso, una nueva red se genera con sus respectivos valores iniciales. Para medir el cambio en el sistema se empleará la matriz Pre que contiene los pesos de los arcos yendo de lugares a transiciones y la matriz Post que contiene los pesos de los arcos yendo de transiciones a lugares y se comparan en formato tabular. El número de moléculas asociadas con cada especie decrecerá conforme a los valores en Pre y crecerá de acuerdo a Post cada vez que ocurre una reacción. Al tomar la diferencia entre Post y Pre se obtiene el cambio de estado en el sistema.

3.3 Prueba de la aplicación

La meta de una aplicación con una finalidad educativa es la de contar con una herramienta práctica y funcional. La aplicación es bastante robusta pues alerta a la parte usuaria de cualquier error cometido en la captura de información. Esta aplicacióon fue probada con un modelo de regulación génica relativamente pequeño, pero que vincula las principales temáticas de biología en EMS para desarrollar la discusión alrededor de modelos básicos en biomatemáticas y su implementación en shiny. Las líneas de ingreso en texto para correr la aplicación se proveen en el material suplementario y son descritas en el Apéndice. En breve, aplicación fue probada con especificaciones que cumplieran:

• Delimitación por coma entre distintos elementos ingresados como parámetros.
• No dejar espacios en blanco entre los elementos.
• Delimitar un horizonte de tiempo para modelos estocásticos donde se generen probabilidades inválidas.
• Editar cada línea de ingreso y especificar la pestaña de salida gráfica, para luego refresacar la aplicación.

El formato de cada una de estas restricciones son archivos de cadenas de texto simple (.txt), por lo que su ingreso se hace directamente en cada widget o copiando y pegando las cadenas de texto. Para conveniencia, al implementar el caso de estudio de la red genética se reportan tablas con estos valores en formato .txt y el archivo completo del modelo se encuentra en el Apéndice, también como archivo de texto plano.

Después de verificar que cada valor de ingreso en l aplicación tenga el formato correcto, las salidas gráficas son utilizadas para analizar el modelo propuesto tomando como guía el análisis en Wilkinson (2020) para corroborar los resultados. Al evaluar este modelo mediante la implementación en shiny se obtienen resultados similares.

4 Resultados y discusión

4.1 Aplicación en shiny

La aplicación está completa y funciona según lo previsto. El código de las aplicaciones se colocó en línea en el servicio de alojamiento del repositorio de GitHub, donde está disponible gratuitamente para su descarga bajo una licencia MIT (lo que significa que no hay restricciones en el uso del código, y tampoco responsabilidad para el autor original del código relacionado con su uso). Las guía para la parte usuaria también se ha colocado en GitHub.

4.2 Uso de la aplicación

4.2.1 Modelo determinista

La modelación cinética de reacciones bioquímicas trata de entender la evolución temporal de un sistema especificado por un conjunto de reacciones. Reconsiderando el modelo de la Ecuación (11), aunque las ecuaciones de reacción capturan las interacciones clave entre las especies biológicas, por sí solas no son suficiente para determinar la dinámica completa del comportamiento del sistema. Para eso, se necesitan conocer las velocidades a las que ocurre cada una de las reacciones (junto con algunas condiciones iniciales adecuadas).

El modelo (11) anima a pensar en el número de de las cinco especies biológicas involucradas (\(g·P_2\), \(g\), \(r\), \(P\) y \(P_2\)) como números enteros, que pueden cambiar solo en cantidades discretas (enteros) cuando ocurre un evento de reacción. Sin embargo, el estudio de la cinética pensando en un escenario determinista, consiste en un conjunto de reacciones bioquímicas de cantidades macroscópicas con las especies biológicas que reaccionan en una solución. Allí, la cantidad de cada especie biológica se considera generalmente como una concentración que puede variar continuamente a medida que avanza la reacción. Convencionalmente, la concentración de una especie \(X\) se denota \([X]\).

Suele ocurrir como una ley empírica que la velocidad instantánea de una reacción es directamente proporcional a la concentración (a su vez directamente proporcional a la masa) de cada reactivo elevado a la potencia de su estequiometría. Entonces, para la red genética auto-regulatoria (AR), la segunda reacción procederá a una velocidad proporcional a \([g][P_2]\). En consecuencia, por efecto de esta reacción, \([g]\) y \([P_2]\) disminuirán a una razón instantánea \(k_2 [g][P_2]\) (donde \(k_2\) es la constante de proporcionalidad para esta reacción), \(k_2[g][P_2]\) se conoce como la ley de velocidad de la reacción, y \(k_2\) es la constante de velocidad. Considerando las ocho reacciones de la red AR, se puede escribir un sistema de EDO como en la Ecuación (12).

\[\begin{equation} \begin{aligned} \frac{d[g\cdot P_2]}{dt} = k_1 [g] [P_2] - k_2 [g\cdot P_2] \\ \frac{d[g]}{dt} = k_2[g\cdot P_2] - k_1[g][P_2] \\ \frac{d[r]}{dt} = k_3[g]-k_7[r] \\ \frac{d[P]}{dt} = 2k_6[P] - 2k_5[P]^2 + k_4[r] - k_8[P] \\ \frac{d[P_2]}{dt} = k_2[g\cdot P_2]-k_1[g][P_2] + k_5[P]^2 - k_6[P]. \end{aligned} \tag{12} \end{equation}\]

Se deben especificar las constantes de velocidad de reacción, \(k_i, \,\ i=1,...,8\) (medidas en las unidades apropiadas), así como las concentraciones iniciales de cada especie. Una vez hecho esto, toda la dinámica del sistema está completamente determinada y puede revelarse resolviendo el sistema de EDO, ya sea analíticamente o numéricamente. Se puede reescribir el sistema ODE en (12) en forma de matriz como en la Ecuación (13).

\[\begin{equation} \begin{aligned} \frac{d}{dt} \begin{pmatrix} [g\cdot P_2] \\ [g] \\ [r] \\ [P] \\ [P_2] \end{pmatrix} = \begin{pmatrix} 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 2 & 0 & -1 \\ -1 & 1 & 0 & 0 & 1 & -1 & 0 & 0 \end{pmatrix} \begin{pmatrix} k_1 [g] [P_2] \\ k_2 [g\cdot P_2] \\ k_3[g] \\ k_4[r] \\ k_5[P]^2 \\ k_6[P] \\ k_7[r] \\ k_8[P] \end{pmatrix} \end{aligned} \tag{13} \end{equation}\]

donde la matriz \(5 \times 8\) es la matriz de estequiometrías, S. Se escribe \(r([X]) = (r_1 ([X]), r_2 ([X]), . . . , r_v ([X]))^T\) para el vector de \(v\) reacciones diferentes, entonces el modelo EDO puede actualizarse coforme a le Ecuación (14).

\[\begin{equation} \begin{aligned} \frac{d}{dt} [X] = Sr([X]). \end{aligned} \tag{14} \end{equation}\]

La implementación en la aplicación requiere de la especificación de los parámetros del modelo mostrados en la Tabla 4.

Table 4: Componentes en shiny de la red genética auto-regulatoria. Modelo determinista.

Componente	Entrada
Cuerpo de la Función	\(\texttt{th1*x2*x5-th2*x1,th2*x1-th1*x5*x2,th3*x2-th7*x3,th4*x3+2*th6*x4-2*th5*x4**2-th8*x4,th2*x1-th1*x2*x5-th6*x4+th5*x4**2}\)
Argumentos de la Función	\(\texttt{x1,x2,x3,x4,x5,th1,th2,th3,th4,th5,th6,th7,th8}\)
Valores Iniciales	\(\texttt{0,1,2,10,12,1,10,0.01,10,1,1,0.1,0.01}\)
Matriz Pre	\(\texttt{0,1,0,0,1,1,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,2,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0}\)
matriz Post	\(\texttt{1,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0}\)
Reacciones	\(\texttt{8}\)
Especies	\(\texttt{5}\)
Nombres de las Especies	\(\texttt{g⋅P2,g,r,P,P2}\)
Nombre del Modelo	\(\texttt{Auto-regulatory genetic network}\)

Esto produce la representación gráfica de la Figura 11.

Figure 11: Versión determinista de la red genética auto-regulatoria

El enfoque determinista no logra capturar la naturaleza discreta y estocástica de la cinética bioquímica a bajas concentraciones. Dado que muchos procesos intracelulares implican reacciones a concentraciones extremadamente bajas, estos efectos estocásticos discretos suelen ser relevantes para los modelos de biología de sistemas, de ahí una de las principales razones de emplear modelos estocásticos.

4.2.2 Modelo estocástico

En una red estocástica en la que el estado (representado por el número de tokens en cada nodo la Figura 10 ) cambia dinámica y aleatoriamente eligiendo activaciones de eventos de una manera aleatoria. En la red genética AR el número de tokens en un nodo representa el número de moléculas de un tipo dado, entonces las tasas de disparos de eventos vienen dadas por leyes de tasas estocásticas, y el algoritmo de Gillespie se puede utilizar para determinar el tiempo hasta el próximo evento de activación y qué evento para disparar.

En la Tabla 5 se encuentran las equivalencias para cada tipo de reacción a considerar en la modelación estocástica.

Table 5: Tabla de equivalencias para las reacciones de un modelo estocástico.

Orden	Reacción	Equivalencia
0	\(R_i : ∅ → X\)	\(h_i(x,k_i)=c_i\)
1	\(R_i : X_j → ?\)	\(h_i(x,k_i)=c_ix_j\)
>1	\(R_i: X_{i_1} + ... + X_{i_n} → ?\)	\(h_i(x,k_i)=f(c_i,x_j,p_{ij})\)

En las reacciones de orden superior, la ley de tasa estocástica se define con \(f(c_i,x_j,p_{ij})=c_i\prod_{i=1}^n \binom{x_j}{p_{ij}}\), donde\(p_{ij}\) es el \(ij\)-ésimo elemento de la matriz Pre.

Considerando las equivalencias de la Tabla 5 en el modelo en (11), este conduce a las leyes de tasas estocásticas de la Ecuación (15).

\[\begin{equation} \begin{aligned} h_1(x,c_1) &= c_1gP_2 \\ h_2(x,c_2) &= c_2(g\cdot P_2) \\ h_3(x,c_3) &= c_3g \\ h_4(x,c_4) &= c_4r \\ h_5(x,c_5) &= c_5(0.5)P(P-1) \\ h_6(x,c_6) &= c_6P_2 \\ h_7(x,c_7) &= c_7r \\ h_8(x,c_8) &= c_8P \end{aligned} \tag{15} \end{equation}\]

Para iniciar el algoritmo de Gillespie se supone un estado inicial \(M = (0,1,2,10,12)^T\) y constantes de tasa estocástica \(c = (1,10,0.01,10,1,1,0.1,0.01)^T\). La salida del algoritmo de Gillespie consiste en una lista que contiene dos elementos. El primer elemento, es un vector de tiempos de eventos, y el segundo elemento, es una matriz cuyas filas representan el estado del sistema inmediatamente antes del tiempo de evento correspondiente. Esta matriz tiene una fila adicional, correspondiente al estado del sistema inmediatamente después del evento final simulado. Como se muestra en la Figura 12, una sola realización de este proceso se muestra para cada especie. Esta salida puede contrastarse con la cinética determinista correspondiente a la Figura 11. Está claro que aunque la solución estocástica sigue aproximadamente la tendencia de la solución determinista, no la sigue extrictamente, sino que es libre en cada punto del timpo de moverse a observaciones cercanas de manera estocástica.

Figure 12: Una simulación de la dinámica estocástica discreta de la red genética auto-regulatoria

Para mejorar la visualización, la Figura 13 muestra la realización del proceso por especie para entender mejor su comportamiento.

Figure 13: Vista por especies de la red genética auto-regulatoria

También es integral la implementación a través de la aplicación puesto que el úncio cambio en los parámetros del modelo es en el cuerpo de la función, el cual se modifica como en la Tabla 6.

Table 6: Componentes en shiny de la red genética auto-regulatoria. Modelo estocástico.

Componente	Entrada
Cuerpo de la Función	\(\texttt{th1*x2*x5,th2*x1,th3*x2,th4*x3,th5*0.5*x4*(x4-1),th6*x5,th7*x3,th8*x4}\)
Trayectorias	\(\texttt{1000}\)
Seleccionar especie	\(\texttt{P}\)
Horizonte temporal	\(\texttt{10}\)

Ahora que se ha convertido un modelo determinista a un modelo estocástico, se procede a analziar las visualizaciones y el resumen de la simulación de la salida de la aplicación, así como el análisis de sensibilidad e incertidumbre para la representación estocástica e introducir la idea de modelar poblaciones celulares.

Se debe considerar que se obtendrán realizaciones diferentes cada vez que se ejecute la simulación. Comparando las salidas de las Figuras 11 y 13 es claro que el comportamiento cualitativo es algo similar, pero que las fluctuaciones estocásticas son muy pronunciadas. Una vez más, debe enfatizarse que estas fluctuaciones son intrínsecas al sistema y no tienen nada que ver con el error de medición experimental.

Sin embargo, llevar a cabo pocas realizaciones no son realmente suficientes para dar una buena información sobre el rango de comportamiento que es probable que exhiba el modelo. Esta información suele obtenerse ejecutando el modelo de simulación muchas veces y resumiendo el resultado de manera sensata. Al enfocarse en la especie \(P\) (\(P\) y \(P_2\) están relacionadas de manera determinista, por lo que no es necesario considerar ambas), se comienza a comprender el rango de dinámica que exhibe ejecutando un número relativamente pequeño de simulaciones y visualizando el conjunto ensamble (trayectoria media y sus bandas de confianza) de la Figura 14 (derecha).

Figure 14: Histograma y trayectoria media de la especie P (proteina) junto a sus bandas de confianza aproximadas basadas en 1,000 realizaciones

Este enfoque consiste en realizar muchas ejecuciones de la simulación a través de la aplicación y resumir la distribución del nivel de la especie de interés mediante el cálculo de las estadísticas adecuadas (como la media y la varianza de la muestra). Nótese que aunque es el caso aquí que los resultados del análisis determinista son consistentes con la media del estudio cinético estocástico esto no es cierto en general (e incluso aquí, la media del proceso estocástico no es exactamente la solución determinista, pero proporciona una aproximación razonable). En general, el análisis determinista no proporciona información útil sobre el análisis cinético estocástico. Si la distribución marginal en cada punto de tiempo fuera normal (que claramente no es el caso, pero a menudo es una aproximación razonable), entonces se esperaría que más del \(99\%\) de las realizaciones se encuentren dentro de las \(3\) desviaciones estándar de la media.

Por lo tanto, parece razonable utilizar la media de la muestra más/menos \(3\) desviaciones estándar de la muestra como guía para el rango de valores probables en cada momento. Para una comprensión más detallada de la distribución de valores en puntos de tiempo particulares, es posible estudiar el conjunto completo de realizaciones en un momento dado. Por ejemplo, la Figura 14 muestra la distribución de probabilidad empírica para la distribución de \(P\) en el tiempo \(t = 7\).

Se toma en cuenta la naturaleza discreta de la distribución y el hecho de que solo son posibles los valores impares, pues el proceso se inicia en un número impar de moléculas y solo puede cambiar por un múltiplo de dos. También se considera lo “normal” que parece ser la distribución de realizaciones en la Fgura 14 (izquierda), lo que justifica el uso de un enfoque de “media más/menos 3 DE” para resumir el resultado. También, la Figura sugiere que el sistema no parece haber convergido a una distribución estacionaria en el tiempo \(t = 7\) (tendencia creciente) y, por lo tanto, no se puede hablar aún de alguna convergencia en el sistema o de estacionareidad en su distribución al tiempo \(t=7\).

Es importante tener en cuenta que, aunque el análisis analítico de procesos simples puede dar una idea de problemas más complejos, la clase de modelos en los que son posibles los enfoques analíticos es muy restringida (normalmente interesas sistemas con interacciones complejas entre varios mecanismos intrincados). Por lo tanto, el estudio computacionalmente intensivo basado en la simulación y el análisis estocástico puede resultar una alternativa más realista de obtener información sobre la dinámica de un sistema en general.

Los modelos en biomatemáticas son interesantes para estudiar desde un punto de vista estocástico, aunque pueden resultar algo insatisfactorios en el sentido de que, a menos que las concentraciones y los volúmenes involucrados sean realmente muy pequeños, las fluctuaciones estocásticas no son particularmente significativas. En estos casos, el tratamiento determinista continuo, aunque es claramente una aproximación a la verdad, en realidad capta razonablemente bien los aspectos más importantes de la dinámica. Sin embargo, para un procesos ruidoso como la expresión génica y su regulación en la red genética AR, sólo pueda haber un número (v.gr., menos de diez) de moléculas de cualquiera de las especies reactivas clave. Por esta dinámica pueden dominar las fluctuaciones estocásticas y los modelos pueden exhibir un comportamiento que es imposible predecir a partir del análisis determinista continuo asociado, debido a su intractabilidad analítica.

Para la red genética AR se cuenta con tres especies clave del modelo: \(P,P_2\) Y \(r\). Su dinámica volátil estocástica discreta es clara en la realización de la Figura 13. Los eventos de transcripción de ARN son comparativamente raros y aleatorios en su ocurrencia. El número de monómeros proteicos oscila enormemente entre \(10\) y \(30\) moléculas, y la cantidad de dímeros de proteínas salta abruptamente en momentos aleatorios y luego decae gradualmente. Mirando más de cerca los primeros \(800\) segundos de la misma realización (Figura 15), está claro que los saltos en los niveles de dímero de proteína coinciden con los eventos de transcripción de ARN.

Figure 15: Vista de una realización del proceso de la red genética AR durante los primero 800 segundos

Esto ilustra un punto importante con respecto a la variación estocástica en modelos complejos, pues a pesar de que hay una cantidad relativamente grande de moléculas de dímero de proteína, su comportamiento es fuertemente estocástico debido al hecho de que se ven afectadas por la cantidad de transcritos de ARN, y hay muy pocas moléculas de transcripción de ARN en el modelo. En consecuencia, incluso si el interés principal radica en una especie con un número relativamente grande de moléculas, un modelo determinista continuo no capturará adecuadamente su comportamiento si se ve afectado por una especie que puede tener un número pequeño de moléculas.

Está claro que incluso durante este corto período de tiempo, las fluctuaciones estocásticas son muy significativas. Para comprender esta variación con más detalle, se toma ahora el número de moléculas de \(P\) en el momento \(t = 17\), pero con sólo una molécula de la especie \(g\) en la marca inicial del proceos. Estas condiciones iniciales son interesantes (y para cualquier problema de conteo en general) para estudiar las fases iniciales de algún proceso y sólos cuente con información limitada de las especies que reaccionan en el sistema. Al ejecutar muchas simulaciones del proceso, se construye la distribución de probabilidad para el número de moléculas, y esto se muestra en la Figura 16 (derecha) (basado en 10,000 simulaciones).

Figure 16: Evolución temporal de las moléculas de P con un horizonte temporal de 17 segundos
El histograma muestra la distribución de probabilidad empírica al tiempo t=17 segundos, basada en 10,000 realizaciones.

Esto muestra claramente que existe una probabilidad casi uniforme de que no haya moléculas de \(P\) en el tiempo \(17\) (y la explicación más probable de esto es que todavía no habrá habido un evento de transcripción). La distribución está claramente lejos de ser normal, por lo que es poco probable que un resumen medio de la distribución de más/menos tres DE sea adecuado en este caso.

Una vez que un modelo presenta este comportamiento como este, es probable que haya cierta incertidumbre con respecto a algunos aspectos cuantitativos de la especificación del modelo. Esto podría ser, por ejemplo, la incertidumbre sobre las condiciones iniciales o las constantes de tasa estocástica adoptadas. En concretoco, se supone que existe un grado de incertidumbre con respecto al valor de la tasa de transcripción del gen, \(c_2\). En el modelo se especificó un valor de \(c_2 = 0.01\), pero se puede suponer que cualquier valor entre \(0.005\) y \(0.03\) es plausible. Por lo tanto, es natural querer investigar la sensibilidad de la dinámica del modelo ante esta especificación particular.

Para los modelos deterministas continuos, está bien establecida una metodología sofisticada para el análisis de sensibilidad. Sin embargo, las técnicas de esta área no se transfieren bien al paradigma de la modelación estocástica. Una de las principales motivaciones para pensar en la sensibilidad del modelo es el deseo de comprender la incertidumbre en la verdadera dinámica del proceso. Sin embargo, en el contexto de la modelación estocástico, la incertidumbre sobre la dinámica del proceso es parte integral de todo el enfoque. En consecuencia, incluso en el caso de una certeza completa con respecto a la estructura del modelo, las leyes de velocidad, las constantes de velocidad y las condiciones iniciales, la evolución temporal del proceso es incierta o estocástica. Una forma natural de incorporar la incertidumbre con respecto a los parámetros del modelo es adoptar una interpretación bayesiana subjetiva de la probabilidad. En el paradigma bayesiano, la incertidumbre con respecto a los parámetros del modelo no es fundamentalmente diferente de la incertidumbre con respecto a la evolución temporal del proceso debido a la dinámica cinética estocástica. Ya hemos construido mecanismos para manejar la incertidumbre en la dinámica del proceso tratando de comprender la distribución de probabilidad de los resultados, en lugar de simplemente observar una realización particular del proceso. No hay ninguna razón por la que estas distribuciones de probabilidad no deban incluir incertidumbre con respecto a los parámetros del modelo además de la incertidumbre inducida por la cinética estocástica.

5 Conclusión

A medida que las técnicas cuantitativas e intensivas en datos van ganando terreno en la práctica y en la enseñanza de áreas profesionales biológicas y de la salud, como lo es la modelación estocática, se convierte en una práctica estándar incrementar la tecnología que ayude a la población de estudiantes de EMS el obtener una formación inicial en estás áreas emergentes. \(\texttt{R}\), si bien es un lenguaje de programación potente y versátil, tiene una curva de aprendizaje pronunciada que limita su utilidad para muchas partes usuarias a este nivel. La aplicación presentada aquí elimina esta barrera y brinda el beneficio adicional de trabajar con visualizaciones dinámicas e interactivas. Dada la extensa colección de paquetes diseñados en \(\texttt{R}\) específicamente para su uso en la modelación estocástica y la enseñanza de las biomatemáticas, este enfoque puede usarse para aumentar su accesibilidad a un grupo más amplio y diverso de estudiante.

Una progresión natural sería crear aplicaciones que pudieran atacar las siguientes temáticas:

• Una aplicación o una modificación al \(\texttt{randomverse}\) que permita reducir la dimensionalidad de este tipo de modelos y trabajar con leyes de conservación en sistemas biológicos.
• Una aplicación o una modificación al \(\texttt{randomverse}\) que utilice un algoritmo más robusto y eficiente para la generación de trayectorias. La aplicación puede permitir alternar entre distintas metodologías e hipótesis en la modelación.
• Una tercera posible aplicación futura sería una que se alimente de datos de muestras de especies biológicas en algún sistema dinámico y lleve a cabo una estimación de los parámetros requeridos para llevar a cabo su representación en el \(\texttt{randomverse}\) y poder simular escenarios. Además, esta aplicación puede incluir una especificación más para hacer una estimación desde el paradigma bayesiano para capturar mejor la incertidumbre sobre los parámetros de un modelo basado en datos y experiencia de área.

Como se indicó anteriormente, es posible realizar todas estas tareas directamente en \(\texttt{R}\), pero el proceso de pasar de un conjunto de datos sin procesar al resultado final deseado no es nada sencillo. Al cambiar estas tareas de problemas de codificación a interacciones con una interfaz gráfica de usuario, shiny hace posible que estas iniciativas puedan implementarse al nivel EMS, al generar herramientas internas de mayor accesibilidad.

6 Apéndices

6.1 Apéndice A: Guía para la parte usuaria

6.1.1 Introducción

Esta aplicación web utiliza el paquete \(\texttt{shiny}\) del software libre \(\texttt{R}\) para contar con una interfaz gráfica en la que la parte usuaria pueda intteractuar con las componentes de la modelación en biomatemáticas desde el enfoque determinista como del estocástico. La aplicación produce las gráficas de las trayectorias de cada modelación permitiendo una alta reactividad al manipular los parámetros del modelo y produciendo visualizaciones de alto nivel, inclusive sin haberse familiarizado antes con el lenguaje \(\texttt{R}\). Además la programación en \(\texttt{shiny}\) da una flexibildad mayor que las rutinas de código en la programación tradicional. Específicamente, la facilidad que tiene la GUI para moverse entre modelaciones y entre los parámetros de un modelo en particular, hacen que el análisis en la aplicación sea inclusive más eficiente aún para la parte usuaria experta en su área.

6.1.2 Configuración de la estación de trabajo en R

6.1.2.1 Instalación de R, RStudio y librerías necesarias

Para poder correr la aplicación localmente es necesario contar con la instalación del ambiente estadístico \(\texttt{R}\) y el ambiente integrado de desarrollo (IDE) \(\texttt{RStudio}\). Las descargas para \(\texttt{R}\) se encuentran disponibles para plataformas Unix, Windows y Mac OS. Las instrucciones de instalación se pueden encontrar en las siguientes ligas:

• The R Project for Statistical Computing: https://www.r-project.org/
• RStudio IDE: https://www.rstudio.com/

Una vez que se cuenta con el software instalado, se necesitan instalar las siguientes librerías de \(\texttt{R}\):

• \(\texttt{shiny}\): herramientas para el desarrollo de aplicaciones web en \(\texttt{R}\)
• \(\texttt{tidyverse}\): manejo y visualización de datos (Wickham et al. 2019)

La instalación de librerías desde \(\texttt{RStudio}\) es relativamente simple:

1. Se abre la IDE de \(\texttt{RStudio}\)
2. Se da clic en la pestaña \(\texttt{Packages}\) para ver una lista de las librerías instaladas disponibles.
3. Todos los paquetes necesarios para ejecutar la aplicación se encuentran en el repositorio de \(\texttt{R}\), Comprehensive R Archive Network (CRAN). En consecuencia, se pueden instalar fácilmente desde el IDE de \(\texttt{RStudio}\) siguiendo esta secuencia de pasos:

1. Hacer clic en la pestaña \(\texttt{Install}\). Aparecerá una ventana emergente.
2. Escribir el nombre de la librería se que desea instalar.
3. Marcar la casilla etiquetada como \(\texttt{install dependencies}\) para asegurarse de que todos las librerías requeridas por la librería que se está solicitando también se instala automáticamente.

A continuación se proporciona una captura de pantalla que ilustra el proceso de instalación de librerías:

Figure 17: Instalación de librerías en R desde la IDE RStudio

6.1.2.2 Descarga de archivos de la aplicación

Todos los archivos de código en \(\texttt{R}\) necesarios para ejecutar la aplicación están disponibles públicamente en línea en GitHub. El trabajo se comparte bajo una Licencia MIT, lo que significa que puede usarse sin limitaciones, siempre que se atribuya la autoría y la autoría esté liberada de cualquier responsabilidad relacionada con uso del código.

Una vez descargados, los archivos de la aplicación pueden almacenarse en cualquier lugar del sistema de archivos local. Para que la aplicación funcione de manera más eficiente, los archivos descargados para la aplicación deben mantenerse juntos en una sola carpeta.\(\texttt{R}\) requiere que cada aplicación \(\texttt{shiny}\) se llame \(\texttt{app.R}\) para que se reconozca y se ejecute como una aplicación, sin necesidad de seleccionar y resaltar manualmente el código. Además, se incluyen los archivos que contienen material que la aplicación \(\texttt{shiny}\) necesita para poder ejecutarse. Mantener cada archivo \(\texttt{app.R}\) y sus carpetas asociadas en un directorio propio con un nombre informativo hace posible la organización de esta y otras aplicaciones \(\texttt{shiny}\). La siguiente captura de pantalla muestra un ejemplo de ruta el sistema local junto al archivo \(\texttt{app.R}\) de la aplicación y sus carpetas asociadas:

Figure 18: Organización de carpetas en aplicaciones shiny

6.1.2.3 Ejecutando la aplicación

Ya que se tiene la estación de trabajo configurada, se puede correr la app desde \(\texttt{RStudio}\). Para hacerlo, se abre \(\texttt{RStudio}\) y se abre el archibo \(\texttt{app.R}\) de la aplicación. Una vez que el script aparece en \(\texttt{RStudio}\), la aplicación se puede desplegar al hacer clic en el butón \(\texttt{Run app}\) como lo ilustra la Figura 19. Una vez que se haya mostrado la aplicación, se puede accecer a su funcionalidad exclusivamente a través de la GUI.

Figure 19: Despliegue de la aplicación
Se abre el archivo app.R en RStudio y se despliega la aplicaciónd esde RStudio al oprimir el botón Run app, ubicado en el panel en la esquina superior derecha (delineado en rojo en la imagen).

6.1.3 Descripción general de la aplicación

6.1.3.1 Barra lateral

Al iniciar la aplicación, se abre la GUI. La Figura 20 señala las diversas características de la especificación de los parámetros del modelo y la Figura 21 muestra las especificaciones del resumen de simulación. Los campos a ingresar son:

1. Cuerpo de la función.
2. Argumentos de la función.
3. Valores iniciales.
4. Matriz Pre.
5. Matriz Post.
6. Número de reacciones.
7. Número de especies.
8. Nombres de las especies.
9. Nombre del modelo.
10. Trayectorias.
11. Seleccionar especie.
12. Histograma.
13. Horizonte temporal.
14. Facet.
15. Graficar.

Figure 20: Widgets de la barra lateral especifican los parámeteros del modelo a estudiar: 1) Cuerpo de la función, 2) Argumentos de la función, 3) Valores iniciales, 4) Matriz Pre, 5) Matriz Post, 6) Número de reacciones, 7) Número de especies, 8) Nombres de las especies, 9) Nombre del modelo

Figure 21: Widgets inferiores de la barra lateral especifican el resumen de la simulación: 10) Trayectorias, 11) Seleccionar especie, 12) Histograma, 13) Horizonte temporal, 14) Facet, 15) Graficar

6.1.3.2 Panel principal

Una vez ingresadda cada especificación de los parámetros del modelo y/o del anaálisis de la simulación, se selecciona un horizonte temporal y se refrescal a aplicación con el botón de graficar. La Figura 22 muestra la captura de pantalla de las pestañas disponibles para pŕoducir una salida gráfica especificada con la aplicación. Las opciones disponibles son:

1. Modelación determinista
2. Modelación estocastica
3. Resumen de la simulación

Figure 22: La selección de pestaña determina la salida gráfica del panel principal: 1) Modelación determinista, 2) Modelación estocástica, 3) Resumen de la simulación

6.1.4 Corriendo una simulación

Cuando se tenga la aplicación con todas sus entradas en formato apropiado, se presiona el botón de graficar para refresacar la información de la aplicación yproducir las salidas gráficas correspondientes. Una vez iniciada la simulación, la aplicación puede tardar unos segundos en cargar si el horizonte temporal es muy amplio o el número de realizaciones es grande. El panel principal se mostrará en pausa antes de graficar el resultado pero se espera a que termine de llevar a cabo los cálculos, como en la Figura 23.

Figure 23: Aplicación en espera durante el procesamiento de valores grandes de los parámetros del modelo o de la simulación

Las gráficas se generan una vez que la aplicación termina su tiempo de espera. A medida que se seleccionan características diferentes en la modelación, la visualización se actualiza dinámicamente. Si alguna de las especificaciones de la barra lateral es errónea la aplicación identifica el error y lo notifica como un mensaje de error en el panel principal.

Figure 24: Mensaje de error al hacer una mala especificación en las líneas de entrada de la barra lateral

Para producir una nueva gráfica simplemente se corrige la o las líneas de entrada que producen el error y se vuelve a refrescar la aplicación. La aplicación nuevamente entrará en estado de espera si los parámetros son grandes o mostrará inmediatamente las visualizacione para parámetros menos restrictivos.

6.1.5 Formatos de ingreso permitidos

La aplicación acepta cadenas de texto simples ingresadas delimitando cada elemento de la cadena por coma y sin dejar espacios en blanco. La aplicación es integral en el sentido de que los parámetros del modelo especificados para un tipo de modelación se mantienen al cambiar de enfoque, salvo el cuerpo de la función de cada modelación, que debe cambiarse en su línea de ingreso respectivo, situarse en la nueva pestaña de modelación y hacer clic en el botón de graficar para actualizar los resultados.

Así mismo, se debe delimitar un horizonte de tiempo para modelos estocásticos donde no se generen probabilidades inválidas. Por último todo cambio, ya sea corrección de errores o actualización de información, se debe llevar a cabo previo a seleccionar una pestaña de modelación y luego mandar a refrescar la aplicación.

6.1.6 Interpretación general de las salidas gráficas

En general la salida gráfica por un simulador determinista representa una solución que se mantiene estable con los parámetros del modelo dados. Para cualquier otra réplica de este experimento la misma salida gráfica resultará una y otra vez.

Figure 25: Solución estable del simulador determinista para un horizonte temporal y las especificaciones de los parámetros del modelo fijos

En cambioo con el simulador estocástico manteniendo fijas ambas condiciones (tiempo y parámetros), la salida gráfica cambia, representando la naturaleza ruidosa intrínseca de estos prooocesos. Estas variaciones son útiles para entender los procesos de riesgo involucrados en el modelo estudiado, como se muestra las Figuras 26. y 27

Figure 26: Solución variable del simulador estocástico para un horizonte temporal y las especificaciones de los parámetros del modelo fijos
Primera simulación.

Figure 27: Solución variable del simulador estocástico para un horizonte temporal y las especificaciones de los parámetros del modelo fijos
Segunda simulación.

Análogamente, estas visualizaciones pueden presentarse por especie biológica. La interpretación es similar a cada caso anterior, sólo que esta salida produce una salida gráfica por especie. En las Figuras 28 y 29 se da un ejemplo.

Figure 28: Vista por especie del simulador estocástico para un horizonte temporal y las especificaciones de los parámetros del modelo fijos
Primera simulación.

Figure 29: Vista por especie del simulador estocástico para un horizonte temporal y las especificaciones de los parámetros del modelo fijos
Segunda simulación.

Por último, la salida gráfica de resumen produce el conjunto ensamble dada la especificación en el número de trayectorias, con una cinta superior para representar la banda de confianza superior y una conta inferior para la banda de confianza inferior. Junto a esta salida gráfica se muestra el histograma de la distribución emípírica de la especie seleccionada, mostrando un valor central o media y una dispersión hacia los lados o desviación estándar.

Figure 30: Resumen de la simulación para cierta especie con un horizonte temporal y las especificaciones de los parámetros del modelo fijos
La salida gráfica muestra las bandas de confianza y el histograma al mismo tiempo en el panel principal al especificar la casilla del histograma en la especificación de la simulación.

6.1.7 Solución de problemas y cómo obtener ayuda

Se han hecho todos los intentos para probar a fondo esta aplicación. Si las funciones no se comportan como se describe, es probable que la fuente más común del problema sea una mala especificación en alguna línea de ingreso con formato incorrecto. Se puede consultar la sección que muestran los formatos aceptados.

También es importante hacer clic en el botón de graficar una vez especificado por completo el modelo y/o la simulación con cada línea de ingreso apropiada para el formato de archivo que se está usando. La aplicación no procesará una cadena de texto que no esté separada por comas o con espacios en blanco.

Este software se pone a disposición del público para todos los usos, incluida la modificación. Por lo tanto, se ofrece sin garantías ni soporte proporcionado. Sin embargo, el autor del software puede ser contactado en jp_6@ciencias.unam.mx, y brindará asistencia cuando sea posible.

6.2 Apéndice B: Proceso de prueba de la aplicación

El siguiente modelo de prueba se puede ir ingresando en la aplicación al copiar y pegar cada especificación en sus parámetros. Cada elemento de una línea de ingreso va delimitado por coma y sin dejar espacios en blanco (salvo la línea del título, que puede incluir espacios en blanco).

Las siguientes líneas se van ingresando por orden en los parámetros del modelo:

\(\texttt{th1*x2*x5-th2*x1,th2*x1-th1*x5*x2,th3*x2-th7*x3,th4*x3+2*th6*x4-2*th5*x4**2-th8*x4,th2*x1-th1*x2*x5-th6*x4+th5*x4**2}\)

\(\texttt{x1,x2,x3,x4,x5,th1,th2,th3,th4,th5,th6,th7,th8}\)

\(\texttt{0,1,2,10,12,1,10,0.01,10,1,1,0.1,0.01}\)

\(\texttt{0,1,0,0,1,1,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,2,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0}\)

\(\texttt{1,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0}\)

\(\texttt{8}\)

\(\texttt{5}\)

\(\texttt{g⋅P2,g,r,P,P2}\)

\(\texttt{Auto-regulatory genetic network}\)

Si se considera una modelación estocástica, se consideran las siguientes líneas de ingreso en el análisis de la simulación, y la línea del cuerpo de la función en los parámetros se modifica.

\(\texttt{th1*x2*x5,th2*x1,th3*x2,th4*x3,th5*0.5*x4*(x4-1),th6*x5,th7*x3,th8*x4}\)

\(\texttt{1000}\)

\(\texttt{P}\)

\(\texttt{10}\)

Como puede verse, las líneas ingresadas utilizadas para la prueba representa un sistema mucho más pequeños de lo que sería un sistema en bioloogía más realista. Los mdoelos que describen la dinámica de una población o de algún organismo suelen contar con muchas más variables (se necesitan más ecuaciones para describir un fenḿeno más completo) y varios parámetros más a considerar (diferentes velocidades de reacción dentro del mismo proceso), sin embargo este trabajo presenta la unificación de la modelación son software libre y los desarrollos en biomatemáticas disponibles para una audiencia de nivel EMS. Además la aplicación analizada cierra la brecha entre la pedagogía tradicional en biología separada de la curricula cuantitativa y permite trabajar con conceptos clave en temas de biología célular y molecular para preparar para la educación superior tanto en áreas biológicas/salud como las cuantitativas, así como las áreas emergentes en esta transdisciplina.

6.3 Apéndice C: Código en R

El código en R se proporciona en este apéndice por completez, pero se ve mejor después de descargarlo de GitHub (consulte el Apéndice D: Material complementario).

# Juan P Acuña
# Student ID: 22253567
# Email: 22253567@gmail.com
# Copyright J Acuña 2023
# Submitted in partial fulfillment of the requirements
# of the degree MADEMS

### Load packages
library(shiny)
library(tidyverse)
library(ggplot2)
library(ggfortify)
library(magrittr)


ui <- fluidPage(

  titlePanel(h4("An Introduction to the randomverse")),

  sidebarLayout(
    sidebarPanel(h3("Model parameters"),
      textInput("text5", label = h6("Input function body (comma delimited, no blanks)"), value = "-(kbound+d)*x1,kbound*x1-(kfuse+dbound)*x2,
                kfuse*x2-(kRT+dRNAcor)*x3,kRT*x3-(kDNAt+dDNAcor)*x4,
                kDNAt*x4-(kint+dDNAnuc)*x5,kint*x5-dDNAint*x6,
                (TRcell+(x16/(thetaTat+x16))*TRTat)*x6-(kssRNAg*(1-beta*x17/(x17+thetaRev))+keRNAg*x17/(x17+thetaRev)+dRNAg)*x7,kssRNAg*(1-beta*x17/(x17+thetaRev))*x7-(kdsRNAss*(1-beta*x17/(x17+thetaRev))+keRNAss*x17/(x17+thetaRev)+dRNAss)*x8,kdsRNAss*(1-beta*x17/(x17+thetaRev))*x8-(keRNAds+dRNAcds)*x9,(keRNAg*x17/(x17+thetaRev))*x7-(ktpRNA+dRNAg)*x10,(keRNAss*x17/(x17+thetaRev))*x8-dRNAss*x11,keRNAds*x9-dRNAcds*x12,ktrans*fgGagPol*x10-(ktpGagPol+dpGagPol)*x13,
                ktrans*fgGag*x10-(ktpGag+dpGag)*x14,ktrans*fssgp160*x11-(ktpgp160+dpgp160)*x15,ktrans*fdsTat*x12-dpTat*x16,ktrans*fdsRev*x12-dpRev*x17,
                ktpRNA*x10-kcomb*NRNA*x18*x19/(x19+KVrel*NGagPol)*x20/(x20+KVrel*NGag)*x21/(x21+KVrel*Ngp160)-dRNAg*x18,ktpGagPol*x13-kcomb*NGagPol*x18*x19/(x19+KVrel*NGagPol)*x20/(x20+KVrel*NGag)*x21/(x21+KVrel*Ngp160)-dmemGagPol*x19,ktpGag*x14-kcomb*NGag*x18*x19/(x19+KVrel*NGagPol)*x20/(x20+KVrel*NGag)*x21/(x21+KVrel*Ngp160)-dmemGag*x20,ktpgp160*x15-kcomb*Ngp160*x18*x19/(x19+KVrel*NGagPol)*x20/(x20+KVrel*NGag)*x21/(x21+KVrel*Ngp160)-dmemgp160*x21,
                kcomb*x18*x19/(x19+KVrel*NGagPol)*x20/(x20+KVrel*NGag)*x21/(x21+KVrel*Ngp160)-(kbud+dcomb)*x22,kbud*x22-(kmat+dbud)*x23,kmat*x23-d*x24"),
      textInput("text6", label = h6("Input function arguments (comma delimited, no blanks)"), value = "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"),
      textInput("text3", label = h6("Input initial markings with hazards (comma delimited, no blanks)"), value = "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"),
      textInput("text4", label = h6("Input Pre matrix (comma delimited, no blanks)"), value = "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0"),
      textInput("text2", label = h6("Input Post matrix (comma delimited, no blanks)"), value = "-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"),
      selectInput('NumRow',label = h6('Reactions'), choices = c(2:100),selected = 52),
      selectInput('NumCol',label = h6('Species'), choices = c(2:100),selected = 24),
      textInput("text", label = h6("Input species names (comma delimited, no blanks)"), value = "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"),
      textInput("title", label = h6("Input model name"), value = "HIV in T CD4 life-cycle"),
      h3("Simulation summary"),
                       numericInput("sims", label = h6("Trajectories"),2,min = 2,max = 10000),                     textInput("spec", label = h6("Select species"), value = "P"),
                       #checkboxInput("do2", "Histogram", value = F),
                       #checkboxInput("do3", "Confidence bands", value = F),
                       #checkboxInput("donum1", "Make #1 plot", value = T),
                       checkboxInput("donum2", "Histogram", value = F),
                       #checkboxInput("donum3", "Confidence bands", value = F),
                       numericInput("time", label = h6("Time horizon"),10,min = 0,max = 10000),
      checkboxInput("facet", "Facet", value = F),
      submitButton("Plot")
    ),

    mainPanel(
      tabsetPanel(
        tabPanel("Deterministic", plotOutput("det")),
        tabPanel("Stochastic", plotOutput("stoch")),
                tabPanel("Summary",
                         #plotOutput("summ")}
                         splitLayout(cellWidths = c("40%", "30%", "30%"), plotOutput("summ"), plotOutput("histo"), plotOutput("conf")))

      )
    )
  )
)

server <- function(input, output, session) {

  text_2 <- reactive({

        input$text2
      })
      text_5 <- reactive({
        input$text5
      })
      text_6 <- reactive({
        input$text6
      })
      text_3 <- reactive({
        input$text3
      })
      text_4 <- reactive({
        input$text4
      })
      title_r <- reactive({
        input$title
      })
      Num_Row <- reactive({
        input$NumRow
      })
      Num_Col <- reactive({
        input$NumCol
      })
      time_h <- reactive({
        input$time
      })
      text_n <- reactive({
        input$text
      })
      text_n <- reactive({
        input$text
      })
      sims_t <- reactive({
        input$sims
      })
      spec_n <- reactive({
        input$spec
      })

      BM <- list()
      BT <-list()
  xmat = list();
  tvec=c();
  x1=c();
  h1=c();

  output$det<-renderPlot({
    BM <- list()
    BT <-list()
    xmat = list()
    tt = 0
    tvec=c();
    h1 <- eval(eval(parse(text = paste('substitute(c(',input$text5,'), setNames(as.list(c(',input$text3,')),strsplit("',input$text6, '",split = ",")[[1]]))', sep=''))))
    x1 <- as.numeric(strsplit(input$text3,split = ",")[[1]])[1:as.numeric(input$NumCol)]
    xmat[[1]] = x1
    tvec[1] = 0
    dt = 0.001
    tt = tt + dt
    x2 <- as.numeric(strsplit(input$text3,split = ",")[[1]])[(as.numeric(input$NumCol)+1):(as.numeric(input$NumCol)+as.numeric(input$NumRow))]
    x1 = x1 + h1*dt
    xmat[[2]] = x1
    tvec[2] = tt
    i=3
    # n=as.numeric(input$time)/dt
    #while(tt<input$time) {
      for (i in 3:input$time) {
      x1 <- paste(c(x1,x2),collapse = ",")
      h1 <- eval(eval(parse(text = paste('substitute(c(',input$text5,'), setNames(as.list(c(',x1,')),strsplit("',input$text6, '",split = ",")[[1]]))', sep=''))))
      x1 <- as.numeric(strsplit(x1,split = ",")[[1]])[1:as.numeric(input$NumCol)]
      x1 = x1 + h1*dt
      tt = tt + dt
      xmat[[i]] = x1
      tvec[i] = tt
      tt=tvec[i]
      i=i+1
    }
    xmat = t(as.data.frame(xmat,row.names = strsplit(input$text,split = ",")[[1]]))
    data_long<-gather(cbind(tvec,as.data.frame(ts(xmat))),species,molecules,colnames(xmat))
    p = ggplot(data_long)+geom_line(aes(x=tvec,y=molecules,color=species))+title(main=input$title)
    if(input$facet) p = ggplot(data_long)+geom_line(aes(x=tvec,y=molecules))+facet_wrap(~species,scales="free")+title(main=input$title)
    #autoplot(ts(xmat),facets = input$facet)
    p
  }
  )
  xmat = list();
  tvec=c();
  x1=c();
  h1=c();

  output$stoch<-renderPlot({
    xmat = list()
    tt = 0
    tvec=c();
    S <- t(matrix(as.numeric(strsplit(input$text2, split = ",")[[1]])-as.numeric(strsplit(input$text4, split = ",")[[1]]),ncol = as.numeric(input$NumCol) ,nrow = as.numeric(input$NumRow),byrow=TRUE))
    h1 <- eval(eval(parse(text = paste('substitute(c(',input$text5,'), setNames(as.list(c(',input$text3,')),strsplit("',input$text6, '",split = ",")[[1]]))', sep=''))))
    x1 <- as.numeric(strsplit(input$text3,split = ",")[[1]])[1:as.numeric(input$NumCol)]
    xmat[[1]] = x1
    tvec[1] = 0
    tt = tt + rexp(1, sum(h1))
    x2 <- as.numeric(strsplit(input$text3,split = ",")[[1]])[(as.numeric(input$NumCol)+1):(as.numeric(input$NumCol)+as.numeric(input$NumRow))]
    j1 <- sample(as.numeric(input$NumRow), 1, prob = h1)
    x1 = x1 + S[, j1]
    xmat[[2]] <- x1
    tvec[2] = tt
    i=3
    #while(tt<input$time) {
      for (i in 3:input$time) {
      x1 <- paste(c(x1,x2),collapse = ",")
      h1 <- eval(eval(parse(text = paste('substitute(c(',input$text5,'), setNames(as.list(c(',x1,')),strsplit("',input$text6, '",split = ",")[[1]]))', sep=''))))
      tt = tt + rexp(1, sum(h1))
      j1 <- sample(as.numeric(input$NumRow), 1, prob = h1)
      x1 <- as.numeric(strsplit(x1,split = ",")[[1]])[1:as.numeric(input$NumCol)]
      x1 = x1 + S[, j1]
      xmat[[i]] <- x1
      tvec[i] = tt
      tt=tvec[i]
      i=i+1
    }
    xmat = t(as.data.frame(xmat,row.names = strsplit(input$text,split = ",")[[1]]))
    data_long<-gather(cbind(tvec,as.data.frame(ts(xmat))),species,molecules,colnames(xmat))
    p = ggplot(data_long)+geom_line(aes(x=tvec,y=molecules,color=species))+
      labs(x = "Time")+title(main=input$title)
    if(input$facet) p = ggplot(data_long)+geom_line(aes(x=tvec,y=molecules))+facet_wrap(~species,scales="free")+
      labs(x = "Time")+title(main=input$title)
    #autoplot(ts(xmat),facets = input$facet)
    p
  }
  )

  pt2 <- reactive({
    if (!input$donum2) return(NULL)
    #qplot(rnorm(500),fill=I("blue"),binwidth=0.2,main="plotgraph2")
        BM <- list()
        xmat <- list()
        tt <- 0
        tvec <- c();
        S <- t(matrix(as.numeric(strsplit(input$text2, split = ",")[[1]])-as.numeric(strsplit(input$text4, split = ",")[[1]]),ncol = as.numeric(input$NumCol) ,nrow = as.numeric(input$NumRow),byrow=TRUE))
        h1 <- eval(eval(parse(text = paste('substitute(c(',input$text5,'), setNames(as.list(c(',input$text3,')),strsplit("',input$text6, '",split = ",")[[1]]))', sep=''))))
        x1 <- as.numeric(strsplit(input$text3,split = ",")[[1]])[1:as.numeric(input$NumCol)]
        xmat[[1]] <- x1
        tvec[1] <- 0
        tt <- tt + rexp(1, sum(h1))
        x2 <- as.numeric(strsplit(input$text3,split = ",")[[1]])[(as.numeric(input$NumCol)+1):(as.numeric(input$NumCol)+as.numeric(input$NumRow))]
        j1 <- sample(as.numeric(input$NumRow), 1, prob = h1)
        x1 <- x1 + S[, j1]
        xmat[[2]] <- x1
        tvec[2] <- tt
        i <- 3
        #while(tt < input$time) {
          for (i in 3:input$time) {
          x1 <- paste(c(x1,x2),collapse = ",")
          h1 <- eval(eval(parse(text = paste('substitute(c(',input$text5,'), setNames(as.list(c(',x1,')),strsplit("',input$text6, '",split = ",")[[1]]))', sep=''))))
          tt <- tt + rexp(1, sum(h1))
          j1 <- sample(as.numeric(input$NumRow), 1, prob = h1)
          x1 <- as.numeric(strsplit(x1,split = ",")[[1]])[1:as.numeric(input$NumCol)]
          x1 <- x1 + S[, j1]
          xmat[[i]] <- x1
          tvec[i] <- tt
          tt <- tvec[i]
          i <- i+1
          }
        xmat <- t(as.data.frame(xmat,row.names = strsplit(input$text,split = ",")[[1]]))
        BM <- list(a=as.integer(tail(xmat[,input$spec],1)))

        j <- 1
        for (j in 1:(as.numeric(input$sims)-1)){
          xmat <- list()
          tt <- 0
          tvec <- c();
          S <- t(matrix(as.numeric(strsplit(input$text2, split = ",")[[1]])-as.numeric(strsplit(input$text4, split = ",")[[1]]),ncol = as.numeric(input$NumCol) ,nrow = as.numeric(input$NumRow),byrow=TRUE))
          h1 <- eval(eval(parse(text = paste('substitute(c(',input$text5,'), setNames(as.list(c(',input$text3,')),strsplit("',input$text6, '",split = ",")[[1]]))', sep=''))))
          x1 <- as.numeric(strsplit(input$text3,split = ",")[[1]])[1:as.numeric(input$NumCol)]
          xmat[[1]] <- x1
          tvec[1] <- 0
          tt <- tt + rexp(1, sum(h1))
          x2 <- as.numeric(strsplit(input$text3,split = ",")[[1]])[(as.numeric(input$NumCol)+1):(as.numeric(input$NumCol)+as.numeric(input$NumRow))]
          j1 <- sample(as.numeric(input$NumRow), 1, prob = h1)
          x1 <- x1 + S[, j1]
          xmat[[2]] <- x1
          tvec[2] <- tt
          i <- 3
          #while(tt < input$time) {
            for (i in 3:input$time) {
            x1 <- paste(c(x1,x2),collapse = ",")
            h1 <- eval(eval(parse(text = paste('substitute(c(',input$text5,'), setNames(as.list(c(',x1,')),strsplit("',input$text6, '",split = ",")[[1]]))', sep=''))))
            tt <- tt + rexp(1, sum(h1))
            j1 <- sample(as.numeric(input$NumRow), 1, prob = h1)
            x1 <- as.numeric(strsplit(x1,split = ",")[[1]])[1:as.numeric(input$NumCol)]
            x1 <- x1 + S[, j1]
            xmat[[i]] <- x1
            tvec[i] <- tt
            tt <- tvec[i]
            i <- i+1
            }
          xmat <- t(as.data.frame(xmat,row.names = strsplit(input$text,split = ",")[[1]]))
          BM <- append(BM,list(a=as.integer(tail(xmat[,input$spec],1))))
          j <- j+1
          }
          #return(
        hist(unlist(BM),freq = FALSE, xlab = "# Molecules", main = paste("Histogram of" , spec_n()))
        #)
    })

  output$summ <- renderPlot({
    BM <- list();
    xmat <- list()
    tt <- 0
    tvec <- c();
    S <- t(matrix(as.numeric(strsplit(text_2(), split = ",")[[1]])-as.numeric(strsplit(text_4(), split = ",")[[1]]),ncol = as.numeric(Num_Col()) ,nrow = as.numeric(Num_Row()),byrow=TRUE))
    h1 <- eval(eval(parse(text = paste('substitute(c(',text_5(),'), setNames(as.list(c(',text_3(),')),strsplit("',text_6(), '",split = ",")[[1]]))', sep=''))))
    x1 <- as.numeric(strsplit(text_3(),split = ",")[[1]])[1:as.numeric(Num_Col())]
    xmat[[1]] <- x1
    tvec[1] <- 0
    tt <- tt + rexp(1, sum(h1))
    x2 <- as.numeric(strsplit(text_3(),split = ",")[[1]])[(as.numeric(Num_Col())+1):(as.numeric(Num_Col())+as.numeric(Num_Row()))]
    j1 <- sample(as.numeric(Num_Row()), 1, prob = h1)
    x1 <- x1 + S[, j1]
    xmat[[2]] <- x1
    tvec[2] <- tt
    i <- 3
    #while(tt < time_h()) {
      for (i in 3:input$time) {
      x1 <- paste(c(x1,x2),collapse = ",")
      h1 <- eval(eval(parse(text = paste('substitute(c(',text_5(),'), setNames(as.list(c(',x1,')),strsplit("',text_6(), '",split = ",")[[1]]))', sep=''))))
      tt <- tt + rexp(1, sum(h1))
      j1 <- sample(as.numeric(Num_Row()), 1, prob = h1)
      x1 <- as.numeric(strsplit(x1,split = ",")[[1]])[1:as.numeric(Num_Col())]
      x1 <- x1 + S[, j1]
      xmat[[i]] <- x1
      tvec[i] <- tt
      tt <- tvec[i]
      i <- i+1
      }
    xmat <- t(as.data.frame(xmat,row.names = strsplit(text_n(),split = ",")[[1]]))
    BM <- list(a=as.integer(xmat[,spec_n()]))

    #xmat = t(as.data.frame(xmat,row.names = strsplit("g⋅P2,g,r,P,P2",split = ",")[[1]]))
    #BT <- list(a=tvec)
    #BM <- list(a=as.integer(xmat[,4]))
    #BM <- list(a=as.integer(tail(xmat[,4],1)))
    j <- 1
    for (j in 1:(as.numeric(sims_t())-1)){
      xmat <- list()
      tt <- 0
      tvec <- c();
      S <- t(matrix(as.numeric(strsplit(text_2(), split = ",")[[1]])-as.numeric(strsplit(text_4(), split = ",")[[1]]),ncol = as.numeric(Num_Col()) ,nrow = as.numeric(Num_Row()),byrow=TRUE))
      h1 <- eval(eval(parse(text = paste('substitute(c(',text_5(),'), setNames(as.list(c(',text_3(),')),strsplit("',text_6(), '",split = ",")[[1]]))', sep=''))))
      x1 <- as.numeric(strsplit(text_3(),split = ",")[[1]])[1:as.numeric(Num_Col())]
      xmat[[1]] <- x1
      tvec[1] <- 0
      tt <- tt + rexp(1, sum(h1))
      x2 <- as.numeric(strsplit(text_3(),split = ",")[[1]])[(as.numeric(Num_Col())+1):(as.numeric(Num_Col())+as.numeric(Num_Row()))]
      j1 <- sample(as.numeric(Num_Row()), 1, prob = h1)
      x1 <- x1 + S[, j1]
      xmat[[2]] <- x1
      tvec[2] <- tt
      i <- 3
      #while(tt < time_h()) {
        for (i in 3:input$time) {
        x1 <- paste(c(x1,x2),collapse = ",")
        h1 <- eval(eval(parse(text = paste('substitute(c(',text_5(),'), setNames(as.list(c(',x1,')),strsplit("',text_6(), '",split = ",")[[1]]))', sep=''))))
        tt <- tt + rexp(1, sum(h1))
        j1 <- sample(as.numeric(Num_Row()), 1, prob = h1)
        x1 <- as.numeric(strsplit(x1,split = ",")[[1]])[1:as.numeric(Num_Col())]
        x1 <- x1 + S[, j1]
        xmat[[i]] <- x1
        tvec[i] <- tt
        tt <- tvec[i]
        i <- i+1
        }
      xmat <- t(as.data.frame(xmat,row.names = strsplit(text_n(),split = ",")[[1]]))
      BM <- append(BM,list(a=as.integer(xmat[,spec_n()])))
      j <- j+1
    }
    gill_tibb<-cbind(tvec,as.data.frame(BM))
    gill_long <- gather(gill_tibb,simulation,molecules,colnames(as.data.frame(BM)))
    gill_summ <- gill_long %>%
      group_by(tvec) %>%
      summarize(mean_molcs = mean(molecules))
    # gill_long %>%
    #   ggplot(aes(x = tvec, y = molecules, color = simulation)) +
    #   geom_point() +
    #   geom_line() +
    #   geom_line(data = gill_summ,
    #             aes(x = tvec, y = mean_molcs), color = "red", size = 2,
    #             inherit.aes = FALSE) +
    #   geom_point(data = gill_summ, aes(x = tvec, y = mean_molcs),
    #              color = "red", size = 2,
    #              inherit.aes = FALSE) +
    #   guides(color = "none") +
    #   labs(x = "Time", y = "Simulation Value")
    quantile_vals <- c(0.025, 0.5, 0.975)
    gill_qntl <- gill_long %>%
      group_by(tvec) %>%
      summarize(
        q_val = quantile(molecules,
                         # x is the column to compute the quantiles
                         probs = quantile_vals
        ),
        q_name = quantile_vals
      ) %>%
      pivot_wider(names_from = "q_name", values_from = "q_val",
                        names_glue = "q{q_name}")
    gill_qntl %>%
      ggplot() +
      geom_line(aes(x = tvec, y = q0.5),
                color = "red", size = 2, inherit.aes = FALSE
      ) +
      geom_ribbon(aes(x = tvec, ymin = q0.025, ymax = q0.975),
                  alpha = 0.2, fill = "red", inherit.aes = FALSE
      ) +
      guides(color = "none") +
      labs(x = "Time", y = "Ensemble average") +
      title(main=input$title)
  #   plot(unlist(BM),type="n",xlim=c(1,max(sapply(BM,length))),xlab="# Reactions",ylab=input$spec, main=title_r())
  #   lines(1:length(colMeans(do.call(rbind,lapply(BM, `length<-`, min(lengths(BM)))))),colMeans(do.call(rbind,lapply(BM, `length<-`, min(lengths(BM))))),type="o",col="red")
  # rect(par("usr")[1], par("usr")[3],
  #        par("usr")[2], par("usr")[4],
  #        col = "#ebebeb")
  #   # Grid blanco
  #   grid(nx = NULL, ny = NULL,
  #        col = "white", lwd = 2)
  #   mapply(lines,BM,col=seq_along(BM),lty=2)
    }
    )

     output$histo = renderPlot({pt2()})
}

shinyApp(ui = ui, server = server)

6.4 Apéndice D: Material complementario

Los archivos de código y la Guía para la parte usuaria para la aplicación se pueden descargar desde:

https://github.com/jp-rgb/randomverse

La aplicación también se puede ver y usar en línea en la actualidad. Sin embargo, el sitio que los aloja tiene “horas activas” mensuales limitadas bajo el plan de alojamiento gratuito, por lo que no se garantiza que la aplicación esté accesible allí en todo momento. La liga para las aplicaciones en línea se proporcionan a continuación:

https://randomverse.shinyapps.io/RV_int/

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

7 Referencias

CAMCM. 2015. Historia. https://camcm.aefcm.gob.mx/wordpress/index.php/historia/.

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.

CUAED-UNAM. 2020. “Curso de Avances y Desarrollos En Biomatemáticas Para Madems Del Suayed.” Issuu. https://issuu.com/pedrodaniellaramaldonado/docs/avances_y_desarrollos_en_biomatem.

García Palacios, Ernesto. 2015. “La Educación a Distancia En México: Una Nueva Realidad Universitaria (Judith Zubieta y Claudio Rama).” Revista Mexicana de Bachillerato a Distancia 7 (14): 2. https://doi.org/10.22201/cuaed.20074751e.2015.14.65264.

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.

Maciel, Francisco Javier Chávez, and Ma. Trinidad Murguía Ángeles. 1970. “La Educación Media Superior a Distancia En México y Sus Efectos Para La Equidad Educativa.” Redalyc.org. Universidad de Guadalajara. https://www.redalyc.org/articulo.oa?id=68813176003.

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

“Sitio Web Oficial Del CBTIS No. 179-a.” n.d. https://www.cbtis179.edu.mx/portal/index.php/informacion-general/oferta-educativa/opcion-educativa-autoplaneada.

“Sitio Web Oficial Del CBTIS No. 179-b.” n.d. https://www.cbtis179.edu.mx/portal/index.php/de-interes/riems.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wilkinson, Darren James. 2020. Stochastic Modelling for Systems Biology. Chapman &amp; Hall/CRC.