Utilidad de las restricciones que reflejan la estabilidad del sistema en los análisis de modelos biológicos

Citación: Kariya Y, Honma M, Tokuda K, Konagaya A, Suzuki H (2022) Utilidad de las restricciones que reflejan la estabilidad del sistema en los análisis de modelos biológicos. PLoS Comput Biol 18(9): e1010441. https://doi.org/10.1371/journal.pcbi.1010441

Editor: James Gallo, Universidad de Buffalo – Universidad Estatal de Nueva York, ESTADOS UNIDOS

Recibió: 4 de diciembre de 2021; Aceptado: 26 de julio de 2022; Publicado: 9 de septiembre de 2022

Derechos de autor: © 2022 Kariya et al. Este es un artículo de acceso abierto distribuido bajo los términos de la Licencia de atribución de Creative Commonsque permite el uso, la distribución y la reproducción sin restricciones en cualquier medio, siempre que se acredite el autor original y la fuente.

Disponibilidad de datos: Todos los datos relevantes se encuentran dentro del manuscrito y sus archivos de información de respaldo.

Fondos: HS recibió JSPS KAKENHI Grant Número 18H04162, 24249034 y 22136001 de la Sociedad Japonesa para la Promoción de la Ciencia (https://www.jsps.go.jp/english/index.html). Los patrocinadores no tuvieron ningún papel en el diseño del estudio, la recopilación y el análisis de datos, la decisión de publicar o la preparación del manuscrito.

Conflicto de intereses: Los autores han declarado que no existen intereses contrapuestos.

Contenidos ocultar

Introducción

Comprender y predecir los comportamientos dinámicos de un sistema biológico complejo, como un sistema metabólico. [1–3] o una señal [4,5], en silico la simulación con un modelo matemático que consta de múltiples ecuaciones diferenciales ordinarias (EDO) puede ser útil. Sin embargo, la utilidad de este enfoque es limitada debido a las dificultades para determinar los valores de los parámetros cinéticos. [6]. En algunos casos, las mediciones experimentales pueden estimar los valores de los parámetros directamente [7], aunque debemos ser cautos a la hora de adoptar los valores medidos. Informes anteriores han indicado que los valores de los parámetros cinéticos determinados en vitro usando enzimas recombinantes purificadas no están en buen acuerdo con los valores estimados bajo en vivo condiciones que utilizan espectrometría de correlación cruzada de fluorescencia, en parte debido a los efectos de hacinamiento molecular en los compartimentos intracelulares [8,9]. Estos hallazgos también indican el riesgo de reutilizar los valores de los parámetros cinéticos de un modelo existente que describe la misma reacción en diferentes condiciones, como un orgánulo, una célula o un tejido diferentes. En ciertos casos, la medición de los valores de los parámetros que describen el comportamiento molecular microscópico es bastante difícil. [8]. Cuando el modelo incluye demasiados parámetros cinéticos, no es práctico estimar todos los valores de los parámetros experimentalmente. Otra forma de determinar los valores de los parámetros cinéticos es ajustando el modelo a los datos observados de las variables del modelo. [5]. Actualmente, varios algoritmos de ajuste están disponibles en casos en los que la flexibilidad de un modelo puede verse restringida por una cantidad suficientemente grande de puntos de datos y rangos de exploración suficientemente estrechos de valores de parámetros. [10,11]. Sin embargo, la biología molecular avanza continuamente y el número de componentes a incorporar en un modelo aumenta. Por lo tanto, se vuelve más difícil y lento determinar los valores exactos de todos los parámetros a medida que aumenta el tamaño del modelo, se expande el espacio de parámetros permitido y aumenta el número requerido de puntos de datos. La dificultad para determinar los parámetros está directamente relacionada con la fiabilidad de los resultados de la simulación. El amplio espacio de parámetros dará como resultado el riesgo de sobreajustar el modelo a los datos y perder otros valores de parámetros posibles que también explican los resultados experimentales, lo que lleva a una mala interpretación del comportamiento del sistema. Sin embargo, una simulación basada en ODE es necesaria en tales casos, ya que las observaciones experimentales de numerosas variables suelen requerir mucho tiempo y mano de obra.

Con base en estas observaciones, nuestra motivación en el presente estudio es obtener la mayor cantidad de información posible sobre la dinámica del modelo, incluso en los casos en que los valores de los parámetros apenas se determinan. Proponemos un proceso alternativo de búsqueda de parámetros adoptando las siguientes dos ideas. Primero, decidimos centrarnos en las características generales compartidas por muchos modelos biológicos. Las concentraciones de varios metabolitos en cada tejido se mantienen en ausencia de un estímulo extrínseco, como la ingesta de alimentos, pero cambian en respuesta al estímulo y luego se restauran nuevamente a los niveles basales después de un cierto período de tiempo en ausencia de el estimulo [12]. De manera similar, una vía de señalización intracelular se activa en respuesta a la estimulación del ligando y luego se restaura al estado inactivo después de la eliminación del estimulador. [13]. Estas propiedades de estabilidad y resiliencia se observan en muchos sistemas biológicos en estado homeostático. [14,15]. Una de nuestras ideas clave es utilizar estas propiedades en el proceso de búsqueda de parámetros. Si se puede suponer que un sistema de interés es biológicamente estable y resistente (BSR), el tamaño del espacio de parámetros que se enfocará se reduce significativamente para satisfacer esta restricción. Aunque hay excepciones, la característica BSR es ampliamente compartida por varios sistemas biológicos y esta restricción es aplicable para resolver muchos problemas (Fig. 1A). En segundo lugar, decidimos evitar la determinación exacta del conjunto de parámetros cinéticos óptimos, una serie de combinaciones de valores de parámetros, y más bien apuntamos a determinar el subconjunto completo en el espacio de parámetros donde el modelo satisface la restricción BSR. A menudo ocurre que la cantidad de datos observados es insuficiente para ajustarse al modelo. Además, algunos parámetros cinéticos no son lo suficientemente sensibles a las variables observables debido a la estructura del sistema, y ​​los valores de los parámetros no se pueden determinar con una precisión de estimación aceptable. [16]. Por lo tanto, adoptamos un enfoque para recopilar tantos conjuntos de parámetros como sea posible, satisfaciendo así las condiciones de las salidas del sistema y las restricciones de BSR. Con este propósito, para permitir el manejo de modelos más complicados, modificamos el método global Cluster-Newton (CNM), que fue desarrollado para obtener múltiples conjuntos de parámetros que satisfacen algunas condiciones de todo el espacio de parámetros sin un alto costo computacional. [17]. En el presente estudio, evaluamos la capacidad de nuestro método propuesto para estimar el espacio de parámetros permitido de un modelo cinético y también evaluamos la utilidad de los conjuntos de parámetros adquiridos para extraer algunas características de un sistema de interés. El método propuesto contribuirá a una mejor comprensión de las propiedades de un sistema biológico con una estructura de red conocida pero valores de parámetros cinéticos desconocidos.

miniatura

Figura 1. Modelo de sistema biológico en estado homeostático.

(A) Muchos sistemas biológicos se comportan en torno a un cierto estado estable. Tal comportamiento estable se puede observar en muchos sistemas biológicos, aunque el tamaño de la región estable alrededor del estado estable difiere según los sistemas. (B) Estructura asumida de la región estable biológica. Si hay regiones convexas hacia arriba en la región enfocada potencial, el tamaño de la cuenca de atracción es limitado, por lo tanto, nuestro concepto es disminuir dichas regiones en la búsqueda de parámetros. (C) Se consideran tres factores para buscar los conjuntos de parámetros que satisfacen la condición BSR y se incorporan en la función objetivo para minimizar.

https://doi.org/10.1371/journal.pcbi.1010441.g001

Resultados

Formalización de la propiedad BSR de un sistema biológico

En esta sección, definimos el sistema BSR como un sistema dinámico que satisface las siguientes dos condiciones: primero, el sistema tiene un punto fijo estable y opera alrededor de él y, segundo, el sistema regresa cerca del punto fijo dentro de un tiempo específico. después de una perturbación externa moderada a las variables de estado (Figura 1A). Consideramos un sistema que satisface la siguiente ecuación diferencial ordinaria:
(1)

dónde representa el estado del sistema, como las concentraciones de metabolitos, y es el parámetro. A continuación, formulamos las dos condiciones mencionadas anteriormente del sistema BSR.

Existe tal que
(2)
(3)

donde ‖ ‖ denota norma, X* es un punto fijo, rmáximo es la distancia que define el rango tolerable (tamaño de un área estable alrededor del punto fijo), y ε es la distancia entre el punto fijo y el estado del sistema, con lo cual se supone que el estado está cerca del punto fijo. Más lejos, es el tiempo mínimo para volver más cerca de ε al punto fijo, a~b indica a y b tienen el mismo orden de magnitud, y Δtconversión especifica la escala de tiempo del sistema para volver significativamente cerca del punto fijo.

Diseño de la función objetivo para lograr la condición BSR

A continuación, diseñamos nuestra función objetivo en el espacio de parámetros para identificar los conjuntos de parámetros con los que el sistema cumplió la condición BSR definida anteriormente.

Primero, consideramos la condición en la que el sistema tenía un punto fijo en X* (ecuación 2). Encontrar la raíz de la siguiente función objetivo produce la condición establecida en la ecuación 2:
(4)

dónde ‖ ‖2 denota norma L2.

A continuación, consideramos las condiciones establecidas en Eq 3. Esta condición tiene dos requisitos: (i) las órbitas que comienzan en las condiciones iniciales en la región especificada deben converger para X*, y (ii) las órbitas deben converger a un punto fijo dentro de una escala de tiempo específica. Menck y otros. conceptualizaron el punto anterior como la “estabilidad de la cuenca”, y propusieron cuantificarlo calculando la relación de las condiciones iniciales que convergen al atractor [18]. En principio, se pueden buscar parámetros que satisfagan estas dos condiciones si se conoce la dependencia de la relación de las órbitas convergentes y el tiempo de convergencia de los parámetros. Sin embargo, la duración de diferentes condiciones iniciales en el espacio de fase de un sistema no lineal de alta dimensión es prácticamente difícil de calcular sin resolver la ecuación diferencial a través de simulaciones numéricas. La incorporación de múltiples simulaciones numéricas para diferentes condiciones iniciales en cada paso de iteración en el proceso de optimización de parámetros implica un alto costo computacional. Por lo tanto, en lugar de evaluar directamente la dependencia de los parámetros de la estabilidad de la cuenca y el tiempo de convergencia, utilizamos un enfoque diferente: construir funciones objetivo que pudieran calcularse sin simulaciones numéricas del sistema. Implementamos la siguiente función objetivo cuya optimización se espera que produzca una expansión de la estabilidad de la cuenca en la región B:
(5)
(6)
(7)

dónde DF(X) es la matriz jacobiana de F(X) a X, METRO(X) es la suma de DF(X) y su transpuesta, B es la región de estabilidad predeterminada, λ(METRO(X)) es el valor propio de la matriz METRO(X), y ReLU(x)≝max(0, x). Más lejos, X(k) es un conjunto finito de puntos elegidos aleatoriamente de B en el número de iteración ka los que en adelante nos referiremos como “puntos de observación”, y norteobservación es el número de puntos de observación. Tenga en cuenta que todos los valores propios de la matriz METRO(X) son números reales porque METRO(X) es simétrico. Los puntos de observación X(k) fueron muestreados repetidamente durante la optimización (Figura S1). La condición max ReLU(λ(METRO(X))≤0 indica que la distancia entre las órbitas vecinas en el espacio de fase no se expande en X. Suponemos que es bastante natural esperar que la estabilidad de la cuenca de un punto fijo dado aumente cuando la contracción local de las órbitas vecinas está asegurada en un gran número de puntos que rodean los puntos fijos, dado que F(X, m) se describe mediante una función paramétrica suave con parámetros . En realidad, podemos probar que si el valor propio máximo de METRO(X) es menor que 0 en un disco que rodea un punto fijo, cualquier órbita que comience en la condición inicial en este disco converge exponencialmente a este punto fijo (ver Información S1 para la demostración).

Finalmente, diseñamos una función para reflejar el intervalo de tiempo requerido para la convergencia del estado del sistema hacia un punto fijo. Nos enfocamos en la escala de tiempo característica de relajación alrededor del punto fijo, que está determinada por la mayor parte real de los valores propios de la matriz jacobiana en el punto fijo. Debido a que nos enfocamos en los sistemas biológicos que tienen un comportamiento suave alrededor del punto fijo, esperamos que la escala de tiempo de relajación alrededor de la región enfocada sea similar a la del punto fijo. Por lo tanto, dada la escala de tiempo de relajación deseada del sistema Δtconversión (Eq 3), nuestra idea es incluir la relajación más lenta en el punto fijo en la función objetivo de la siguiente manera:
(8)

donde max Re(λ*) es la mayor parte real de los valores propios de la matriz jacobiana de F(X) a X*, y λ*objetivo = −1/Δtconversión es un valor constante negativo que define el estándar de la velocidad de relajación más lenta. En los experimentos biológicos, los valores medidos a menudo varían dentro de un cierto rango a pesar de realizar mediciones experimentales para múltiples conjuntos en la misma condición estable. Esto puede deberse a errores de medición; sin embargo, la existencia de múltiples estados estables puede formar una región neutralmente estable alrededor del valor medido. Por lo tanto, permitimos la existencia de valor(es) propio(s) cero e implementamos la siguiente función objetivo utilizando valores propios distintos de cero:
(9)

Esta función objetivo permite optimizar a 0 las partes reales más grandes de los valores propios. Por lo tanto, el cambio de estado del sistema en ciertas direcciones no provoca que regrese hacia el punto fijo. Incluso bajo esta configuración, el sistema seguiría siendo estable según Lyapunov, que a menudo se considera un indicador de estabilidad biológica. [19,20].

Finalmente, la suma ponderada o el vector concatenado de Oarreglar(m), y Orelax(m) se utilizaron como funciones de objetivo único para optimizar y producir conjuntos de parámetros que satisfacen la condición BSR. Considerando que los valores de la función objetivo de Oarreglar(m), y Orelax(m) eran términos diferentes matemáticamente, el peso necesita ajuste.

El enfoque para revelar el espacio de parámetros que satisface BSR y sus ventajas

La función objetivo diseñada es útil para encontrar conjuntos de parámetros que satisfagan BSR. Si podemos encontrar un espacio de parámetros que satisfaga BSR por completo, hay algunas ventajas en los análisis del modelo dinámico enfocado (Figura 2A). Primero, se espera que los conjuntos de parámetros reales para el modelo se incluyan en el espacio BSR. Esta función es potencialmente útil en el proceso de determinación de parámetros, especialmente para modelos a gran escala en los que se incluyen muchos parámetros. Para esos modelos, el volumen del espacio de parámetros que se busca es enorme, mientras que las restricciones de BSR pueden reducir el volumen y mejorar la eficiencia de la búsqueda de parámetros. Villaverde y Banga han demostrado la importancia de las restricciones de la estructura del modelo (red) en el ajuste efectivo de los parámetros [21], lo que respalda el hecho de que reducir el espacio de parámetros de enfoque por restricciones del sistema biológico es útil en el proceso de determinación de parámetros. Tan también ha demostrado la ventaja de reducir el espacio de parámetros permitido. y otros. que se centran en los flujos de estado estacionario como restricciones en los parámetros [22]. En segundo lugar, incluso si no podemos determinar un conjunto de parámetros único para el modelo, BSR puede proporcionar algunas ideas para el modelo. El análisis de la distribución de conjuntos de parámetros que satisfacen BSR proporcionará la conectividad entre parámetros, como la correlación entre parámetros, lo que ayudará a comprender el modelo. Además, realizar múltiples simulaciones utilizando conjuntos de parámetros que satisfacen BSR bajo algunos tipos de perturbaciones puede revelar algunos comportamientos dinámicos comunes, que pueden interpretarse como una propiedad de los sistemas bajo BSR. Teniendo en cuenta estas ventajas potenciales, desarrollamos un algoritmo para encontrar tantos conjuntos de parámetros que satisfagan BSR como sea posible, basado en el método de newton de clúster global (gCNM) porque el método gCNM es un método utilizado para revelar todo el espacio de parámetros que satisface algunas restricciones (Fig. 2B) . Mientras que el gCNM original usa el método de Broyden como algoritmo base para la globalización y paso de optimización, optamos por utilizar el método L-BFGS para reducir el consumo de memoria y permitir una aplicación al modelo, incluidos muchos parámetros. Teniendo en cuenta que recopilar tantos conjuntos de parámetros como sea posible es clave para las ventajas que se muestran arriba, también incluimos otra modificación. Cada conjunto de parámetros encontrado por gCNM es agregado por ruido para expandir la distribución formada por todos los conjuntos de parámetros encontrados; a partir de entonces, se vuelve a someter al paso de optimización y globalización basado en L-BFGS (denominado g-LBFGS). El ciclo repetitivo del proceso de adición y optimización de ruido es una gran modificación del gCNM original (Figs. 2B y S2). Este gCNM modificado está diseñado para iterar hasta satisfacer las dos condiciones siguientes: (1) los valores medianos de los conjuntos de parámetros obtenidos para todos los parámetros no se modifican significativamente por la última iteración de gCNM, y (2) la mayoría de los valores de los parámetros encontrados en el La última iteración de gCNM está dentro del rango de distribución del conjunto de parámetros ya encontrado. Usando este algoritmo junto con las funciones objetivo BSR definidas anteriormente, esperábamos que una gran parte del espacio de parámetros satisfaga BSR y llamamos a este esquema completo la exploración exhaustiva del espacio de parámetros permitido para sistemas biológicos (TEAPS).

miniatura

Figura 2. Ventajas de enfocar BSR y marco de TEAPS.

(A) Se representaron las ventajas de enfocar BSR. Ver detalle en el texto. (B) El marco de exploración exhaustiva de conjuntos de parámetros que satisfacen las restricciones BSR, que llamamos TEAPS.

https://doi.org/10.1371/journal.pcbi.1010441.g002

Evaluación del rendimiento del algoritmo utilizando los modelos de muestra más simples.

A continuación, evaluamos la cobertura de los múltiples conjuntos de parámetros obtenidos usando TEAPS en todo el espacio de parámetros que satisface BSR. Primero, se prepararon modelos de muestra (es decir, T1, T2 y T3) que constan de tres entidades biológicas que reflejan las estructuras del sistema biológico existente (Figura 3A). T1 representa el sistema de retroalimentación positiva y T2 representa el sistema de retroalimentación negativa. Estos sistemas se observan ampliamente en varios procesos biológicos, incluida la señalización. [23,24] y vías metabólicas [25]. Estos modelos no contienen reacciones reversibles que son fundamentales para varios procesos biológicos. [26]. Por lo tanto, preparamos T3 como una reacción reversible con una entidad reguladora. Establecimos las condiciones BSR de estos modelos de la siguiente manera: el punto fijo del sistema se estableció en X* = (1,1,⋯,1)T porque el valor absoluto de la cantidad de una entidad en el sistema depende de la unidad utilizada y la información importante es la alteración relativa del valor. De manera similar, el valor absoluto del tiempo depende de la unidad utilizada. Por lo tanto, el tiempo de relajación hasta el punto fijo después de la perturbación se ajustó configurando λ*objetivo = −0.3 para observar la tendencia de relajación hasta el momento t = 100. Los puntos de observación para confirmar el área estable alrededor de los puntos fijos se generaron aleatoriamente como Xmetro,norte = 1+adónde metro es un índice para los puntos de observación y norte es un índice para el espacio de fase, a fue tomado de una distribución uniforme sobre [−d, d]. El valor d = 0.1 se usó en el paso final del algoritmo porque es probable que un sistema biológico sea estable frente a una perturbación que altere la cantidad de cierta entidad en menos del 10 % y regrese al punto basal después de eliminar la perturbación en mayoria de los casos [27,28]. Para calcular eficientemente en los casos de modelos complicados, valor d está previsto que aumente gradualmente. Es decir, en la etapa CNM y la primera subrutina de la etapa de globalización, el valor se establece en 0,004, seguido de 0,036 y 0,1 en la segunda y tercera subrutina de la etapa de globalización, respectivamente. Usando estos hiperparámetros, obtuvimos múltiples conjuntos de parámetros para cada modelo de muestra para satisfacer la condición BSR. Primero, evaluamos si el algoritmo funciona, como se esperaba, para asegurar un área estable alrededor del punto fijo. Para monitorear la estabilidad del sistema alrededor del punto fijo, se calculó la estabilidad de la cuenca para los conjuntos de parámetros adquiridos después de cada paso de cálculo. Los resultados del cálculo para cada modelo mostraron que la estabilidad estimada de la cuenca alrededor del punto fijo Xss aumentó a medida que avanzaba el paso de cálculo (Fig. 3B). El valor aumentó en gran medida en la etapa de globalización, lo que era consistente con la estructura del algoritmo; es decir, los conjuntos de parámetros candidatos que satisfacían vagamente la condición BSR se recolectaron en la etapa CNM, y se llevaron a cabo ciclos de optimización adicionales para capturar conjuntos de parámetros que satisfacían estrictamente la condición BSR en la etapa de globalización. Luego comparamos los conjuntos de parámetros obtenidos usando este enfoque con los obtenidos usando la búsqueda de fuerza bruta (Fig. 3C–3E). Por búsqueda de fuerza bruta, cada modelo mostró una amplia distribución para algunos parámetros, como v2 y kilómetros2 para T1. Por el contrario, hubo casos de convergencia a una distribución estrecha de menos de 0,3, como k1, k3, y k4 para T1. Dichos perfiles fueron reproducibles usando nuestro enfoque para cada modelo cuando se compararon los histogramas de las distribuciones de los valores de los parámetros. Debido a que la simple comparación de los histogramas no puede delinear claramente las formas de distribución en el espacio de múltiples dimensiones, también evaluamos la densidad en el espacio de múltiples dimensiones mediante el suavizado del núcleo multivariado y las normalizamos por la suma total. Luego comparamos las formas (Fig. 3C–3E, paneles medio e inferior) en función de la densidad. La densidad se calculó en los puntos de la cuadrícula del espacio de búsqueda de parámetros, y cada punto de la cuadrícula se indexa en un orden preestablecido. Por lo tanto, las densidades en los puntos de la cuadrícula se representaron como dentsitogramas. La similitud de los perfiles de densidad normalizados se confirmó utilizando el valor de similitud del coseno y el valor de divergencia de Jensen-Shannon (JS) entre los vectores linealizados de cada tensor de densidad (Tabla 1). Como referencia negativa, generamos vectores aleatorios que se muestrearon para dar las mismas frecuencias de densidad que el vector linealizado para la búsqueda de fuerza bruta. Los índices entre las distribuciones de parámetros obtenidas por nuestro enfoque y la búsqueda de fuerza bruta superan claramente a aquellos entre el vector aleatorio y la búsqueda de fuerza bruta. Estos resultados indican que la estimación del espacio de parámetros que satisface BSR usando TEAPS es práctica para modelos biológicos simples.

miniatura

Fig. 3. La búsqueda de TEAPS expande la estabilidad de la cuenca durante su proceso y encontró los conjuntos de parámetros que satisfacen la condición BSR.

(A) Se utilizaron tres modelos básicos para la implementación de la BSR. (B) La estabilidad estimada de la cuenca se incrementó mediante el proceso de optimización individual (gL-BFGS) de conjuntos de parámetros de búsqueda que satisfacen la condición BSR utilizando el algoritmo de búsqueda TEAPS. (C–E) Comparación de las distribuciones de parámetros encontradas por los TEAPS y las búsquedas de fuerza bruta como histogramas o perfiles de densidad multidimensional para los modelos T1 (C), T2 (D) y T3 (E). Los conjuntos de parámetros que satisfacían la condición BSR se buscaron mediante TEAPS (arriba a la izquierda) o mediante la búsqueda de fuerza bruta (arriba a la derecha), y se compararon los histogramas de distribución de los conjuntos de parámetros encontrados. Los perfiles de densidad de los conjuntos de parámetros obtenidos se compararon siguiendo la estimación de la densidad del kernel (centro a la izquierda: TEAPS, centro a la derecha: búsqueda de fuerza bruta). Se superpusieron dos perfiles de densidad (abajo).

https://doi.org/10.1371/journal.pcbi.1010441.g003

Evaluación del desempeño de TEAPS utilizando varios modelos de estructura simple

A continuación, evaluamos el rendimiento de TEAPS utilizando un modelo de tres variables ligeramente complicado que refleja variaciones en sistemas biológicos distintos de T1-3 (Figura 4A). Cambio de la vía activada, que a menudo se observa en las vías de señalización. [29], se puede simplificar como se expresa en T4, es decir, la activación de una vía suprime la activación de otra vía y viceversa. T5 es un modelo en el que la regulación inhibitoria en T4 se reemplaza por una regulación positiva. En este caso, las dos vías no se excluyen mutuamente y cooperan para evitar la sobreacumulación de productos finales. La distribución de los conjuntos de parámetros adquiridos por TEAPS se alineó con los obtenidos por la búsqueda de fuerza bruta (Fig. 4B y 4C, Tabla 1). Estos resultados indican una gran eficiencia de TEAPS en la búsqueda del espacio de parámetros permitido para varios modelos de tres variables. Posteriormente, preparamos modelos de cuatro variables y evaluamos el desempeño de TEAPS. El componente mínimo de la cascada de señalización podría modelarse como se expresa en T6 (Fig. 4A). El algoritmo podría encontrar una gran parte del espacio de parámetros permitido (Fig. 4D). Regulación de retroalimentación positiva, que a veces se incluye en cascadas de señalización. [30]fue añadido a T6 para construir T7 (Fig. 4A). Nuevamente, el algoritmo encontró una gran parte del espacio de parámetros permitido (Fig. 4E). Finalmente, examinamos T8, en el que una reacción particular está regulada por múltiples factores. Esta estructura se informó anteriormente como una red representativa en sistemas biológicos. [31]. Este modelo incluye nueve parámetros, lo que hace que la cantidad de puntos de cuadrícula sea enorme en la búsqueda de fuerza bruta; por lo tanto, fijamos un valor de k0 y comparó los espacios para los otros ocho parámetros. Incluso en presencia de leyes cinéticas complejas por múltiples regulaciones en una reacción particular, TEAPS estimó con éxito el espacio de parámetros que satisface BSR (Fig. 4F), que se entiende intuitivamente mediante la combinación de gráficos tridimensionales para tres de los ocho parámetros (S3 Fig. ). Teniendo en cuenta la posibilidad de un muestreo sesgado por la búsqueda basada en puntos de cuadrícula, también realizamos una búsqueda aleatoria basada en vectores para revelar la distribución de parámetros que satisfacen BSR para todos los modelos de muestra. Como resultado, confirmamos los tres enfoques, la búsqueda basada en puntos de cuadrícula, la búsqueda basada en vectores aleatorios y TEAPS, encontramos una distribución similar para todos los modelos, que respaldan la capacidad de búsqueda de TEAPS (Fig. S4). A continuación, evaluamos el tiempo de cálculo para todo el espacio permitido de búsqueda entre TEAPS y fuerza bruta para el modelo T8, debido a que el número de parámetros es 8, que es el más grande en los modelos de muestra. TEAPS tarda aproximadamente 100 min en completarse, mientras que la búsqueda de fuerza bruta tarda unas 16 h en el mismo entorno computacional. Considerando la fuerza bruta buscada para un espacio formado solo por cinco parámetros ya que las restricciones de punto fijo fueron dadas analíticamente de antemano, se sugirió el muestreo efectivo por TEAPS. Además, realizamos simulaciones utilizando conjuntos de parámetros de TEAPS para confirmar que el sistema regresa cerca del punto fijo después de agregar la perturbación a las variables de estado (informes computacionales en S3 Information). Una serie de pruebas confirmaron el desempeño de TEAPS en la exploración del espacio de parámetros permitido para los sistemas BSR.

miniatura

Figura 4. Evaluación del rendimiento de TEAPS utilizando varios modelos de estructura simple.

(A) Se prepararán cinco modelos de muestra adicionales para la prueba de rendimiento. (B–F) Comparación de los resultados de TEAPS y los resultados de búsqueda de fuerza bruta para los modelos T4 (B), T5 (C), T6 (D), T7 (E) y T8 (F). Los conjuntos de parámetros que satisfacían la condición BSR se buscaron mediante TEAPS (arriba a la izquierda) o mediante la búsqueda de fuerza bruta (arriba a la derecha), y se compararon los histogramas de distribución de los conjuntos de parámetros encontrados. Los perfiles de densidad de los conjuntos de parámetros encontrados se compararon siguiendo la estimación de la densidad del kernel (centro a la izquierda: TEAPS, centro a la derecha: búsqueda de fuerza bruta). Se superpusieron dos perfiles de densidad (abajo).

https://doi.org/10.1371/journal.pcbi.1010441.g004

Modelo de vía NF-κB

A continuación, examinamos la aplicabilidad de TEAPS en búsquedas de espacio de parámetros en un modelo de red biológica real. Se centró en el modelo de la vía de señalización de NF-κB. Este sistema incluye inhibición por retroalimentación contra NF-κB aguas abajo de su activación; por lo tanto, este sistema puede mostrar un comportamiento oscilatorio continuo [32–34] (Figura 5A). En observaciones reales, este sistema muestra una oscilación amortiguada una vez que se activa la vía de señalización, y este comportamiento es consistente con la restricción BSR. Para analizar matemáticamente este comportamiento característico, se han desarrollado varios modelos ODE basados ​​en varios parámetros, muchos de los cuales se determinaron ajustando las observaciones reales. [33,34] (Fig. 5B, línea roja). También analizamos la propiedad de estabilidad de la cuenca para el conjunto de parámetros informado. Los estados iniciales del sistema se generaron aleatoriamente alrededor del punto fijo; luego, se evaluó la probabilidad de convergencia alrededor del punto fijo. Por lo tanto, se confirmó que existía un cierto tamaño de estabilidad de la cuenca alrededor del punto fijo (Fig. 5C, actual). Por lo tanto, examinamos si TEAPS puede explorar el espacio de parámetros, incluidos los valores de parámetros informados. El punto fijo objetivo se fijó en el estado en el que convergía el modelo con el conjunto de parámetros notificados, y la relajación objetivo se fijó para que convergiera en aproximadamente 20 h para garantizar la coherencia con las observaciones reales. [33,34]. El tamaño objetivo de la cuenca se fijó para que fuera tan grande como el de los modelos de muestra. Confirmamos que el conjunto de parámetros informado es consistente con esta configuración para BSR (informes computacionales en información S4). Con esta configuración, finalmente obtuvimos 1095 conjuntos de parámetros por TEAPS dentro de las 2 h del tiempo de cálculo en nuestro entorno y confirmamos que la estabilidad estimada de la cuenca alrededor del punto fijo también se elevó en este caso después del g-LBFGS (Fig. 5D). Para muchos parámetros, los valores informados se incluyeron en el rango de valores de parámetros encontrados por TEAPS (Fig. 5E). Sin embargo, los valores de los parámetros encontrados están ampliamente distribuidos, lo que dificulta juzgar si las características de los parámetros críticos están realmente determinadas por TEAPS. Para evaluar aún más si el conjunto de parámetros informado está incluido en una variedad formada por los conjuntos de parámetros obtenidos teniendo en cuenta la conectividad entre parámetros, se realizó un análisis de componentes primarios (PCA) contra los conjuntos de parámetros obtenidos y se examinó la inclusión en los ejes de componentes primarios. El número de componentes se fijó para que fuera el mismo que las dimensiones del espacio de parámetros para identificar las direcciones más anchas y más estrechas. En todos los componentes primarios (PC), excepto PC 1, el valor informado estaba claramente dentro del rango del espacio de parámetros obtenido. En particular, vale la pena centrarse en los ejes de PC con distribuciones más estrechas (PC 40 y más). TEAPS recogió conjuntos de parámetros limitados cuando se proyectó a estos ejes; por lo tanto, se espera que estos ejes definan restricciones estrictas para BSR. Se confirma claramente que el conjunto de parámetros informado está incluido en el rango encontrado en aquellos ejes con restricciones estrictas (Fig. 5F y 5G). A continuación, examinamos si los conjuntos de parámetros que dan un comportamiento oscilatorio no amortiguado estaban realmente fuera del espacio de parámetros encontrado. Para este propósito, buscamos un conjunto de parámetros que cambiara ligeramente del conjunto de parámetros informado. Encontramos que el modelo genera una oscilación no amortiguada cuando los valores de dos parámetros que describen la velocidad de traducción de IκBα y A20 se cambian de 0.5 a 5.0 (Fig. 5B, línea verde). Para este conjunto de parámetros, evaluamos la estabilidad de la cuenca alrededor de la región de estabilidad enfocada, que estaba alrededor del punto fijo dado por el conjunto de parámetros original. Por lo tanto, se confirmó que la estabilidad de la cuenca desaparecía cuando se usaba el conjunto de parámetros oscilatorios (Fig. 5C, oscilación). Cuando este conjunto de parámetros oscilatorios se proyectó en los ejes de los componentes primarios, se confirmó que este conjunto de parámetros existía fuera de la distribución dominante en PC 47 (Fig. 5G), lo cual es consistente con el hecho de que el comportamiento oscilatorio no satisface la restricción BSR. Estas observaciones sugieren que TEAPS es aplicable para buscar conjuntos de parámetros que satisfagan el BSR para modelos biológicos reales.

miniatura

Figura 5. Aplicación de TEAPS para el modelo de vía de señalización de NF-κB.

(A) La estructura modelo de las vías de señalización de NF-κB. Esta vía de señalización incluye la inhibición por retroalimentación de la activación de NF-κB a través de IκBα y A20. (B) Cursos temporales de activación de NF-κB. El conjunto de parámetros que se ajusta a la observación in vitro genera una oscilación amortiguada que es consistente con la observación experimental (rojo), mientras que el modelo genera una oscilación continua cuando se cambia el parámetro relacionado con las inhibiciones de retroalimentación (verde). (C) La estabilidad estimada de la cuenca para el conjunto de parámetros informados (real) y el conjunto de parámetros que genera una oscilación continua (oscilación). (D) La estabilidad estimada de la cuenca aumentó mediante el proceso de optimización individual (gL-BFGS) cuando se aplicaron TEAPS al modelo NF-κB. (E) La distribución de conjuntos de parámetros encontrados por TEAPS para el modelo NF-κB. Los valores de los parámetros encontrados se distribuyeron ampliamente (azul). Para muchos parámetros, el valor informado anteriormente (rojo) se incluyó en el rango encontrado por TEAPS. (F, G) La distribución de los valores de los parámetros encontrados por TEAPS en el espacio de la PC. Los valores informados (rojo) se incluyeron en la distribución encontrada para la mayoría de los ejes, mientras que el conjunto de parámetros oscilatorios estaba fuera de la distribución dominante para las PC posteriores. (F) Se muestran las distribuciones para todas las PC. (G) Las distribuciones de PC posteriores se amplían.

https://doi.org/10.1371/journal.pcbi.1010441.g005

Modelo de vía metabólica del ácido araquidónico

Luego evaluamos la utilidad de los conjuntos de parámetros encontrados por TEAPS en un sistema biológico más complicado. Un modelo informado previamente con el cual se han analizado bien las propiedades dinámicas del sistema es adecuado para este propósito. Por lo tanto, seleccionamos una ruta metabólica del ácido araquidónico (AA) en células polimorfonucleares (PMN), células endoteliales (EC) y plaquetas (PLT). [1,35]. Estos modelos se usaron para simular los efectos de 11 medicamentos antiinflamatorios no esteroideos (AINE) (es decir, aspirina, ibuprofeno, naproxeno, ácido 6-metoxi naftaleno acético (6-MNA), paracetamol, indometacina, meloxicam, nimesulida, diclofenaco , celecoxib y rofecoxib) sobre las concentraciones de prostaglandinas. Combinamos los tres modelos separados de vías AA en PMN, EC y PLT en un solo modelo que consta de 75 entidades biológicas, 60 reacciones y 109 parámetros cinéticos (Fig. 6A). En este modelo, se adoptó un valor de parámetro único para la constante de Michaelis-Menten para describir la misma reacción enzimática en diferentes tipos de células. Las reacciones que incluyen tales parámetros cinéticos compartidos se colorean de negro, y las reacciones específicas del tipo de célula se colorean de manera diferente como en la Fig. 6A. En este sentido, los modelos de diferentes tipos de células se vincularon a un solo modelo. Aplicamos TEAPS a este modelo y evaluamos su capacidad de búsqueda. El valor de cada variable en el punto fijo se fijó en 1, y se supuso que era permisible una perturbación del 10 % alrededor del punto fijo, como en los análisis del modelo de muestra. Dentro de un período de tiempo realista (alrededor de 3 días), obtuvimos 1391 conjuntos de parámetros de 4000 conjuntos iniciales que convergieron cerca del punto fijo objetivo. Además, durante el proceso de TEAPS, se confirmó que aumentó la estabilidad de la cuenca (Fig. 6B). Estas observaciones sugieren que TEAPS se puede aplicar a una cierta escala de un modelo biológico real, similar a este ejemplo.

miniatura

Figura 6. Utilidad de TEAPS en los análisis del modelo de ruta del ácido araquidónico.

(A) El modelo que contiene vías de ácido araquidónico en células polimorfonucleares (PMN), células endoteliales (EC) y plaquetas (PLT) se reconstruyó de acuerdo con la estructura del modelo informado anteriormente. Las constantes de Michaelis-Menten se establecieron en los mismos valores si había las mismas reacciones en diferentes celdas. Las reacciones que solo aparecen en cada tipo de célula están coloreadas en cian (PMN), rosa (EC) y marrón (PLT). Las reacciones con flechas negras se comparten con al menos dos tipos de células. Para modificaciones en las reacciones, las líneas que terminan en un círculo indican un efecto de promoción en una reacción, mientras que las líneas que terminan en una barra indican un efecto inhibitorio. Los efectos modificativos por metabolitos pero no por enzimas se describen como líneas discontinuas. (B) La estabilidad estimada de la cuenca aumentó mediante el proceso de optimización individual (gL-BFGS) cuando se aplicaron TEAPS al modelo de vía del ácido araquidónico. (C) Los conjuntos de parámetros obtenidos para el modelo de vía del ácido araquidónico se sometieron al análisis de componentes principales. Para determinar la forma completa de la distribución de parámetros, la cantidad de componentes se estableció en la cantidad de parámetros en el modelo. (D) Se buscaron las reacciones con una regulación estricta y robusta. Cada flecha corresponde a cada parámetro en cada reacción. Las reacciones relacionadas con los parámetros que aparecían como alta contribución a los ejes PCA en los que las distribuciones de parámetros son estrechas están coloreadas en rojo, mientras que las de cada uno de los ejes PCA en los que las distribuciones de parámetros son amplias están coloreadas en azul.

https://doi.org/10.1371/journal.pcbi.1010441.g006

A continuación, analizamos las propiedades de los conjuntos de parámetros obtenidos. Debido a que una de las características de TEAPS es la recopilación de múltiples conjuntos de parámetros del espacio de parámetros permitido, la forma de este espacio podría incluir la información de las restricciones del modelo para satisfacer la condición BSR. Los análisis de sensibilidad generalmente se realizan repitiendo simulaciones para determinar los parámetros o reacciones que afectan de manera crucial el comportamiento del sistema. Por el contrario, en los casos en que tenemos información sobre el espacio de parámetros permitido para satisfacer las condiciones de BSR, la identificación de la dirección más estrecha en el espacio de parámetros posiblemente conducirá a la identificación de parámetros cruciales que afectan el comportamiento del sistema. Dado que una simple evaluación de la distribución del valor de cada parámetro sería inapropiada debido a la conectividad entre parámetros, adoptamos PCA. El número de componentes se fijó para que fuera el mismo que las dimensiones del espacio de parámetros para identificar las direcciones más anchas y más estrechas (Figura 6C). También calculamos un índice denominado índice del eje central principal (PCI) para cada parámetro (consulte la sección Método); es un promedio ponderado calculado multiplicando cada número de eje de PC por la carga relativa de cada parámetro en su eje (Tabla S1). Este índice está diseñado para representar un número de eje de PC hipotético entre 1 y 109, donde cada parámetro contribuye en gran medida. Por lo tanto, se supone que los parámetros con un valor de PCI pequeño contribuyen a los ejes de PC en los que la distribución de parámetros encontrada es amplia. En este contexto, buscamos parámetros con valores de PCI pequeños, que se espera que varíen significativamente incluso si el sistema cumple con los requisitos de BSR. Estos parámetros se muestran en azul en el modelo (Fig. 6D). Tenga en cuenta que un parámetro, pero no una reacción, se indica con una flecha en la figura 6D. Entre los parámetros de color azul, los relacionados con la producción de AA en cada celda fueron distintivos. Dado que la producción de AA corresponde al proceso de iniciación de las reacciones metabólicas en cada célula, el sistema puede alcanzar un estado estacionario equilibrando los flujos de las reacciones aguas abajo con los flujos de producción de AA. Por lo tanto, es razonable que los parámetros relacionados con la producción de AA puedan tomar una amplia gama de valores para satisfacer las restricciones de BSR. Por el contrario, los parámetros que tienen un valor de PCI grande son los principales contribuyentes a los ejes de PCA en los que la distribución de parámetros encontrados es estrecha. En estos ejes PCA, la combinación lineal de varios valores de parámetros puede tomar un rango estrecho. Por lo tanto, los parámetros que contribuyen a esos ejes pueden tomar un rango de valores estrechamente limitado por otros parámetros que contribuyen al mismo eje PCA para cumplir con las restricciones de BSR, lo que lleva a la restricción de la flexibilidad de los parámetros. En la figura 6D, estos parámetros se muestran en rojo. Podemos notar fácilmente que las reacciones aguas abajo son menos flexibles que las otras reacciones. Estas reacciones son en gran medida procesos no enzimáticos cuyos flujos se describen mediante un parámetro, y el flujo de cada reacción debe equilibrarse con el flujo de la reacción subsiguiente en estado estacionario. Por lo tanto, es bastante razonable que estas reacciones se concentren en las direcciones más estrechas, lo que indica que los conjuntos de parámetros adquiridos por TEAPS parecen delinear con éxito el espacio de parámetros permitido que satisface la restricción BSR incluso en un sistema que consta de más de 100 parámetros cinéticos.

Aplicaciones de conjuntos de parámetros obtenidos por TEAPS en análisis biológicos de sistemas

En informes anteriores, se utilizó el modelo AA para comparar los efectos farmacológicos de los AINE utilizando parámetros ajustados a las observaciones experimentales. Por lo tanto, a continuación examinamos si se podían extraer características farmacológicas similares utilizando conjuntos de parámetros basados ​​en TEAPS. En informes anteriores, las acciones farmacológicas de los AINE se expresaron como factores de corrección dependientes de la concentración sobre parámetros cinéticos que describen las actividades de la ciclooxigenasa (COX)-1 y COX-2, y se analizaron las concentraciones de prostaglandinas en estado estacionario. [35]. Esto significa que los efectos farmacológicos se evaluaron como la alteración del punto fijo del sistema asociado con los cambios en los valores de los parámetros particulares. TEAPS fue diseñado para encontrar la condición donde el sistema se relajaba al mismo punto fijo cuando los valores iniciales de las variables de estado estaban a cierta distancia del punto fijo; por el contrario, el punto fijo posiblemente podría cambiar cuando se cambiaron los parámetros cinéticos particulares incluidos en cada conjunto de parámetros obtenidos por TEAPS. Por lo tanto, aunque TEAPS encuentra conjuntos de parámetros basados ​​en la estabilidad de los sistemas, los análisis farmacológicos como los del informe anterior se pueden realizar utilizando conjuntos de parámetros basados ​​en TEAPS. Como paso inicial, realizamos simulaciones utilizando cada conjunto de parámetros obtenido por TEAPS para identificar conjuntos de parámetros con los que se puede reproducir la reducción de la concentración de PGE2 por exposición a AINE, ya que la reducción de los niveles de PGE2 en PMN se utilizó como indicador del efecto farmacológico. de AINE en un informe anterior [35]. Examinamos si el nivel de PGE2 en PMN se puede reducir a alrededor del 10 % aumentando la concentración de cada fármaco, como se ha realizado en un informe anterior. Específicamente, la concentración inicial de cada fármaco se fijó en una centésima parte de la IC50 valor de cada fármaco frente a la COX-2, y la concentración del fármaco se incrementó en 1,5 veces por pasos hasta que el nivel de PGE2 simulado en los PMN se redujo a menos del 15 % del nivel sin exposición a los AINE. Los conjuntos de parámetros con los que no se logró la reducción de PGE2 incluso cuando la concentración del fármaco excedió cien veces la IC50 valor fueron excluidos de los siguientes análisis. En consecuencia, determinamos la concentración de cada fármaco para los 234 conjuntos de parámetros seleccionados (Fig. 7A). Posteriormente, se simuló la alteración en la relación entre la concentración de prostaciclina y la concentración de tromboxano A2 (relación PT) y se comparó con los valores informados. [35]dado que este índice suele utilizarse como indicador de reacciones adversas cardiacas [36,37]. Los valores medios de los cocientes de PT calculados en 234 conjuntos de parámetros disminuyeron en el siguiente orden: aspirina, ibuprofeno, naproxeno, 6-MNA, paracetamol, indometacina, meloxicam, nimesulida, diclofenaco, celecoxib y rofecoxib (Fig. 7B). Esto fue consistente con las observaciones en el informe original, en el que se estimó el valor de cada parámetro independientemente ajustando el modelo a las observaciones experimentales. A continuación, examinamos si era posible clasificar los medicamentos en función de los resultados simulados. Si las propiedades farmacológicas de los fármacos A y B son similares, es plausible suponer que las proporciones de PT para los fármacos A y B serán similares en cualquier simulación que utilice cada uno de los conjuntos de parámetros obtenidos. Por lo tanto, para cada fármaco, se construyeron vectores consistentes en proporciones de PT para todos los conjuntos de parámetros, y se evaluó la similitud de los vectores para todos los fármacos mediante análisis de agrupamiento (Fig. 7C). Un grupo consistió en ibuprofeno, naproxeno, 6-MNA y paracetamol. Estos medicamentos comparten una característica de baja selectividad de COX-2 [35]. Otro grupo consistió en fármacos de alta selectividad como meloxicam, nimesulida, diclofenaco, celecoxib y rofecoxib. Esta clasificación es consistente con las propiedades farmacológicas conocidas de los AINE. Estos resultados sugieren que los conjuntos de parámetros obtenidos usando TEAPS incluyen aquellos que reflejan algunos eventos biológicos farmacológicamente importantes.

miniatura

Figura 7. Parte de los conjuntos de parámetros obtenidos por TEAPS pueden reproducir la en vivo fenotipo asociado a la administración de AINE.

(A, B) Las situaciones en las que se administró cada AINE para suprimir los niveles de PGE2 en PMN hasta aproximadamente el 10 % del control se compararon con el artículo informado anteriormente. Cada punto azul corresponde a un resultado de un conjunto de parámetros obtenido por TEAPS. Los asteriscos rojos indican los valores informados por los resultados de la simulación utilizando el modelo informado anteriormente. Se muestran el nivel de PGE2 en PMN (A) y la proporción de PGI2 a TXA2 (relación de PT) (B), un marcador de rendimiento fisiológico. (C) El patrón de resultados farmacológicos mediante análisis de agrupamiento mostró la diferencia de las acciones farmacológicas. Se realizó un análisis de agrupamiento basado en el patrón de las proporciones de PT simuladas por los conjuntos de parámetros obtenidos bajo cada exposición al fármaco. La selectividad de COX-2 fue descifrada por el árbol filogenético.

https://doi.org/10.1371/journal.pcbi.1010441.g007

Discusión

En el presente estudio, intentamos obtener la mayor cantidad de información posible sobre la dinámica del modelo, incluso en los casos en que los valores de los parámetros apenas están determinados. Para analizar el comportamiento dinámico evitando el proceso de determinación de parámetros, un enfoque es la construcción de modelos basados ​​en datos, que a veces se logra mediante el aprendizaje automático. Por ejemplo, el enfoque de Henriques y otros. no requiere determinación de parámetros individuales [38]. Este enfoque es bastante poderoso cuando se dispone de datos de observación a gran escala a lo largo del tiempo. Cuando la atención se centra en los eventos intracelulares, como el sistema de transducción de señales, lamentablemente, en muchos casos, tales datos de series temporales para muchas moléculas biológicas no están disponibles. Creemos que también se requieren enfoques complementarios para tales casos para analizar los modelos biológicos. Propusimos un enfoque que evita la determinación precisa de parámetros y sin ningún requisito de datos de observación a gran escala. Para lograr esto, adoptamos un enfoque para recoger todos los conjuntos de parámetros permitidos bajo algunas restricciones biológicas.

Nos centramos en los sistemas biológicos que satisfacen las condiciones BSR en este estudio. El concepto de estabilidad dinámica de los modelos biológicos ha sido mencionado como robustez alrededor de 2004 [39]. En esta literatura, el autor explicó que muchos sistemas biológicos consisten en componentes de recepción de señales, componentes de procesamiento central y componentes de salida, y que en muchos casos hay bucles de retroalimentación desde los componentes de salida hasta los componentes de recepción de señales, lo cual es una estructura razonable para que los sistemas la mantengan. un cierto estado Además, el autor también mencionó que una robustez demasiado fuerte provoca la pérdida de una oportunidad en la transición de estado por cualquier estimulación; por lo tanto, se pensó que la robustez existía en un rango limitado alrededor de un cierto estado. Nuestro concepto de BSR refleja bien la robustez, por lo tanto, el método TEAPS será aplicable a la mayoría de los sistemas metabólicos. [1–3] y varios sistemas de señalización, como respuestas celulares a estímulos externos [4,5]. Para obtener conjuntos de parámetros que satisfagan este concepto de estabilidad, diseñamos funciones objetivo utilizando estrategias de aproximación lineal para el modelo, como calcular los valores propios de la matriz jacobiana alrededor del punto fijo. Este enfoque es razonable ya que la utilidad de los valores propios de la matriz jacobiana en la evaluación de la estabilidad del estado ha sido bien reconocida en el campo de la teoría de control. [40,41]. En nuestra función objetivo diseñada, necesitamos definir la región de interés para evaluar la estabilidad del estado. Establecimos tentativamente el tamaño del área estable para permitir un cambio del 10% en cada variable, según los dos informes anteriores. Un informe mostró que el nivel de expresión del ARNm tenía aproximadamente un 30 % de variabilidad entre clones aislados de la misma línea celular. [42]. La variabilidad del nivel de expresión puede ser el resultado de fluctuaciones alrededor de un valor particular, o puede indicar que hay muchos puntos estables en ese rango. En cualquier caso, el sistema se mantuvo estable en ese rango. Otro informe se centró en la personalización del modelo de metabolismo celular [43]. En este informe, los parámetros para el modelo cinético de la ruta metabólica de los eritrocitos se determinaron individualmente para 24 personas. La diferencia en los parámetros cinéticos entre individuos fue del 14,5 % como valor medio de las diferencias para todos los parámetros. Las cantidades de moléculas reguladas por dichos parámetros cinéticos varían potencialmente en una escala similar si el sistema es lineal. A pesar de algunas excepciones, como los sistemas no lineales complicados, asumimos que muchos sistemas biológicos pueden tolerar al menos un 10% de cambios en cada variable. En los casos en que esta suposición no sea realista, este valor se puede cambiar considerando información previa sobre el sistema de interés. Por ejemplo, en el caso de un evento biológico, incluidas las transiciones del estado de un sistema a otros estados, como la diferenciación celular de varios pasos. [44], el rango estable de cada estado individual puede ser más estrecho que el de los estados en un sistema uniestable y se debe aplicar la configuración correspondiente en TEAPS. También hay que tener cuidado con la posibilidad de que el sistema tenga una serie de puntos fijos en estos casos. Conceptualmente, cuando identificamos varios estados estacionarios para un determinado sistema e intentamos encontrar conjuntos de parámetros que satisfagan BSR para todos los estados estacionarios, podemos buscar por separado el espacio de parámetros que satisfaga BSR para cada estado estacionario, luego podemos encontrar el espacio superpuesto, pero el desempeño real debe evaluarse como tareas adicionales. Por el contrario, el método TEAPS no será aplicable a sistemas biológicos con comportamientos periódicos autónomos, como el sistema de conducción eléctrica del corazón. [45]sistema generador de ritmo circadiano [46]y regulación de los ciclos celulares [47]. Si se puede diseñar una función objetiva para reflejar los comportamientos periódicos, entonces el algoritmo CNM global modificado utilizado en el presente estudio podría ser aplicable para los análisis de los sistemas sin determinar los valores exactos de los parámetros cinéticos. Además, si hay sistemas que satisfacen BSR pero no son continuos ni diferenciables, TEAPS no se puede aplicar a tales sistemas.

Modificamos el CNM global para permitir el muestreo de múltiples conjuntos de parámetros de todo el espacio de la solución de una manera relativamente imparcial (Figs. 3 y 4). Sin embargo, hay algunos casos con riesgo de búsqueda insuficiente. Uno de esos casos es cuando el verdadero espacio de parámetros permitido es grande, similar al modelo T2, mientras que la búsqueda TEAPS utiliza un número finito de conjuntos de parámetros para la búsqueda espacial. Si el número de conjuntos de parámetros sembrados es inadecuado, el conjunto de parámetros que satisface la restricción BSR puede encontrarse escasamente en todo el espacio de parámetros, lo que dificulta delinear un verdadero espacio de parámetros permisible. Además, el proceso de globalización puede ser insuficiente si el verdadero espacio permitido es grande. Debido a que TEAPS incluye pasos de adición de ruido durante la optimización para minimizar el riesgo mencionado anteriormente, la modificación del proceso de adición de ruido puede mejorar aún más la capacidad de búsqueda en algunos casos. Además, debemos tener cuidado con la posibilidad de que el espacio de parámetros permitido no esté necesariamente conectado. El algoritmo de optimización utilizado en este estudio expande el área de búsqueda adyacente al espacio de solución encontrado en los pasos anteriores del proceso de optimización. Por lo tanto, si se permiten varios espacios de parámetros aislados para que el sistema satisfaga las condiciones de BSR, existe el riesgo de que el algoritmo pierda una parte del espacio de solución. El concepto de encontrar una solución adicional utilizando la solución o los valores ya encontrados en sus regiones vecinas como valor inicial se implementó para superar el problema de los mínimos locales en varios informes. [48,49]. Dado que nuestro enfoque es similar a este concepto, nuestra implementación actual puede tener potencial para encontrar un espacio de solución aislado. Teniendo en cuenta que se incluyen varios procesos estocásticos en el método TEAPS, la repetición de ensayos independientes con diferentes conjuntos de parámetros iniciales podría reducir el riesgo de perder espacios de solución aislados.

Para analizar los comportamientos dinámicos de un sistema, la simulación a menudo se realiza utilizando un conjunto de parámetros determinado ajustando el modelo a las observaciones experimentales [5]; sin embargo, este enfoque tiene limitaciones. Según la estructura del sistema y los datos experimentales disponibles, es posible que conjuntos de parámetros distintos al seleccionado también puedan reproducir las observaciones experimentales. En este caso, la simulación que utiliza el conjunto de parámetros seleccionado no puede reproducir los comportamientos del sistema que no se reflejan en los datos experimentales utilizados para la determinación de parámetros. Este problema se ha categorizado como un problema de sobreajuste que aún no se ha resuelto. [16,50]. Las simulaciones paralelas que utilizan múltiples conjuntos de parámetros podrían ser útiles para evitar este problema, como se propone en este estudio. Aunque es difícil determinar el conjunto de parámetros más cercano a los valores reales de otros conjuntos de parámetros, las características o tendencias comunes observadas en los resultados de la simulación paralela reflejan las propiedades esenciales de los comportamientos del sistema. Este enfoque es similar a la estrategia de embolsado que se usa a menudo en estadística, donde se extrae la tendencia común en las predicciones de una gran cantidad de modelos diferentes (con diferentes parámetros) para evitar el sobreajuste. [51]. Desde una perspectiva diferente, puede ser razonable predecir la dinámica de los sistemas como una superposición de múltiples posibilidades. El valor del parámetro cinético para una molécula en particular cambia según el entorno externo, como el estado de hacinamiento molecular. [9]. En ese sentido, lo que observamos puede ser una superposición de múltiples posibilidades. Además, el análisis de la forma del espacio de parámetros permitido basado en los múltiples conjuntos de parámetros, con los cuales los comportamientos del modelo son plausibles, es útil para comprender las características del comportamiento y la estructura del sistema, como se muestra en los modelos de NF-κB. vía de señalización y vía AA (Figs. 5F y 6D). El análisis de la forma generará parámetros candidatos que influirán significativamente en el comportamiento del sistema. La determinación exacta de dichos parámetros mediante observaciones experimentales puede aumentar la previsibilidad de la dinámica.

Finalmente, proponemos una posible utilidad de nuestro enfoque para predecir el comportamiento de los sistemas. Las simulaciones dinámicas de muchos eventos biológicos a menudo son difíciles debido a la disponibilidad de valores de parámetros cinéticos confiables u observaciones biológicas relacionadas. Incluso en tales casos, nuestro método puede recopilar muchos conjuntos de parámetros bajo la restricción de la BSR. Usando estos conjuntos, puede ser posible encontrar posibles alteraciones en el estado del sistema que pueden ocurrir bajo algunas perturbaciones. Podríamos especular el efecto de las perturbaciones analizando las tendencias de las alteraciones simuladas bajo el supuesto de que el sistema obedece a la restricción de la BSR. Por ejemplo, cuando consideramos la situación en la que la información sobre los efectos de los AINE en la PGE2 no está disponible, a diferencia del caso de Fig. 7, no podemos usar una restricción de concentración de PGE2 para seleccionar conjuntos de parámetros biológicamente plausibles. Incluso en este caso, la simulación de los efectos de los AINE utilizando todos los conjuntos de parámetros adquiridos puede tener algunas implicaciones biológicas. Al simular todos los conjuntos de parámetros, primero encontramos que las simulaciones con muchos conjuntos de parámetros arrojan pequeños cambios en la proporción de PT en la concentración clínica de cada fármaco. Esto es razonable porque una parte de los conjuntos de parámetros aún satisface la condición BSR incluso después de agregar las modificaciones de parámetros correspondientes a las perturbaciones, lo que puede implicar que es poco probable que la exposición de cada fármaco a la concentración clínica cambie la proporción de PT. Luego, cuando nos enfocamos en los conjuntos de parámetros con los que se simuló la proporción de PT para alterar más de cierto nivel (p. ej., 10 %) a una concentración clínica de al menos un fármaco, los resultados simulados fueron consistentes con las observaciones clínicas (Fig. S5) . Bajo la terapia con aspirina en dosis bajas, se encontró que la reducción en el nivel de PGE2 era débil y la proporción de PT era más alta en comparación con las condiciones no tratadas, lo que era consistente con el uso clínico de aspirina en dosis bajas como terapia antiplaquetaria. [52]. En el caso de rofecoxib, el valor mediano de los cocientes de PT fue el más bajo, lo que es consistente con el alto riesgo de eventos cardíacos con rofecoxib [53]. Estos resultados de simulación suposicionales respaldan que nuestro enfoque puede ser útil para extraer los comportamientos dinámicos de los sistemas biológicos cuyos parámetros cinéticos no están disponibles.

Colectivamente, la suposición de las condiciones BSR y la adquisición de tantos conjuntos de parámetros plausibles como sea posible es una nueva estrategia efectiva para comprender los comportamientos dinámicos de un sistema biológico.

materiales y métodos

Preparación del archivo modelo

Todos los modelos se prepararon como archivos XML compatibles con SBML utilizando Cell Designer 4.4 [54].

Modelo T1.

Las siguientes son las ecuaciones que indican las velocidades de las reacciones:

El sistema ODE se compone de la siguiente manera:

Modelo T2.

Las siguientes son las ecuaciones que indican las velocidades de las reacciones:

El sistema ODE se compone de la siguiente manera:

Modelo T3.

Las siguientes son las ecuaciones que indican las velocidades de las reacciones:

El sistema ODE se compone de la siguiente manera:

Modelo T4.

Las siguientes son las ecuaciones que indican las velocidades de las reacciones:

El sistema ODE se compone de la siguiente manera:

Modelo T5.

Las siguientes son las ecuaciones que indican las velocidades de las reacciones:

El sistema ODE se compone de la siguiente manera:

Modelo T6.

Las siguientes son ecuaciones que indican las velocidades de las reacciones:

El sistema ODE se compone de la siguiente manera:

Modelo T7.

Las siguientes son las ecuaciones que indican las velocidades de las reacciones:

El sistema ODE se compone de la siguiente manera:

Modelo T8.

Las siguientes son las ecuaciones que indican las velocidades de las reacciones:

El sistema ODE se compone de la siguiente manera:

Modelo de vía NF-κB

La estructura del modelo se refiere a la de un informe anterior [33,34]. También se proporciona el archivo de modelo SBML (archivo S1).

Algoritmo de TEAPS

Se muestra un esquema en Figura S1 y algoritmo detallado en Información S3. Diseñamos el algoritmo para una búsqueda exhaustiva de conjuntos de parámetros que satisfacen las restricciones de BSR cuando se les da un rango de búsqueda para cada parámetro. Optimizamos las funciones objetivo de BSR utilizando el método de Newton de clúster global (gCNM) [17] con varias modificaciones. El CNM global consta de dos etapas, la etapa CNM y la etapa del método cuasi-Newton globalizado. En la etapa CNM, en lugar de realizar múltiples iteraciones hasta la convergencia estricta, modificamos para recopilar los conjuntos de parámetros que cumplieron con las condiciones moderadas en cada iteración e iteramos hasta que el número de conjuntos de parámetros recopilados alcanzó el número preestablecido. En la etapa del método cuasi-Newton globalizado, se implementa originalmente el método de Broyden con modificación para la búsqueda global. En este estudio, adoptamos la optimización basada en BFGS de memoria limitada (LBFGS) para reducir el consumo de memoria durante la optimización considerando la situación de aplicación a modelos con muchos parámetros. También modificamos el método de globalización. En el método original, la dirección de cambio en cada paso de iteración en el método cuasi-Newton se torció para permitir la búsqueda global. Describimos LBFGS con esta modificación como m-LBFGS. Además de m-LBFGS, agregamos la adición de ruido de paso a los conjuntos de parámetros obtenidos y el bucle de paso de reoptimización (g-LBFGS). Repetimos g-LBFGS expandiendo el tamaño de una región para muestrear puntos de observación para evaluar Ocuenca. Este gCNM modificado, que consta de la etapa CNM y la etapa g-LBFGS, se repitió usando diferentes conjuntos de parámetros iniciales hasta que la distribución de los conjuntos de parámetros no cambió significativamente por la inclusión de los conjuntos de parámetros encontrados por la última iteración de gCNM. Para juzgar la importancia, implementamos dos condiciones. Una es la prueba de suma de rangos de Wilcoxon entre una distribución de todos los valores para cada parámetro hasta la iteración anterior y otra obtenida al incluir la última iteración. En la implementación actual, el nivel de significación (valor alfa) fue de 0,1. La otra es que para todos los parámetros, la mayoría de los valores de parámetros encontrados en la última iteración de gCNM (específicamente, más del 99 % en la configuración actual) están dentro de la distribución del conjunto de parámetros ya encontrado. Para obtener conjuntos de parámetros de manera efectiva, se requiere ajustar los pesos para la función objetivo. Encontramos el tratamiento equivalente de (Oarreglar)2(Ocuenca)wy (Orelax)w donde el peso w = 1~2 es aplicable en muchos casos como regla general.

Ajuste de TEAPS utilizado en este estudio

El espacio de búsqueda de parámetros se determinó arbitrariamente. Para los modelos T1 a T8, el rango de 5 × 10−3 a 5 × 104 se utilizó a excepción de un flujo de entrada. El rango de flujo de entrada se estableció en 5 × 10−1 a 5 × 104. Para el modelo de la ruta del ácido araquidónico, el número de reacciones secuenciales en el modelo, que es una especie de estructura de red profunda, es mayor que en los modelos T1 a T8, por lo tanto, se permitió la mayor diferencia relativa de flujos en la ruta del ácido araquidónico y el rango de 1 × 10−3 a 1 × 1024 se utilizó. En las funciones objetivo se utilizaron la siguiente función e hiperparámetros. Para la restricción de punto fijo, Oarreglar(m) = ‖F(X*, m)‖2usamos punto fijo objetivo X* = (1,⋯,1)T como se menciona en el texto principal. en el calculo de se generaron 20 puntos por X(k). La distribución de X(k) se amplió a medida que TEAPS progresa como se describe en 7) en el algoritmo de TEAPS. en el calculo de se usó en lugar de Re(λ*) ≠ 0 como errores de cálculo tolerables. Finalmente, el parámetro que satisface ‖F(X*, m)‖2−6 y fue juzgado como conjuntos de parámetros que convergen a X* y se utiliza en análisis posteriores.

Estimación de la estabilidad de la cuenca

Muestreamos aleatoriamente 500 conjuntos de parámetros en cada paso de TEAPS. Para cada conjunto de parámetros, la integración ODE se realizó desde 500 estados iniciales alrededor del punto fijo objetivo hasta el tiempo t = 100. Para cada estado inicial, usando el estado en t = 100 como valor inicial, un estado que minimiza la distancia entre y (0,…,0)T se buscó utilizando la función «fsolve» en Matlab y se trató como un punto fijo. En el presente estudio, se generaron estados iniciales para Xi incluido Xyo = 1+adónde a se toma de una distribución uniforme sobre [−0.15, 0.15] para modelos pequeños. Para el modelo NF-κB, considerando un aumento de la dimensión de la variable de estado, se generaron estados iniciales en una región menos pequeña, tal que a se toma de una distribución uniforme sobre [−0.105, 0.105]. Para el modelo AA, los estados iniciales se generaron en una región menos pequeña, tal que a se toma de una distribución uniforme sobre [−0.101, 0.101]. Si todos los componentes de un punto fijo obtenido están en rangos entre el 90% y el 110% de los valores objetivo, la combinación del estado inicial y el conjunto de parámetros se consideró estable. Para cada conjunto de parámetros, se calculó una relación de estados iniciales estables. Las medias y los errores estándar de las proporciones de todos los conjuntos de parámetros muestreados se calcularon en cada paso de TEAPS.

Búsqueda de fuerza bruta

Los puntos de cuadrícula para la búsqueda de fuerza bruta se generaron a partir de una distribución uniforme que va desde 5 × 10−4 a 5 × 104 en una escala logarítmica. Para reducir el número de puntos, las restricciones para dar el punto fijo Xss = 1 y se redujo el número de parámetros libres. El número de puntos de cuadrícula generados fue de 3 × 108 a 3 × 109 puntos. Para cada punto de cuadrícula generado, los valores de Ocuenca(m) y Orelax(m) se calcularon para juzgar si el punto satisface la condición BSR. Para el cálculo del modelo T8, para reducir el número de puntos de cálculo y evitar la dispersión global, solo se fijó el flujo de entrada en 0,5, que fue el valor observado con mayor frecuencia en los resultados de TEAPS. Además de la búsqueda de puntos de cuadrícula, también comparamos la capacidad de búsqueda utilizando vectores aleatorios muestreados del espacio de parámetros bajo la restricción de que los vectores deben satisfacer Oarreglar(m) = 0.

Cálculo de la matriz jacobiana y sus valores propios

Se extrajeron las descripciones de ODE en modelos SBML. Usando las ecuaciones extraídas, se generaron matrices jacobianas usando la función “jacobiana” en Maxima. Luego, la salida se transformó al formato de función de MATLAB para poder calcular la matriz jacobiana para cada modelo de manera analítica. Los valores propios de la matriz jacobiana se calcularon utilizando la función «eig» en MATLAB.

Evaluación de la similitud de distribución

La estimación de la densidad del kernel se realizó utilizando la función «mvksdensity» en MATLAB. El ancho de banda de la función kernel se calculó mediante la regla de Silverman [55]. Se vectorizaron los tensores que componen la densidad estimada en los puntos de la cuadrícula en un espacio de parámetros multidimensionales. La similitud entre los vectores por TEAPS y la búsqueda de fuerza bruta se evaluó en los siguientes dos indicadores. Para evaluar la similitud de las dos distribuciones de conjuntos de parámetros, se calcularon la similitud del coseno y la divergencia JS. La similitud del coseno se calculó de la siguiente manera: dónde Ai y Bi son las componentes de los vectores a comparar. Este valor oscila entre 0 y 1 y toma el valor de 1 cuando se comparan vectores completamente iguales. La divergencia JS se calculó de la siguiente manera: dónde Ai y Bi son las componentes de los vectores a comparar. Este valor puede ser 0 o más y toma un valor de 0 cuando se comparan vectores completamente iguales. Para evaluar la superioridad del muestreo TEAPS frente al muestreo aleatorio, los vectores de comparación se prepararon de la siguiente manera: (1) se prepararon histogramas de densidad del vector de densidad de fuerza bruta para cada modelo, (2) cada componente de los vectores de comparación se muestreó aleatoriamente utilizando la densidad histograma de cada modelo como una distribución de probabilidad. Al repetir el muestreo aleatorio, se calcula el intervalo de confianza del 95% de la relación de cada índice de similitud en una comparación entre fuerza bruta y TEAPS con respecto a una comparación entre fuerza bruta y muestreo aleatorio para cada modelo.

Análisis de agrupamiento

Se realizó un agrupamiento jerárquico de las proporciones de PT simuladas. Se preparó un vector que constaba de una serie de resultados simulados utilizando múltiples conjuntos de parámetros para cada fármaco. Se calcularon los coeficientes de correlación entre los vectores de todos los pares de fármacos. El análisis de conglomerados basado en los coeficientes de correlación se realizó mediante el análisis de conglomerados de varianza mínima de Ward utilizando funciones de «vinculación» en MATLAB.

Cálculo de índices del eje central principal (PCIs)

Los PCI se calcularon de la siguiente manera: dónde yok,metro es la carga de la metroth parámetro en el kel componente principal, y nortePCA es la dimensión del espacio de parámetros.

Integración ODE

Las ODE se resolvieron utilizando el solucionador SUNDIALS CVode implementado en MATLAB o los solucionadores rígidos de ODE de MATLAB.

Plataforma de cómputo

Los cálculos simbólicos se realizaron con Maxima. Todos los cálculos numéricos se realizaron con MATLAB utilizando un Mac Pro (CPU: Xeon E5-2697 v2 2,7 GHz × 1, SO: mac OS Catlina, RAM: 64 GB, versión de MATLAB: R2019a), o una fuente de computación en la nube (AMAZON EC2, tipo de instancia: c5.12xlarge, número de CPU virtuales: 48, RAM: 96 GB, sistema operativo: AlmaLinux 8.5 de 64 bits, versión de MATLAB: R2019a, R20221b o R2022a).

Información de soporte

S1 Fig. Diagrama esquemático del algoritmo TEAPS.

TEAPS es un algoritmo para buscar el espacio de parámetros globales que satisface BSR iterando método de Newton de clúster global modificado (CNM). Nuestro CNM global modificado compuesto por el paso CNM y el paso global-LBFGS (g-LBFGS). El método g-LBFGS compuesto por el método L-BFGS con modificación de búsqueda globalizada (m-LBFGS) y pasos de adición de ruido. Para optimizar eficientemente las funciones objetivo de BSR en el paso g-LBFGS, la función objetivo se cambia en cada subproceso como se muestra.

https://doi.org/10.1371/journal.pcbi.1010441.s001

(TIF)

S2 Fig. Método CNM global modificado utilizado en nuestro estudio.

(A) (ac) El CNM se utilizó para muestrear los conjuntos de parámetros cerca de las situaciones objetivo. Los conjuntos de parámetros iniciales generados aleatoriamente (círculos rojos) se utilizaron como conjuntos de parámetros iniciales para CNM (a). Durante la iteración de CNM, algunos conjuntos de parámetros se acercan a los óptimos, y dichos conjuntos de parámetros (círculos verdes con contorno rojo) se recopilaron (b, c). (B) (dg) Los conjuntos de parámetros muestreados se utilizaron como conjuntos de parámetros iniciales para las etapas de expansión y optimización. Los conjuntos de parámetros muestreados se aplicaron al método cuasi-Newton modificado y se obtuvieron nuevos conjuntos de parámetros (círculos naranjas) (d). Para expandir la distribución, se agregaron ruidos a los nuevos conjuntos de parámetros (círculos azules con contorno negro) (e). Estos conjuntos de parámetros se aplicaron aún más al método cuasi-Newton modificado y se obtuvieron los conjuntos de parámetros renovados (círculos naranjas con contorno negro) (f). Para expandir la distribución de los conjuntos de parámetros, las subrutinas que consisten en la adición y optimización de ruido se repitieron varias veces (g). Finalmente, se utilizó el método normal de cuasi-Newton para optimizar el punto de solución más cercano de cada conjunto de parámetros. Por lo tanto, la distribución de conjuntos de parámetros se amplió a nivel mundial.

https://doi.org/10.1371/journal.pcbi.1010441.s002

(TIF)

S3 Fig. Visualización de la estructura de conjuntos de parámetros que cumplieron la condición BSR para el modelo T8.

Aunque las leyes cinéticas para el modelo T8 son complicadas debido a las múltiples regulaciones sobre una reacción, TEAPS ha logrado buscar la forma completa. Cada conjunto de parámetros se indica con un círculo. Los círculos azules indican los conjuntos de parámetros encontrados por TEAPS y los círculos rojos indican los conjuntos de parámetros encontrados por la búsqueda de fuerza bruta.

https://doi.org/10.1371/journal.pcbi.1010441.s003

(TIF)

S4 Fig. Comparación de la distribución entre el espacio de parámetros BSR encontrado por búsqueda basada en puntos de cuadrícula o basada en vectores aleatorios y uno por TEAPS.

Los histogramas de Grid point y TEAPS son los mismos datos que el texto principal. La búsqueda basada en vectores aleatorios se agrega en esta figura. Cuando realizamos una búsqueda usando vectores aleatorios, cada elemento de los cuales fue generado por distribución uniforme en el espacio de búsqueda, y se obtuvieron resultados similares para todos los modelos. En la generación de vectores aleatorios, se fijaron algunos parámetros ya que la restricción de punto fijo para cada modelo fue dada analíticamente de antemano. El número de vectores generados fue 3 × 108 a 3 × 109 puntos. Para todos los modelos de muestra, confirmamos que los tres histogramas mostraban distribuciones similares.

https://doi.org/10.1371/journal.pcbi.1010441.s004

(TIF)

S5 Fig. Los conjuntos de parámetros obtenidos por TEAPS ayudan a predecir la en vivo fenotipo asociado a la administración de AINE.

(A, B) Como un experimento in silico hipotético, asumimos la situación en la que los efectos de la administración de AINE en las concentraciones de PGE2 eran completamente desconocidos. Esta suposición se estableció para imitar modelos biológicos con información limitada o nula sobre el comportamiento cuantitativo de los sistemas. Se simularon las situaciones en las que cada AINE se administró a la concentración promedio de uso clínico normal. Luego, para detectar tantas alteraciones significativas por los AINE como fuera posible, nos enfocamos en conjuntos de parámetros con los que la relación PT simulaba alterar más del 10 % de al menos un fármaco. Cada punto azul corresponde a un resultado de un conjunto de parámetros obtenido por TEAPS. Se muestran el nivel de PGE2 en PMN (A) y la proporción de PT (B), un marcador del gasto fisiológico. Incluso en simulaciones con información limitada, la tendencia de la relación PT es comparable a las observaciones clínicas, lo que respalda que TEAPS ayuda a predecir el comportamiento de los sistemas biológicos. Los recuadros verdes indican los valores medianos.

https://doi.org/10.1371/journal.pcbi.1010441.s005

(TIF)

Fuente del artículo

¿Que te ha parecido?

Deja un comentario