Introducción

La diabetes mellitus (DM) es una enfermedad crónica derivada de la incapacidad del páncreas al no secretar suficiente insulina o la incapacidad del cuerpo para utilizarla eficazmente. La insulina es la hormona que regula la concentración de glucosa en la sangre (glucemia). La consecuencia común asociada a la DM no controlada es la hiperglucemia, es decir, altos niveles de glucemia. En el caso particular de la diabetes tipo 1 (DMT1) es el resultado de un ataque autoinmune a las células beta pancreáticas (células β) que secretan insulina^{[1] [2]}.

Actualmente, el modelado matemático aplicado en el análisis de patologías humanas es una herramienta utilizada, en el área de la medicina, para mejorar la comprensión de un sistema fisiológico con el objetivo de prevenir futuras situaciones de enfermedades^[3]. Sus revisiones aportan una descripción más detallada de esta enfermedad y son una herramienta que contribuye a su diagnóstico y prueba de tratamientos, gracias a las simulaciones basadas en estos^[4]. A causa de esto, en las últimas décadas se han realizado esfuerzos para el diseño y análisis de modelos matemáticos que describan su dinámica.

Es conocido que la DM ha alcanzado cifras preocupantes a nivel mundial y al no tener cura, ha motivado al desarrollo de investigaciones sobre métodos para su diagnóstico, control y tratamiento. Sin embargo, no todos estos modelos cumplen con el principio de positividad, por lo que es necesario determinar si un modelo cumple con él para asegurar que el comportamiento del sistema siempre se mantendrá en los valores positivos para todo tiempo futuro, que generalmente ocurre en los sistemas fisiológicos (por ejemplo, no existen valores negativos de glucosa en sangre).

La actividad emocional, consecuencia a eventos de estrés, es considerada como uno de los factores asociados en la desestabilización del control metabólico, actuando sobre los niveles de glucemia a través de la movilización de las “hormonas del estrés” (como la epinefrina) e indirectamente por su impacto sobre las conductas de adherencia al tratamiento^[5], la cual se considera común en pacientes con esta enfermedad y provoca que sus niveles de glucemia aumenten, haciéndola compleja de controlar. En^[6], se reporta que la pandemia de enfermedad por coronavirus 2019 (COVID-19), ha impactado en DM con una mayor tasa de incidencia, así como formas de exacerbación de la hiperglucemia y sus complicaciones en pacientes. Aunado a ello, se reporta hiperglucemia inducida por estrés en muchos pacientes no diabéticos hospitalizados después de contraer COVID-19. Implicando que nuestro análisis matemático contribuye en continuar estudiando el comportamiento de la epinefrina en un modelo no lineal mediante un conjunto de ecuaciones diferenciales ordinarias.

Es importante destacar que recientemente en^[7], se ha reportado un modelo que conduce a la predicción del nivel de glucemia, donde tales predicciones tienen el potencial de ser utilizadas como parte de sistemas de control las cuales pueden ser robustos para modelar imperfecciones y datos ruidosos. Lo interesante de este modelo es que su capacidad predictiva se compara con modelos enfocados a diabetes mellitus tipo 2 (DMT2). Así mismo, en^[8] se reporta un análisis global enfocado a las biomatemáticas para continuar estudiando las hipótesis con las que algunos modelos de DM se definen con base a datos experimentales o in-silico. Es por ello, la inquietud de continuar estudiando los modelos con el que se predice DMT1 se convierte en un desafío vigente y, por ende, el análisis matemático mediante el método de Localización de Conjuntos Compactos Invariantes fortalece las hipótesis en el desarrollo de cada modelo que se ha propuesto en la literatura. Un avance reportado en^[9], se discute un esquema de control activo para forzar que el régimen caótico de DMT1 siga un estado oscilatorio periódico, es decir, un equilibrio libre de desorden.

Por último, cabe destacar que el estudio de su comportamiento no lineal no ha sido reportado en la literatura y es donde impactará este análisis, el cual permite observar y estudiar la evolución del modelo a través del tiempo. Determinar la importancia de la influencia de la epinefrina dentro esta dinámica asociada a la DMT1; así como también determinar condiciones que sirvan como base para proponer y diseñar nuevos modelos matemáticos relacionados con la DMT1, considerando estos parámetros.

El estudio de la dinámica a largo plazo es de gran aportación dentro del campo de investigación ya que ayudan a encontrar o probar nuevas formas de diagnóstico y tratamiento de esta enfermedad. Comprueba si los modelos cumplen con el principio de positividad y determina las condiciones que permiten establecer un preámbulo matemático donde se evidencia la importancia de la epinefrina en esta dinámica asociada a la DMT1.

Materiales y métodos

En esta sección se presentan los conocimientos necesarios para la aplicación del método de Localización de Conjuntos Compactos Invariantes (LCCI) y su aplicación para determinar el dominio de localización que contiene todos los conjuntos invariantes del sistema en estudio. Este método permite conocer el comportamiento en el corto y largo plazo de sistemas dinámicos descritos por Ecuaciones Diferenciales Ordinarias (EDOs) de primer orden^[13]. En años recientes, ha sido aplicado en el estudio de diferentes sistemas dinámicos tales como sistemas biológicos que describen el desarrollo y evolución de enfermedades como el cáncer o la diabetes^{[14] [15] [16]}. El método es estrictamente analítico, lo que implica que es posible encontrar la solución del problema sin la necesidad de resolver el sistema de EDOs por algún método numérico. Sin embargo, la solución numérica es utilizada para ilustrar los resultados matemáticos obtenidos, es decir, el dominio de localización y los conjuntos compactos invariantes del sistema.

Resultados y discusión

En la Tabla 1, se muestra el tipo de estabilidad correspondiente a cada punto de equilibrio del modelo (1)-(3) después de haber aplicado el método indirecto de Lyapunov. La columna relacionada a Signos de Reλ_i, implica el signo correspondiente al valor propio obtenido por la ecuación característica.

La complejidad con el que se describe el modelo^[2] es discutido solo con el análisis de los puntos de equilibrio por lo que, la aplicación del método LCCI permite conocer las concentraciones máximas de G, I y β derivado de la influencia de la epinefrina para proponer algún tipo de controlador.

Aunado a ello, se realiza la prueba de positividad al modelo (1)-(3), es importante mencionar que en las ecuaciones (8) y (9) el valor de G*=dGdtG=0 ; satisfaciendo el Lema 1.

Ahora, se aplica el método LCCI proponiendo a la función localizadora ℎ₁ = 𝐺, donde después de aplicar el Teorema 1 se obtiene el siguiente conjunto de localización:

definiendo un G_mín = 0 y un G_máxR0+GeEG0. Posteriormente se propone una segunda función localizadora h₂ = q₁G + q₂β, el cual mediante la aplicación del Teorema 2 se obtiene el conjunto de localización definido como:

donde h₂|_S(h2) ⋂K₁ = q₁G_máx + q1R0+Ged0, implicando que β_mín = 0 y β_máx≤q1R0+GeEG0+q1R0+Ged0.

Es importante mencionar que dentro de la función se encuentra los parámetros libres q₁ y q₂, siendo ambos parámetros positivos y puede tomar cualquier valor. Finalmente, se propone la última función localizadora h₃ = I, debido a que se cuenta con G_máx y β_máx, esto facilita definir, mediante el Teorema 2, una región acotada para I, la cual queda definida por el siguiente conjunto:

identificando como I_mín = 0 y un I_máx≤βmaxσGmáx2α+ρ+k.

Por lo tanto, derivado de la aplicación del método de LCCI, se establece el siguiente teorema:

Teorema 3. La localización de todos los conjuntos compactos invariantes se obtiene a través de funciones localizadoras dependientes de las variables de estado del sistema, denotado por (1)-(3), definiendo el dominio de interés del conjunto K = {K₁ ∩ K₂ ∩ K₃}, donde K = G_máx ∩ I_máx ∩ β_máx debido a las funciones localizadoras h_i para i = 1,2,3. Donde la región de localización se satisface bajo la restricción de R³_+,0. Dado que R³_+,0 = G(t) ≥ 0, I(t) ≥ 0, β(t) ≥ 0.

Simulaciones numéricas

Considerando al punto de equilibrio P₁ de la Tabla 1, se observa en la Figura 1 que todas las trayectorias convergen a dicho punto cuando se considera la condición inicial para G(0) =2 00, I(0) = 4.5, β (0)= 65; además, se puede ver como sus trayectorias convergen al punto de equilibrio P₁, validando que se encuentra en el ortante positivo y es asintóticamente estable en forma local, de acuerdo con el método indirecto de Lyapunov. En la Tabla 2, se muestra el valor numérico relacionado a cada parámetro en el modelo (1)-(3). Para el desarrollo de la etapa de simulación se aplica el comando de ODE45 de MatLab para la comprobación de resultados, el cual representa el uso del algoritmo Runge-Kutta de cuarto y quinto orden, utilizando un paso de integración de 0.001.

Figura 1 Comportamiento de las trayectorias hacia el punto de equilibrio P₁, cuando se considera la condición inicial para G(0) = 200, I(0) = 4.5, β(0) = 65. 

Ahora, para pasar a la validación de los resultados del método de LCCI. En la Figura 2 se muestra el dominio de localización K = {K₁ ∩ K₂ ∩ K₃} cuando se cumple la condición d₀ + r₂G₂ > r₁G, de acuerdo con los valores establecidos en la Tabla 2. Además, se presenta la primera vertiente del modelo, donde se considera G = G_mín en la concentración máxima de β, es decir βmáx≤q1R0+GeEG0+q1R0+Ged0.

Figura 2 Dominio de localización para cada variable de estado cuya trayectoria se encuentra contenida, es decir, se muestra el dominio de localización del modelo (1)-(3) implicando que se tiene un G_máx, I_máx y β_máx; satisfaciendo el Teorema 3 mediante el Teorema 1 y Teorema 2. 

Utilizando las mismas condiciones iniciales que en la Figura 1, se muestran las trayectorias de G(t), I(t) y β(t), que se encuentran por debajo de los valores obtenidos para cada concentración máxima que en este caso son G_máx ≤ 697.22, I_máx ≤ 42326.15 y β_máx ≤ 17430.55. Así mismo, en el la Figura 3 se muestra que todas las trayectorias se encuentra en el ortante positivo y es asintóticamente estable de forma local; además, se observa el dominio de localización que contienen a las trayectorias del modelo (1)-(3), satisfaciendo el Teorema 3 cuando se considera el punto de equilibrio P₃.

Figura 3 Comportamiento de las trayectorias hacia el punto de equilibrio P₃, cuando se considera la condición inicial para G(0) = 500, I(0) = 1.5 , β(0) = 18. 

Conclusiones

El modelo analizado en este estudio contiene tres puntos de equilibrio en el ortante positivo, los cuales son influenciados por el aumento de glucosa debido a la secreción de epinefrina (Ge) y la eficacia de la epinefrina en suprimir la secreción de insulina (𝜌); a través de las funciones localizadoras se obtienen las condiciones de su dominio de localización definido por K = {K₁ ∩ K₂ ∩ K₃}. Este tiene dos vertientes para la condición de β_máx en presencia de epinefrina. Donde las condiciones de este modelo se ven influenciadas por los parámetros Ge y 𝜌, por lo que se muestra un aumento de sus valores con respecto a los establecidos en el modelo antecesor^[6].

Es importante mencionar que la ecuación que corresponde a la masa de células β (dβ) tiene cambios de signo y es igual para ambos modelos, porque representa la muerte de células β y su replicación (considerando la reducción cuando hay hiperglucemia) mediante un polinomio que no incluye la influencia de la epinefrina. Esta característica particular provoca que haya limitaciones al proponer y seleccionar la función localizadora h₂, esto debido a que la mayoría llevan a un resultado donde la condición de β_máx es negativa, que biológicamente no es posible.

De acuerdo a lo establecido en el marco teórico del modelo^[2], un suministro de epinefrina (causado por activación emocional) produce un incremento en la concentración de glucosa y disminuye la secreción de insulina por las células β. Sin embargo, al obtener los puntos de equilibrio y las condiciones de su dominio de localización de este modelo (sus concentraciones máximas de G, I y β) se demuestra que hay un aumento en los niveles de tres variables en comparación con los obtenidos del modelo de Topp donde, se esperaría que los valores de I fueran menores.

Por último, para llevar a cabo las simulaciones numéricas se consideran condiciones iniciales que son cercanas a sus puntos de equilibrio. Se observa que el comportamiento de las variables con respecto al tiempo siempre converge al mismo punto de equilibrio, el cual se considera asintóticamente estable en forma local (P₁ o P₃) con esas condiciones y las trayectorias son similares entre modelos por lo que, la epinefrina no provoca grandes cambios en esto. Sin embargo, la influencia de la epinefrina en esta dinámica provoca un aumento en el valor de los puntos de equilibrio y las concentraciones máximas de G, I y β, como puede ser observado en las Figuras 1-2.

Como trabajo a futuro se pretende analizar otra variante del modelo, considerando la propuesta en^[19]. Además, plantear otros valores positivos para los parámetros q₁ y q₂ de la función localizadora h₂, los cuales se consideran como q₁ = q₂ = 1 y determinar su influencia en los niveles máximos de cada concentración. Adicionalmente, se buscará analizar los modelos considerando la presencia de una entrada de control en la variable de insulina (I), para determinar si existe la posibilidad de diseñar una ley de control considerando el parámetro de epinefrina. Este trabajo de investigación se encuentra distante de implementar hipotéticamente un páncreas artificial mediante una discusión matemática, sin embargo, las implicaciones matemáticas obtenidas contribuyen para el desarrollo de algún entrenamiento computacional^[20], el cual permita diseñar una ley de control considerando el parámetro de epinefrina. Así mismo, se establece un preámbulo matemático el cual permite ser una referencia para el simulador UVA/Padova^[21]; analizando el comportamiento de la epinefrina considerando a un conjunto de células pancreáticas.

Contribución de los autores

D. G. conceptualización, análisis formal del modelo, estado del arte, escritura y edición del borrador del manuscrito. P. J. C. validación, organización de información, edición, revisión y correcciones del manuscrito.