Numerical solution of biological reaction diffusion models on fixed domains by the finite element method

Libardo Andrés González, Juan Carlos Vanegas^* , Diego Alexander Garzón

Grupo de Modelamiento Matemático y Métodos Numéricos GNUM. Facultad de Ingeniería. Universidad Nacional de Colombia, Sede Bogotá, Carrera 30 Calle 45, Ciudad Universitaria, Bogotá, Colombia

Resumen

Múltiples fenómenos biológicos se han descrito mediante modelos matemáticos formulados a partir de ecuaciones de reacción difusión. La solución de este tipo de ecuaciones da lugar a la formación de patrones espacio-temporales que se ajustan a la realidad biológica del fenómeno modelado. En este artículo se describe la implementación numérica de tres modelos de reacción difusión bien referenciados: el modelo de morfogénesis de Schnakenberg, y los modelos de reacción cinética de Gierer-Meinhardt y Thomas. El objetivo es analizar el conjunto de parámetros asociados con la formación de los patrones espaciotemporales. La implementación numérica se realiza utilizando el método de los elementos finitos en dominios unidimensionales y bidimensionales. Se concluye que la formación de patrones espacio-temporales en modelos de reacción difusión depende de los parámetros constantes del modelo, de las condiciones iniciales y de la técnica de implementación. El análisis de estas dependencias es útil para la formulación y validación de nuevos modelos matemáticos que describan fenómenos biológicos.

Palabras clave: Modelos de reacción difusión, formación de patrones, método de elementos finitos, biología matemática

Abstract

Keywords: Reaction diffusion models, pattern formation, finite element method, biomathematics

Introducción

La formación de patrones biológicos espacio temporales es descrita a través de diversos tipos de modelos matemáticos [1-3]. Uno de estos tipos de modelos corresponde a los modelos basados en ecuaciones de reacción difusión. Estos modelos describen las interacciones químicas que generan patrones complejos en el espacio y/o el tiempo como la interacción de términos de transporte, síntesis y degradación que dependen de las sustancias químicas presentes en el dominio de análisis [2, 4]. En 1952 A. Turing [5] fue el primero en observar y atribuir a las interacciones químicas la auto-formación de patrones en la naturaleza y estudió las soluciones de los modelos biológicos descritos por ecuaciones de reacción difusión. En ellas encontró que pueden existir tres tipos de inestabilidades: a) oscilatorias en el tiempo y uniformes en el espacio, relacionadas con las inestabilidades de Hopf independientes del espacio, b) estacionarias en el tiempo y periódicas en el espacio, y c) oscilatorias en el espacio y en el tiempo [6]. Turing demostró además que un sistema químico de reacción difusión puede evolucionar hacia patrones espaciales heterogéneos desde un estado estacionario uniforme en respuesta a pequeñas perturbaciones [1, 2, 3, 7]. Conforme a esto, estableció que la difusión puede llevar un sistema químico a la inestabilidad, induciendo así la formación de un patrón espacial donde antes no existía. Este tipo de inestabilidad, estacionaria en el tiempo, es más conocida como inestabilidad de Turing [3, 7, 8]. En el caso biológico, la formación de patrones está relacionada con la distribución de sustancias que reaccionan y se difunden en cierta geometría [4, 8, 9]. En el desarrollo tisular, por ejemplo, estas sustancias son consideradas como marcadores que contienen alguna información que el tejido necesita para su crecimiento, formación y maduración [10-12]. Tal es el caso de un organismo en etapa embrionaria, donde sus extremidades están en crecimiento y determinadas sustancias pueden indicarle en qué zonas y en cuáles direcciones debe crecer [13]. La complejidad de estos modelos deriva en la dificultad de obtener soluciones analíticas, debido condiciones de no linealidad y la dependencia respecto a la geometría del dominio de solución [14]. Debido a esto, en diversas aplicaciones se obtiene una solución mediante el uso de métodos numéricos implementados en computador [2, 3, 9]. En este artículo se describe la implementación numérica del conjunto de ecuaciones de reacción difusión que componen tres modelos biológicos conocidos: el modelo de morfogénesis de Schnakenberg, y los modelos de reacción cinética de Gierer-Meinhardt y Thomas [2, 8, 9, 15]. Las soluciones se obtuvieron utilizando la discretización en espacio a través del método de los elementos finitos, y la discretización en tiempo por el método Backward Euler, analizando tanto el caso unidimensional como el bidimensional. En la siguiente sección se describe el modelamiento por ecuaciones de reacción difusión. Luego se hace mención a los modelos biológicos analizados y a la técnica de implementación utilizada. Finalmente, se presentan y discuten los resultados obtenidos, y se formulan algunas conclusiones.

Ecuaciones de reacción difusión

De la misma manera que una gota de tinta evoluciona en el tiempo en un recipiente con agua, ciertas sustancias se mueven en sus medios mediante un proceso que se denomina difusión. Algo similar ocurre cuando el extremo de una barra, inicialmente con una distribución de temperatura homogénea, se lleva a una temperatura superior a la de cualquier otra zona de la barra. Se puede decir que la temperatura se difunde a lo largo de la barra. Por otro lado, cuando dos o más sustancias se encuentran en un medio cualquiera, se puede producir una reacción entre las mismas caracterizada por procesos de consumo, destrucción o transformación. Si se consideran los dos fenómenos conjuntamente, es decir, la difusión y la reacción, cada sustancia se puede modelar a partir de ecuaciones de reacción difusión, y al conjunto de dos o más sustancias se le denomina sistema de reacción difusión [7, 14]. Sea u(t,x) la concentración de una sustancia, con t el tiempo y x la variable espacial. El movimiento de u(t,x), o término difusivo, indica los cambios en la concentración desde puntos de mayor concentración hacia puntos de menor concentración [15-17]. Este principio es conocido como Ley de Fick y se expresa de la siguiente forma (1), donde J es el vector de flujo de u(t,x), y D es el coeficiente de difusión.

A su vez, la reacción entre dos o más sustancias establecen un término reactivo adicional en la función de concentración u(t,x) denotado por f(t,x,u). De acuerdo al principio de conservación, la razón de cambio de la cantidad de materia contenida en un volumen V debe ser igual al flujo neto de materia a través de la superficie S que la delimita, más la cantidad de materia transformada al interior de V debido al término reactivo. Esto expresado matemáticamente es (2):

]]> En (2) n es el vector normal a la superficie S. Utilizando el teorema de la divergencia en el término difusivo y combinando (1) y (2) se obtiene (3):

La ecuación (3) corresponde a una ecuación diferencial definida en el dominio Ω=V con unas condiciones de contorno definidas por la superficie ‎ =S que rodea al volumen V. Para garantizar que el patrón espacial formado se deba únicamente a la organización al interior del contorno, y no a flujos externos, se deben asumir condiciones de flujo en el contorno iguales a cero [14]. Expresando (3) en forma diferencial se obtiene (4):

La ecuación (4) se conoce como ecuación de reacción difusión y permite predecir la evolución de la sustancia denotada como u(t,x). Si se tienen dos sustancias, el sistema de ecuaciones de reacción difusión puede escribirse en su forma adimensional como (5):

En la ecuación (5) u(t,x) y v(t,x) son las concentraciones de las dos sustancias, f(u,v) y g(u,v) son los términos reactivos y d es el coeficiente de difusividad dado por la adimensionalidad del sistema.

Inestabilidad de Turing

Supóngase que en (5) los términos difusivos son cero y que existe un estado estacionario estable. Si al incluir nuevamente los términos difusivos el sistema no alcanza dicho estado estacionario estable, se dice entonces que el sistema de ecuaciones en (5) presenta inestabilidades por difusión o inestabilidades de Turing [14, 16]. El análisis matemático de la inestabilidad de Turing predice el tipo y forma de los patrones emergentes de la solución temporal de las ecuaciones de reacción difusión [2, 7, 18]. El primer paso del análisis es desacoplar el término de variación espacial para garantizar la estabilidad temporal. Luego se incorpora el término difusivo y se determina el espacio de parámetros que producen la inestabilidad espacial [1, 2, 7, 16]. Esta metodología aplica para un sistema de reacción difusión como el de la ecuación (5). En ausencia de difusión, el sistema en (5) queda reducido a (6):

]]> El estado uniforme de (6) es (u,v)=(u₀,v₀) tal que (7):

La inestabilidad de Turing tiene lugar cuando el estado estacionario (u₀,v₀) en (7) es estable en ausencia de difusión y se torna inestable en presencia del término difusivo. Se puede demostrar que para que un sistema de reacción difusión como el mostrado en (5) presente formación de patrones espacio-temporales se debe restringir el espacio de parámetros del sistema [9, 16]. Estas restricciones dan lugar al denominado espacio de Turing o espacio matemático al interior del cual el sistema exhibe la formación de patrones espacio- temporales. Este espacio está descrito por el siguiente conjunto de desigualdades (8) [2, 16]:

Los términos ƒ_u, ƒ_v, g_u y g_v en (8) representan las primeras derivadas de los términos de reacción respecto de las concentraciones u y v, evaluadas en el estado estacionario. D₁ y D₂ son los coeficientes de difusión de las sustancias u y v, y k es el autovalor de la solución de (5) usando el método de separación de variables [2, 16]. Una condición adicional permite aislar un patrón de formación espacial o modo de onda mediante la selección adecuada de los parámetros del modelo [18, 19]. Para la solución de un sistema de reacción difusión el modo de onda es el número de medio ciclos de onda coseno que describen los patrones espaciales en la dirección de cada uno de los ejes coordenados [16, 18]. Como se muestra más adelante, una adecuada selección de parámetros en cada modelo permite obtener diferentes modos de onda y, en general, diferentes patrones espacio-temporales.

Modelos biológicos

Entre la gran variedad de modelos biológicos basados en ecuaciones de reacción difusión [2, 3, 4, 6, 13, 15], existen tres modelos ampliamente documentados y analizados: el modelo de Schnakenberg, el modelo de Gierer-Meinhardt y el modelo de Thomas.

Modelo de Schnakenberg

El modelo de Schnakenberg es conocido por ser uno de los modelos de reacción difusión más sencillo y utilizado en morfogénesis. En su forma adimensional, el modelo está descrito por las siguientes ecuaciones (9) [2, 7, 15]:

]]>

Este modelo determina el comportamiento de un químico activador en presencia de un químico inhibidor. Si u es el químico activador, la reacción cinética es tal que en la ecuación (9a) el término u²v representa la producción de u en presencia de v, en tanto que en la ecuación (9b) el mismo término representa el consumo de v en presencia de u [18]. Las constantes a, b, d y Υ son todas parámetros positivos, donde a y b corresponden a valores de producción, Υ es una constante adimensional y d es el coeficiente de difusión [2, 15, 18].

Modelo de Gierer-Meinhardt

El modelo de Gierer-Meinhardt es un modelo fenomenológico de reacción cinética en el que una de las sustancias químicas, el activador, inicia la producción de la segunda sustancia, el inhibidor, que a su vez detiene la producción del activador. En su forma adimensional, el modelo está dado por las siguientes ecuaciones (10) [8, 9, 18]:

En (10) u es la sustancia activadora y v es la sustancia inhibidora. Las constantes a, b, d y Υ son todas parámetros positivos adimensionales y k es una medida de la concentración de saturación. En la ecuación (10a) el término u² / v (1 + ku²) representa tanto auto-catálisis en u con saturación para altos valores de concentración, como inhibición de u mediante la producción de v [9].

Modelo de Thomas

El modelo de Thomas se basa en la reacción de inhibición entre el oxígeno como sustrato y el ácido úrico en presencia de la enzima uricasa [4, 9]. En su forma adimensional, el modelo está descrito por las ecuaciones dadas en (11), donde u representa al oxígeno y v representa la uricasa.

En (11) las constantes a, b, Υ, α, y d son todos parámetros positivos. El término h(u,v) representa la inhibición por sustrato, inversamente proporcional al valor de u.

]]> En la siguiente sección se presenta la técnica de solución numérica utilizada para implementar los modelos anteriores mediante el método de los elementos finitos.

Metodología

Como los sistemas de ecuaciones de reacción difusión están expresados en términos de dos variables independientes, a saber, el tiempo (t) y la posición (x), su implementación numérica se suele realizar tratando independientemente estas dos variables [9, 18, 19]. En este trabajo se utiliza el método de elementos finitos para la discretización espacial y un esquema de diferencias finitas implícito del tipo Backward Euler para la discretización temporal. El sistema no lineal resultante, como consecuencia de la no linealidad de los términos reactivos, se resuelve usando el método de Newton-Raphson para la determinación de los ceros de una función [2, 20]. En seguida se describen en detalle las técnicas usadas.

Método de los elementos finitos

El empleo del método de los elementos finitos permite transformar el sistema de ecuaciones diferenciales parciales mencionado en (5), en un sistema de ecuaciones diferenciales ordinarias respecto del tiempo [21]. Su implementación se hace a partir de la forma débil del sistema, para lo cual inicialmente se realiza el producto interno de las ecuaciones por una función de prueba:

En (12) L₂ (Ω) es el conjunto de funciones cuadrado integrables, es decir:

Luego, se restringe el espacio de posibles soluciones de u al espacio de Sobolev H¹ (Ω) y se realiza el producto interno en H¹ (Ω), definido como (a, b)= ∫_Ω a •b • dΩ. Tomando inicialmente la ecuación (5a), a partir de (12) y (13) se obtiene (14) [20, 21]:

]]> Mediante el Teorema de Green y reemplazando las condiciones de frontera homogéneas, la expresión (14) puede expresarse como (15):

A partir de (15), el sistema de ecuaciones en (5) se discretiza usando el método de Bubnov-Galerkin, aproximando u y v como sigue (16) [2]:

En (16) N es la cantidad de nodos de la malla usada y las funciones u^h y w^h pertenecen a un espacio de dimensión finita H_h¹ (Ω), cuya base está formada por las funciones N_i, conocidas como funciones de forma [20]. Usando las expresiones en (16), la ecuación en (15) se transforma en (17):

Al aplicar las propiedades del producto interno, a partir de (17) se obtiene (18):

De esta manera, la ecuación en (18) es la transformación de (5a) de su forma continua en una expresión equivalente discreta en el espacio. El mismo procedimiento se sigue para la ecuación (5b). Adicionalmente, para obtener un sistema equivalente al expresado en (5) que pueda ser implementado y solucionado en un computador, es necesario complementar el procedimiento hasta ahora mencionado con la discretización de la segunda variable independiente, es decir, el tiempo.

Discretización en el tiempo

]]> Para discretizar en el dominio del tiempo la expresión (18) se utiliza el método de las diferencias finitas en atraso (Backward Euler). Asumiendo que u^k_i representa la concentración en el instante k para el nodo i, que u^h u_i^k+1 - u^k_i y que t (18) puede escribirse como (19):

Al reordenar los términos en (19) se obtiene (20):

La expresión (20) corresponde al equivalente discreto en tiempo de la expresión (18). De manera similar, el equivalente discreto en tiempo para la expresión (5b) es de la forma (21):

Método de Newton-Raphson

Siguiendo los procedimiento anteriores, el problema de expresar el sistema en (5) en su forma discreta, conforme a las expresiones (20) y (21), se transforma en la determinación de los ceros de las funciones F(u,v) y G(u,v) . Es decir, el problema se convierte en determinar los valores u_i ^k+1 y v_i ^{k +1} para los cuales las formas F y G del sistema en (5) son cero [20, 21].

El método de Newton-Raphson permite aproximar los valores buscados mediante un proceso iterativo que converge rápidamente (de forma cuadrática) a los valores de las raíces buscadas [2], a través del sistema lineal de ecuaciones dado por (22):

]]> Utilizando la forma matricial abreviada, el sistema de ecuaciones de la expresión (22) puede escribirse como (23):

En (23) A es la matriz de rigidez tangente, RHS es el vector de residuo, y u es el vector de incrementos. Este vector u corresponde además al vector de incógnitas que satisfacen las expresiones (20) y (21). La matriz A y el vector RHS son conocidos para un paso de tiempo i y resuelven calculando para cada paso de tiempo i+1. La solución del sistema (22) permite determinar de manera aproximada el valor de los incrementos u^h y v^h, que se considera es alcanzado cuando los elementos del vector RHS son menores que cierta tolerancia elegida (en nuestro caso 1x10^-6). Con el valor de los incrementos se calculan las concentraciones en el instante actual u_i ^k+1 = u^k_i + u^h y v_ik + v^h. A partir de estos valores se puede iniciar un nuevo proceso iterativo para calcular u_i ^k+2 y v_i ^k+2 [2, 20, 21].

Resultados

Los siguientes resultados corresponden a la implementación numérica de los modelos biológicos mencionados utilizando la técnica numérica propuesta en la sección anterior y una rutina de usuario programada en Fortran. El hardware utilizado es un PC de escritorio con procesador AMD Athlon 64 de 2.4 GHz y 1 GB de memoria RAM.

Caso unidimensional

Utilizando un segmento de longitud unitaria, se estudió la formación de patrones en una dimensión. En este caso, los tres modelos biológicos mencionados fueron implementados utilizando una segmentación del dominio en 300 elementos de tres nodos cada uno, para un total de 601 nodos. Las figuras 1, 2 y 3 muestran el estado estable en una dimensión para los tres modelos de reacción utilizados. Para cada modelo se presenta la formación de patrones en los cuatro primeros modos de onda disponibles, a saber, modo 1, modo 2, modo 3 y modo 4.

Caso bidimensional

Para el caso bidimensional se implementó el modelo de Schnakenberg en un dominio cuadrado de lado unitario. Se utilizaron mallas de 36, 144, 625, 2.500 y 10.000 elementos cuadráticos isoparamétricos. Los patrones obtenidos se muestran en las figuras 4, 5 y 6. Se presentan los resultados en una gráfica bidimensional donde los tonos más claros representan valores de concentración baja y los oscuros valores de concentración alta. Los medios tonos representan estados intermedios de concentración. Se analizan las distintas mallas y varios modos de oscilación en el espacio. En todos los casos la respuesta evoluciona gradualmente desde su condición inicial hacia los patrones de estado estable descritos por la teoría y los valores de Υ y d determinaron los modos de oscilación en el espacio como se referencia en la literatura [2, 19].

Tabla 1 Características físicas y químicas de los materiales utilizados

]]>

Figura 1 Patrones de estado estable de la concentración de u en una dimensión, usando el término de reacción de Schnakenberg con d=10. Las curvas de menor amplitud representan las condiciones iniciales. Las curvas restantes representan el valor de la concentración en el dominio considerado. a. Modo 1, Υ=35; b. Modo 2; Υ=130; c. Modo 3; Υ=250; d. Modo 4; Υ=465. Se usaron pasos de tiempo t =0,01 y se requirieron 100 iteraciones en todas las simulaciones

Figura 2 Patrones de estado estable de la concentración de u en una dimensión, usando el término de reacción de Thomas con d=30 Las curvas de menor amplitud representan las condiciones iniciales. Las curvas restantes representan el valor de la concentración en el dominio considerado. a. Υ=32; b. Υ=130; c. Υ=250; d. Υ=460. Se usaron pasos de tiempo t =0,01 y se requirieron 100 iteraciones en todas las simulaciones

Figura 3 Patrones de estado estable de la concentración de u en una dimensión, usando el término de reacción de Gierer-Meinhardt con d=72. Las curvas representan la concentración a lo largo del dominio considerado. a. Υ=70; b. Υ=270; c. Υ=650; d. Υ=1100. Se usaron pasos de tiempo t =0,01 y se requirieron 100 iteraciones en todas las simulaciones

Figura 4 Estado estable en un sistema de reacción difusión usando el término de reacción de Schnakenberg, en el modo (1,0) (Υ=29; d=10). Concentración de u. a. 36 elementos, b. 144 elementos, c. 625 elementos, d. 2.500 elementos. Se usaron pasos de tiempo t=0,01 y se requirieron 1.000 iteraciones en todas las simulaciones

Figura 5 Estado estable en un sistema de reacción difusión usando el término de reacción de Schnakenberg, en el modo (2,2) (Υ=230,82; d=8,6676). Concentración de u. a. 36 elementos, b. 144 elementos, c. 625 elementos, d. 2.500 elementos. Se usaron pasos de tiempo t=0,01 y se requirieron 1.000 iteraciones en todas las simulaciones

]]>

Figura 6 Estado estable en un sistema de reacción difusión usando el término de reacción de Schnakenberg, en el modo (3,3) (Υ=435,99; d=8,6076). Concentración de u. a. 36 elementos, b. 144 elementos, c. 625 elementos, d. 2.500 elementos. Se usaron pasos de tiempo t=0,01 y 1.000 en todas las simulaciones

Discusión

La implementación numérica de tres modelos biológicos conocidos descritos por ecuaciones de reacción difusión para el caso unidimensional y bidimensional permite verificar la formación de patrones espaciales propia de cada modelo. La metodología utilizada se presenta como una herramienta flexible para la implementación de diferentes sistemas de reacción difusión en dominios fijos en una, dos e incluso tres dimensiones [18, 19, 20]. Para el caso unidimensional y bidimensional se observa que la formación de patrones espacio-temporales esta relacionada con los valores de los parámetros del modelo (figuras 1-6). Estos parámetros, a su vez, pertenecen al espacio de Turing al interior del cual se cumplen las condiciones dadas por (8) y se garantiza tanto la inestabilidad espacial como la estabilidad temporal [7, 8]. Un análisis detallado del espacio de Turing permite establecer rangos específicos para los parámetros que condicionan tanto la aparición como la evolución de la característica espaciotemporal de la solución [16].

Los resultados obtenidos concuerdan con los resultados reportados por otros autores [2, 18]. Biológicamente, los patrones espacio-temporales encontrados son similares a los patrones de pigmentación en la piel de algunos animales [1, 15], al comportamiento de activación e inhibición que gobierna la interacción entre diferentes sustancias [4, 7] y a la aparición de zonas de formación y crecimiento de organismos vivos [2, 10], entre otros. Se determinó además que el estado estable es independiente de las condiciones iniciales, incluso cuando éstas se encuentran alejadas del estado estacionario. En el caso bidimensional se observa que para cada modo de onda se obtienen resultados cualitativamente similares entre las diferentes mallas usadas, siempre que la geometría del patrón espacial sea mayor que las dimensiones de uno de los elementos de la malla (véase los puntos oscuros de las figuras 5 y 6). Cuando las dimensiones del elemento son mayores a las del patrón espacial esperado se obtienen patrones erróneos. La figura 7 muestra el caso del modo de onda (4,4), resuelto usando mallas de 36, 144, 625 y 2.500 elementos. En los dos primeros casos, se observa que el patrón no coincide con el esperado modo de onda (4,4), esto es, un total de cuatro cambios en el valor de la amplitud del patrón en cada una de las direcciónes x y y. Este hecho pone en evidencia el problema de distorsión del patrón espacial por deficiencias del mallado. En general, el uso de mallas con mayor número de elementos permite aumentar la exactitud de la solución numérica y a su vez obtener un mejor ajuste del patrón buscado [2]. No obstante, la semejanza entre los dos últimos casos permite concluir que un aumento exagerado del mallado (2.500 elementos), desde el punto de vista cualitativo, puede no representar una mejora sustancial en la solución, siempre que el mallado con 650 elementos ya permite obtener el patrón espacio-temporal esperado [21].

Figura 7 Estado estable en un sistema de reacción difusión usando el término de reacción de Schnakenberg, en el modo (4,4) (Υ=909,66; d= 8,6076). Concentración de u. a. 36 elementos, b. 144 elementos, c. 625 elementos, d. 2.500 elementos. Se usaron pasos de tiempo t=0,01 y 1.000 en todas las simulaciones

Los patrones bidimensionales no en todos los casos se presentan con la misma orientación para todas las mallas usadas (figura 4). Con la metodología usada y por la simetría de las soluciones analíticas para el sistema linealizado en las direcciones x y y, no se puede garantizar la aparición de las oscilaciones en direcciones específicas, aunque es posible predecir la presencia de oscilaciones en dirección horizontal y en dirección vertical. La técnica utilizada permite incluir otros términos en el modelo que hacen referencia, por ejemplo, a las deformaciones del dominio debido a fuerzas externas. Trabajos futuros podrían extenderse al estudio de la formación de patrones en dominios crecientes, como en el caso de un organismo vivo donde la formación de patrones se puede presentar durante la etapa de crecimiento [2, 3].

La implementación numérica utilizada proporciona una solución a problemas complejos que requiere de menor costo computacional y permite obtener una mejor aproximación numérica, siempre que el dominio, el mallado y las características temporales sean bien especificados [2]. De esta manera, se espera que la técnica de solución empleada sea de utilidad en la formulación e implementación de modelos matemáticos biológicos complejos que describan fenómenos de crecimiento y desarrollo celular y tisular.

Referencias

]]>

1. P. K. Maini."Mathematical models in morphogenesis". Mathematics Inspired by Biology. Ed. Springer Berlin- Heidelberg. 1999. pp. 151-189.        [ Links ]

2. D. A. Garzón. Simulación de procesos de reacción difusión: aplicación a la morfogénesis del tejido óseo. Tesis de Doctorado. Zaragoza. España: Centro Politécnico Superior de la Universidad de Zaragoza. 2007. pp. 38-70.        [ Links ]

3. D. Benson, P. Maini, J. A. Sherratt."Unravelling the Turing bifurcation using spatially varying diffusion coefficient". J. Math. Biol. Vol. 37. 1998. pp. 381- 417.        [ Links ]

4. P. K. Maini, K. Painter, H. Chau."Spatial pattern formation in chemical and biological systems". J. Chem. Soc. Faraday Trans. Vol. 93. 1997. pp. 3601- 3610.        [ Links ]

5. A. Turing."The chemical basis of morphogenesis". Philos. Trans. Roy. Soc. Vol. 237. 1952. pp. 37-72.        [ Links ]

6. M. Dolnik, A. M. Zhabotinsky, A. B. Rovinsky, I. R. Epstein."Spatio-temporal patterns in a reaction diffusion system with wave instability". Chemical Engineering Science. Vol. 55. 2000. pp. 223-231.        [ Links ]

7. J. D. Murray. Mathematical Biology II. Spatial models and biomedical applications. 3^a. ed. Ed. Springer- Verlag. New York. 2003. pp. 75-97.        [ Links ]

8. K. Page, P. K. Maini, N. Monk."Patter formation in spatially heterogeneous Turing reaction diffusion models". Physics D. Vol. 181. 2003. pp. 80-101         [ Links ]

9. A. Madzvamuse, A. J. Wathen, P. K. Maini."A moving grid finite element method applied to a model biological patter generador". Journal of Computational Physics. Vol. 190. 2003. pp. 478-500.        [ Links ]

10. T. Miura, K. Shiota."TGFß-2 acts as an activator molecule in reaction-diffusion model and is involved in cell sorting phenomenon in mouse limb micromass culture". Dev. Dyn. Vol. 217. 2000. pp. 241-249.        [ Links ]

11. S. Kondo, R. Asai."A reaction diffusion wave on the skin of the marine angelfish Pomacanthus". Nature. Vol. 376. 1995. pp. 765-768.        [ Links ]

12. H. Jung, P. Francis-West, R. Widelitz, T. Jiang, S. Ting- Berreth, C. Tickle, L. Wolpert, C. Chuong."Local inhibitory action of BMPs and their relationships with activators in feather formation: implications for periodic patterning". Dev. Biol. Vol. 196. 1998. pp. 11-23.        [ Links ]

13. T. Miura, K. Morriss, P. K. Maini."Mixed-mode pattern in doublefoot mutant mouse limb: Turing reaction diffusion model on a growing domain during limb development". J. Theor Biol. Vol. 240. 2006. pp. 562-573.        [ Links ]

14. J. D. Murray. Mathematical Biology I. An introduction.3^a ed. Ed. Springer-Verlag. New York. 2002. pp. 405- 509.        [ Links ]

15. K. Painter. Chemotaxis as a mechanism for morphogenesis. PhD Thesis. United Kingdom. Oxford University. 1997. pp. 3-25.        [ Links ]

16. J. C. Vanegas A., N. S. Landinez, D. A. Garzón- Alvarado."Análisis de la inestabilidad de Turing en modelos biológicos". Revista DYNA, Universidad Nacional de Colombia - Medellín. Aprobado para publicación. 2009.        [ Links ]

17. J. D. Murray, G. F. Ester."Cell traction models for generation pattern and form in morphogenesis". J. Math. Biology. Vol. 19. 1984. pp. 265-279.        [ Links ]

18. A. Madzvamuse. A numerical approach to the study of spatial pattern formation. PhD Thesis. Oxford, UK: Computing Laboratory. University of Oxford. 2000. pp. 10-40.        [ Links ]

19. E. Crampin. Reaction diffusion patterns on growing domains. PhD Thesis. Oxford, UK: Magdalen College. University of Oxford. 2000. pp.1-34.        [ Links ]

20. E. Oñate, J. Miquel, F. Zárate."Stabilized solution of the multidimensional advection-diffusion-absorption equation using linear finite elements". Computers and Fluids. Vol. 36. 2007. pp. 92-112.        [ Links ]

21. A. Madzvamuse."Time-stepping schemes for moving grid finite elements applied to reaction diffusion systems on fixed and growing domains". J. Comput Phys. Vol. 214. 2006. pp. 239-263.        [ Links ]

(Recibido el 3 de junio de 2008. Aceptado el 12 de marzo de 2009)

^*Autor de correspondencia: teléfono: +57 +1 +316 53 39, fax: +57 +1 +316 53 43, correo electrónico: jcvanegasa@unal.edu.co (J. Vanegas).

]]> 1 1999 151-189 2 38-70 3 1998 37 381- 417 4 1997 93 3601- 3610 5 1952 237 37-72 6 2000 55 223-231 7 2003 3 75-97 8 2003 181 80-101 9 2003 190 478-500 10 2000 217 241-249 11 1995 376 765-768 12 1998 196 11-23 13 2006 240 562-573 14 2002 3 405- 509 15 3-25 16 2009 17 1984 19 265-279 18 10-40 19 1-34 20 2007 36 92-112 21 2006 214 239-263