Turing patterns formation on surfaces under deformation: A total lagrangian method approach

Resumen

En este artículo se desarrollan varios ejemplos numéricos sobre ecuaciones de reacción-difusión con dominio creciente. Para este fin se utiliza el modelo de reacción de Schnakenberg, con parámetros en el espacio de Turing. Por tanto se realizan ensayos numéricos sobre la aparición de los patrones de Turing en superficies que tienen alta tasa de deformación. Para la solución de las ecuaciones de reacción difusión se presenta un método de solución en superficies en 3 dimensiones mediante el método de los elementos finitos bajo el uso de la formulación lagrangiana total. Los resultados muestran que la formación de los patrones de Turing depende de las funciones de deformación de la superficie y la tasa a la cual se presenta el cambio de posición de cada punto del dominio donde se lleva a cabo la solución numérica. Estos resultados pueden esclarecer algunos fenómenos de cambio de patrón en la superficie de la piel de aquellos animales que exhiben manchas características.

Abstract

In this work we have developed several numerical examples of reaction-diffusion equations with growing domain. For this purpose we have used the Schnakenberg reaction model with parameters in space Turing. Therefore numerical tests are performed on the appearance of Turing patterns on surfaces that have high deformation rate. For the solution of reaction diffusion equations is presented a solution method on surfaces in three dimensions using the finite element method under the use of the total Lagrangian formulation. The results show that the formation of Turing patterns depends on the features of surface deformation and the rate at which change in position of each point of the domain. These results can explain some phenomena of change of pattern on the surface of the skin of animals that exhibit characteristic spots.

Palabras clave

Reacción-Difusión ; Turing ; Lagrangiano total ; Elementos finitos ; Deformación de superficies

Keywords

Reaction-Diffusion ; Turing ; The total Lagrangian formulation ; Finite elements ; Deformation of surfaces

1. Introducción

Las ecuaciones de reacción-advección-difusión [1] , [2] , [3] , [4] , [5] , [6] , [7] , [8] , [9] , [10] , [11] , [12]  and [13] (RAD) (1) y otros modelos más complejos, donde intervienen más especies o reactantes, poseen la habilidad para crear patrones espacio-temporales. Un caso particular de estos patrones son las inestabilidades de Turing[14]  and [15] , las cuales se caracterizan por la aparición de distribuciones de especies (patrones) estables en el tiempo e inestables en el espacio. Dada su especial respuesta, este tipo de modelos matemáticos han inspirado modelos para el estudio de problemas en muy diversos campos, tales como dinámica de fluidos [16] , transferencia de calor [17]  and [18] , física de semiconductores [19] , ingeniería de materiales [20] , química [21] , biología [7] , [8] , [14]  and [22] , dinámica de poblaciones [9] , [23]  and [24] , astrofísica [25] , ingeniería biomédica [3] , [5]  and [26] y matemáticas financieras entre otros.

El análisis de estos sistemas de reacción-difusión (RD), que presentan inestabilidad de Turing, se ha desarrollado tradicionalmente desde 2 marcos de trabajo: mediante análisis matemático [27] y mediante simulación numérica [14] , [15]  and [28] .

Desde el punto de vista analítico, los esfuerzos por entender el comportamiento de los sistemas de RD se han centrado en el estudio de la relación entre las bifurcaciones del espacio de parámetros y la formación de patrones. Desde esta perspectiva, se han estudiado los sistemas de RD mediante comparaciones de sub y super soluciones, teoría de grado, índice de Conley, teoría de puntos críticos y perturbaciones singulares para varios tipos de máximos principales [27] . Estos métodos han sido efectivos para el análisis de soluciones estacionarias y ondas viajeras [29]  and [30] . También se han estudiado escenarios de bifurcaciones complejas en sistemas de RD aplicando métodos de teoría de grupo para problemas con simetrías [31]  and [32] . Los esfuerzos en esta área del análisis matemático y, específicamente, de la dinámica de sistemas han permitido construir un gran conocimiento, el cual se ha comprobado y ampliado con el uso de la simulación numérica.

La simulación numérica de sistemas de RD ha permitido corroborar el conocimiento obtenido analíticamente sobre la formación de patrones, como en: Madzvamuse, [15]  and [28] y Painter et al. [32]  and [33] , donde se han desarrollado ejemplos numéricos sobre la formación de patrones en dominios bidimensionales bajo consideraciones de dominio creciente. En [28] se reporta la aparición de diferentes estructuras que pueden variar entre sistemas de bandas, puntos y combinaciones de estos patrones en dominios con decrecimiento exponencial. En [34] se reporta la formación de patrones en presencia de campos convectivos con divergencia nula. Es así como se ha dedicado un gran esfuerzo al estudio de los patrones de Turing bajo la deformación del dominio de la solución y cuando esta sometido a campos de flujo. Es importante notar que el trabajo de Madzvamuse et al. [5] , [15] , [22] , [28]  and [35] se ha centrado en el estudio del efecto que tiene el crecimiento uniforme sobre la dinámica de los patrones de Turing. Adicionalmente, importantes trabajos y aplicaciones se han desarrollado en el área de morfogénesis y movimiento celular [36] . A diferencia de estos trabajos, en este artículo se hace especial énfasis en el crecimiento no uniforme del dominio, con el objeto de estudiar el efecto de los gradientes de desplazamiento en la aparición y formación de los patrones de Turing. Además, en este trabajo se hacen ensayos numéricos con crecimiento perpendicular al plano, por lo que se puede controlar la evolución de cada punto del dominio.

Siguiendo un planteamiento similar al utilizado en [35]  and [36] , se utiliza el modelo de reacción de Schnakenberg [14] con parámetros en el espacio de Turing para simular la aparición de los patrones en superficies que tienen alta tasa de deformación. Para la solución de las ecuaciones de reacción difusión se presenta un método de solución en superficies en 3 dimensiones mediante el método de los elementos finitos bajo el uso de la formulación lagrangiana total.

El artículo está organizado como sigue: en primer lugar se muestran las ecuaciones de reacción difusión y las condiciones necesarias para generar patrones de Turing, luego se muestra el procedimiento de solución de las ecuaciones no lineales resultantes. Posteriormente se presentan los ejemplos numéricos para diferentes tipos de patrones. Por último se presenta la discusión y conclusiones.

2. Materiales y métodos

2.1. Sistema de reacción-difusión

Un sistema de reacción difusión (RD) puede controlar la formación de los patrones siempre que sus parámetros reactivos y difusivos se encuentren en el espacio de Turing [14] . Para estudiar los patrones de Turing, se inicia con la definición de un sistema de RD, para 2 especies, dado por (1):

${\begin{array}{c}{\frac {\partial u_{1}}{\partial t}}-{\nabla }^{2}u_{1}=\gamma \cdot f(u_{1},u_{2})\\{\frac {\partial u_{2}}{\partial t}}-d{\nabla }^{2}u_{2}=\gamma \cdot g(u_{1},u_{2})\end{array}}$

( 1)

donde u₁y u₂ son las concentraciones, en su forma adimensional, de las especies químicas presentes las ecuaciones (1), d es el coeficiente de difusión adimensional y ( es una constante de adimensionalización. Por su parte, f y g , son las funciones de reacción que se presentan en cada una de las ecuaciones (1). Estas funciones tienen, en este artículo, un carácter no lineal y tienen una participación importante en la formación de los patrones de Turing.

Turing [37] en su libro «The chemical basis of morphogenesis» desarrolló las condiciones necesarias para la formación de patrones espaciales a partir de la ecuación (1). Las condiciones para la formación de patrones determinan el espacio de Turing dado por las siguientes restricciones (2):

${\begin{array}{c}f_{u_{1}}g_{u_{2}}-f_{u_{2}}g_{u_{1}}>0\\f_{u_{1}}+g_{u_{2}}<0\\df_{u_{1}}+g_{u_{2}}>0\\{\left(df_{u_{1}}+g_{u_{2}}\right)}^{2}>4d\left(f_{u_{1}}g_{u_{2}}-f_{u_{2}}g_{u_{1}}\right)\end{array}}$

( 2)

donde f_l y g_l indican las derivadas de las funciones de reacción con respecto a las variables de concentración, por ejemplo ${\textstyle f_{u}={\frac {\partial f}{\partial u}}}$[14] . Estas condiciones (2) están evaluadas en el punto de equilibrio que se obtiene haciendo ${\textstyle f{\left(u_{1},u_{2}\right)}_{0}=g{\left(u_{1},u_{2}\right)}_{0}=}$$0$ .

Las ecuaciones (1) y sus restricciones (2) permitieron el desarrollo de una rama de investigación de los sistemas dinámicos [14] : las inestabilidades de Turing. La teoría acerca de los patrones de Turing ha permitido explicar la formación de patrones biológicamente complejos, como las manchas que se encuentran en la piel de algunos animales [38]  and [39] y en problemas de morfogénesis [40] , entre otros. Además, recientemente se ha comprobado, experimentalmente, que el comportamiento de algunos sistemas RD generan patrones de ondas viajeras y patrones espaciales estables [41] , [42]  and [43] .

Las ecuaciones (1) contienen términos reactivos (f (u₁ , u₂ ) y  g (u₁ , u₂ )) que tienen una importante contribución en la formación de los patrones de Turing. En este artículo, las ecuaciones que se han utilizado para el término reactivo son las de Schnakenberg [14]  and [28] dadas por:

${\begin{array}{c}f(u_{1},u_{2})=a-u_{1}+{\left(u_{1}\right)}^{2}u_{2}\\g(u_{1},u_{2})=b-{\left(u_{1}\right)}^{2}u_{2}\end{array}}$

( 3)

Donde a y b   son parámetros adimensionales del modelo. Los puntos de estado estable están dados por ${\textstyle {\left(u_{1},u_{2}\right)}_{0}=\left(a+\right.}$$\left.b,{\frac {b}{{\left(a+b\right)}^{2}}}\right)$ . Aplicando las restricciones (2) al modelo (1) en el punto de estado estable ${\textstyle {\left(u_{1},u_{2}\right)}_{0}}$ se obtiene un conjunto de restricciones que permiten establecer el sitio geométrico de los parámetros en el espacio de Turing [14]  and [28] .

Haciendo uso del análisis de estabilidad lineal [14] , [15]  and [28] , se pueden calcular los valores de los parámetros adimensionales d y γ , requeridos para la formación de patrones de Turing bajo el modelo de RD. El cálculo de estos parámetros se hace mediante el uso de las ecuaciones (2). En primer lugar se fijan los valores a =0,1 y b =0,9; luego, con las desigualdades descritas en (2) se obtienen intervalos para d y γ . Por último, en estos intervalos se eligen valores particulares, con el objeto de cumplir los requisitos de formación de patrones específicos, como se describe en la tabla 1 . Estos patrones pueden ser identificados por los correspondientes números de onda, tal como los que se mencionan en la tabla 1 .

Tabla 1. Modos de vibración (m , n ) para diferentes valores de d y γ del modelo de Schnakenberg con valores constantes a  = 0, 1 y b  = 0, 9 [28]

Modo (m,n)	d	γ
(1,0)	10,0000	29,0000
(1,1)	11,5776	70,6000
(2,0)	10,0000	114,0000
(2,1)	9,1676	176,7200
(2,2)	8,6676	230,8200
(3,0)	8,6176	265,2200
(3,1)	8,6676	329,2000
(3,2)	8,8676	379,2100
(3,3)	8,6076	535,0900
(4,0)	8,6676	435,9900
(4,1)	8,5876	492,2800
(4,2)	8,7176	625,3500
(4,3)	8,6676	666,8200
(4,4)	8,6076	909,6600

2.2. Deformación de las superficies

Para deformar la superficie donde se lleva a cabo se solución de las ecuaciones (1) se utilizaron las ideas sugeridas en [44] , [45]  and [46] . Con el objeto de estudiar el efecto de la deformación no uniforme del dominio, se desplaza cada punto, con diferente velocidad, en dirección perpendicular al plano xy . Por tanto, la ecuación de desplazamiento normal esta dada por (4):

${\frac {dz}{dt}}=h(x,y,t)$

( 4)

Donde h (x , y , t   ) es una función que determina la velocidad de desplazamiento. De esta forma, la velocidad de cada punto de la superficie esta dada por ${\textstyle {\boldsymbol {\mbox{v}}}=h(x,y,t){\overset {\mbox{ˆ}}{k}}}$ , donde ${\textstyle {\overset {\mbox{ˆ}}{k}}}$ es el vector unitario en dirección z (perpendicular al plano).

Al incluirse el término de crecimiento de la superficie [ecuación (4)], se modifican las ecuaciones (1), donde se presenta un nuevo término que tiene en cuenta la convección y la dilatación del dominio dado por:

${\begin{array}{c}{\frac {\partial u_{1}}{\partial t}}+div(u_{1}{\boldsymbol {\mbox{v}}})-{\nabla }^{2}u_{1}=\gamma \cdot f(u_{1},u_{2})\\{\frac {\partial u_{2}}{\partial t}}+div(u_{2}{\boldsymbol {\mbox{v}}})-d{\nabla }^{2}u_{2}=\gamma \cdot g(u_{1},u_{2})\end{array}}$

( 5)

Donde el término nuevo ${\textstyle div(u_{i}{\boldsymbol {\mbox{v}}})}$ incluye la convección y dilatación que se debe al crecimiento del dominio.

Para solucionar el RD convección descrito en (5) se utiliza el método de los elementos finitos [47] y el método de Newton-Raphson [48] para solucionar el sistema no lineal de ecuaciones en derivadas parciales que se derive de la formulación. La imposición del campo de crecimiento sobre la superficie se hace mediante la solución de la ecuación (4), con lo que se obtiene la nueva configuración (actual) y el campo de velocidades que se incluirá en el problema de RD.

En primer lugar se muestra la solución de las ecuaciones de reacción difusión mediante el método de los elementos finitos.

Para calcular el movimiento de la malla y la velocidad a la cual se deforma el dominio se utiliza la ecuación (4), la cual se integra mediante el método de Euler, dado por [26] :

$z_{t+dt}=z_{t}+h(x,y,t)dt$

( 6)

Donde z_{t +dt} y z_t son la configuración de la superficie en el estado t y t+dt. Por tanto la velocidad se obtiene mediante (7):

${\boldsymbol {\mbox{v}}}={\frac {z_{t+dt}-z_{t}}{dt}}{\overset {\mbox{ˆ}}{k}}$

( 7)

Donde el término de velocidad tiene dirección y magnitud que depende del punto material de la superficie S .

2.3. Solución del sistema de reacción-convección-difusión

La formulación del sistema RD incluyendo el transporte convectivo se puede escribir como (8) [14] :

${\begin{array}{c}{\frac {\partial u_{1}}{\partial t}}+div(u_{1}{\boldsymbol {\mbox{v}}})={\nabla }^{2}u_{1}+\gamma \cdot f(u_{1},u_{2})\\{\frac {\partial u_{2}}{\partial t}}+div(u_{2}{\boldsymbol {\mbox{v}}})=d{\nabla }^{2}u_{2}+\gamma \cdot g(u_{1},u_{2})\end{array}}$

( 8)

Donde u₁y u₂ son las variables químicas del RD. Esta ecuación también se puede escribir en términos de la derivada total (9) [14] :

${\begin{array}{c}{\frac {du_{1}}{dt}}+u_{1}div({\boldsymbol {\mbox{v}}})={\nabla }^{2}u_{1}+\gamma \cdot f(u_{1},u_{2})\\{\frac {du_{2}}{dt}}+u_{2}div({\boldsymbol {\mbox{v}}})=d{\nabla }^{2}u_{2}+\gamma \cdot g(u_{1},u_{2})\end{array}}$

( 9)

Donde se debe tener en cuenta que ${\textstyle {\frac {du}{dt}}={\frac {\partial u}{\partial t}}+}$${\boldsymbol {\mbox{v}}}\cdot grad(u)$[49]  and [50] .

Según como se describe en [49] se llega a que, el RD convección en la configuración inicial, o de referencia Ω₀ [con coordenadas materiales X(x) ], está dada por la siguiente ecuación escrita en términos de las coordenadas materiales y con notación indicial [49] :

${\frac {dU_{1}}{dt}}+U_{1}{\frac {\partial v^{i}}{\partial x^{i}}}=$$\gamma F(U_{1},U_{2})+{\left(F^{-1}\right)}_{i}^{I}{\frac {\partial }{\partial X^{I}}}\left({\frac {\partial u_{1}}{\partial x^{j}}}\right){\delta }^{ij}$

( 10a)

${\frac {dU_{2}}{dt}}+U_{2}{\frac {\partial v^{i}}{\partial x^{i}}}=$$\gamma G(U_{1},U_{2})+d{\left(F^{-1}\right)}_{i}^{I}{\frac {\partial }{\partial X^{I}}}\left({\frac {\partial u_{2}}{\partial x^{j}}}\right){\delta }^{ij}$

( 10b)

Donde U₁y U₂ son las concentraciones de cada una de las especies en la configuración inicial Ω₀ , esto es, U (X , t ) = u (x (X , t ), t   ). Además, ${\textstyle {\left(F^{-1}\right)}_{I}^{i}}$ es el inverso del gradiente de deformación dado por ${\textstyle F_{I}^{i}={\frac {\partial x^{i}}{\partial X^{I}}}}$[3] , xⁱ son las coordenadas actuales (en cada instante de tiempo) y X^I son las coordenadas iniciales (de referencia, donde se realizarán los cálculos) [49]  and [50] .

Multiplicando la ecuación (10) por una función de ponderación W e integrando en el dominio Ω₀ se obtiene la forma débil general (11) [47] :

$\int _{{\Omega }_{0}}W\left({\frac {d\left(U\right)}{dt}}J+\right.$$\left.Jdiv(v)U-\gamma F(U,V)J\right)d{\Omega }_{0}+$$\int _{{\Omega }_{0}}{\frac {\partial W}{\partial X^{I}}}J{\underset {{\left(C^{-1}\right)}^{IJ}}{\underbrace {{\delta }^{ij}{\left(F^{-1}\right)}_{i}^{I}{\left(F^{-1}\right)}_{j}^{J}} }}{\frac {\partial U}{\partial X^{J}}}d{\Omega }_{0}=$$0$

( 11)

Donde U es cualquiera de las 2 especies en estudio (U₁ o U₂ ), W es la función de peso (o ponderación), J es el jacobiano (y es igual al determinante del gradiente de deformación F ) y C^-1 es el inverso del tensor Cauchy – Green por la derecha [47]  and [48] .

En el caso de la formulación lagrangiana total, el cálculo se realiza, siempre, en la configuración inicial de referencia Ω₀ . La integración temporal se lleva a cabo mediante el método de Backward-Euler, el cual es estable para problemas parabólicos [47] . En este caso, todas las variables están calculadas en el nuevo tiempo. Además, gracias a la no linealidad se utiliza el método de Newton Raphson para la solución numérica. Por tanto, la solución del sistema (10a) y (10b) se inicia con la discretización de las variables U₁ y U₂ (calculadas en el tiempo t) mediante (12) [47] :

$U_{l}^{h}={\boldsymbol {\mbox{N}}}_{U}(X,Y){\boldsymbol {\mbox{U}}}_{1}=$$\sum _{p=1}^{nnod}N_{p}U_{1p}$

( 12a)

$U_{2}^{h}={\boldsymbol {\mbox{N}}}_{V}(X,Y){\boldsymbol {\mbox{U}}}_{2}=$$\sum _{p=1}^{nnod}N_{p}U_{2p}$

( 12b)

donde nnod es el número de nodos, U₁ y U₂ son los vectores que contiene los valores de U₁y U₂ en los puntos nodales y el superíndice h indica la discretización de la variable en elementos finitos. Mediante la elección de las funciones de peso igual a las funciones de forma (Galerkin Estándar) se obtienen los vectores residuo del método de Newton-Raphson dado por [47] (13):

$r_{U_{1}}^{h}=\int _{{\Omega }_{0}}N_{p}J{\frac {dU_{1}^{h}}{dt}}d{\Omega }_{0}+$$\int _{{\Omega }_{0}}N_{p}{\frac {dJ}{dt}}UU_{1}^{h}d{\Omega }_{0}-$${\int }_{{\Omega }_{0}}N_{p}J\gamma F\left(U_{1},U_{2}\right)d{\Omega }_{0}+$${\int }_{{\Omega }_{0}}{\frac {\partial N_{p}}{\partial X^{I}}}J{\left(C^{-1}\right)}^{IJ}{\frac {\partial U_{1}}{\partial X^{J}}}d{\Omega }_{0}$

( 13a)

${\boldsymbol {\mbox{r}}}_{{\boldsymbol {\mbox{U}}}_{2}}^{\boldsymbol {h}}=$${\int }_{{\Omega }_{0}}N_{p}J{\frac {dU_{2}^{h}}{dt}}d{\Omega }_{0}+$${\int }_{{\Omega }_{0}}N_{p}{\frac {dJ}{dt}}U_{2}^{h}d{\Omega }_{0}-$${\int }_{{\Omega }_{0}}N_{p}J\gamma G(U_{1},U_{2})d{\Omega }_{0}+$${\int }_{{\Omega }_{0}}{\frac {\partial N_{p}}{\partial X^{I}}}dJ{\left(C^{-1}\right)}^{IJ}{\frac {\partial U_{2}}{\partial X^{J}}}d{\Omega }_{0}$

( 13b)

Con p=1,…, nnod   . Donde ${\textstyle {\boldsymbol {r}}_{{\boldsymbol {U}}_{1}}^{\boldsymbol {h}}}$ y ${\textstyle {\boldsymbol {r}}_{{\boldsymbol {U}}_{\boldsymbol {2}}}^{\boldsymbol {h}}}$ son los vectores residuos que están calculados en el nuevo tiempo. Se debe notar que, ${\textstyle {\frac {dU_{1}^{h}}{dt}}={\frac {U1h\left(t\right)-U1h\left(t-\Delta t\right)}{\Delta t}},{\mbox{ }}{\frac {dU_{2}^{h}}{dt}}=}$${\frac {U2h\left(t\right)-U2h\left(t-\Delta t\right)}{\Delta t}},{\mbox{ }}U_{1}^{h}=$$U1h\left(t\right){\mbox{ }}y{\mbox{ }}U_{2}^{h}=U2h\left(t\right){\mbox{,}}$ donde se han omitido los superíndices izquierdos en las ecuaciones (12) y (13) para facilidad de lectura. Por su parte, cada una de las posiciones (entradas) de la matriz jacobiana está dada por (14):

${\frac {\partial {\boldsymbol {\mbox{r}}}_{U_{1}}^{h}}{\partial {\boldsymbol {\mbox{U}}}_{1}}}=$${\frac {1}{\Delta t}}{\int }_{{\Omega }_{0}}N_{p}JN_{s}d{\Omega }_{0}+$${\int }_{{\Omega }_{0}}N_{p}{\frac {dJ}{dt}}N_{s}d{\Omega }_{0}-$${\int }_{{\Omega }_{0}}N_{p}\gamma J{\frac {\partial F(U_{1},U_{2})}{\partial U_{1}}}N_{s}d{\Omega }_{0}+$${\int }_{{\Omega }_{0}}{\frac {\partial N_{p}}{\partial X^{I}}}J{\left(C^{-1}\right)}^{IJ}{\frac {\partial N_{s}}{\partial X^{J}}}d{\Omega }_{0}$

( 14a)

${\frac {\partial {\boldsymbol {\mbox{r}}}_{U_{1}}^{h}}{\partial {\boldsymbol {\mbox{U}}}_{2}}}=$$-{\int }_{{\Omega }_{0}}N_{p}\gamma J{\frac {\partial F(U_{1},U_{2})}{\partial U_{2}}}N_{s}d{\Omega }_{0}$

( 14b)

${\frac {\partial {\boldsymbol {\mbox{r}}}_{U_{2}}^{h}}{\partial {\boldsymbol {\mbox{U}}}_{1}}}=$$-{\int }_{{\Omega }_{0}}N_{p}\gamma J{\frac {\partial G(U_{1},U_{2})}{\partial U_{1}}}N_{s}d{\Omega }_{0}$

( 14c)

${\frac {\partial {\underline {{\boldsymbol {\mbox{r}}}_{U_{2}}^{h}}}}{\partial {\boldsymbol {\mbox{U}}}_{2}}}=$${\frac {1}{\Delta t}}{\int }_{{\Omega }_{0}}N_{p}JN_{s}d{\Omega }_{0}+$${\int }_{{\Omega }_{0}}N_{p}{\frac {dJ}{dt}}N_{s}d{\Omega }_{0}-$${\int }_{{\Omega }_{0}}N_{p}\gamma J{\frac {\partial G(U_{1},U_{2})}{\partial U_{2}}}N_{s}d{\Omega }_{0}+$${\int }_{{\Omega }_{0}}{\frac {\partial N_{p}}{\partial X^{I}}}dJ{\left(C^{-1}\right)}^{IJ}{\frac {\partial N_{s}}{\partial X^{J}}}d{\Omega }_{0}$

( 14d)

Donde J es el determinante del gradiente de deformación, C^-1 es el inverso del tensor de Cauchy-Green por la derecha y p,s =1,…, nnod e I,J=1 ,…,dim , donde dim representa la dimensión en la cual se resuelve el problema. Por tanto, utilizando las ecuaciones (13) y (14) se puede implementar el método de Newton-Raphson para solucionar el sistema de reacción-difusión utilizando su descripción material. En este sentido se debe anotar que la integración de (13) y (14) se hace sobre la configuración inicial [49] . Una deducción detallada de las ecuaciones se puede encontrar en [49] .

2.4. Aspectos de la implementación computacional

Para implementar el modelo de reacción-difusión en elementos finitos se utiliza la formulación antes descrita. Es importante anotar que, aunque la superficie se encuentra orientada en el espacio 3D, los cálculos numéricos se hacen en 2D. Para este objetivo, se encuentra la normal a cada elemento (Z’) y se ubican los ejes primos (X’Y’) formando un plano paralelo al plano del elemento (fig. 1 ). Para enmallar la geometría se utilizan elementos triangulares de primer orden con 3 nodos. Por tanto, el cálculo se simplifica pasando de un sistema 3D a un sistema bidimensional donde se soluciona los modelos de reacción-difusión en cada instante de tiempo. La relación existente entre los ejes X’Y’Z’ y XYZ se puede obtener mediante una matriz de transformación T [49] .

Figura 1. Orientación espacial de los ejes coordenados.

Para solucionar el sistema de ecuaciones resultantes del método de los elementos finitos con el método de Newton-Raphson se hizo un programa en FORTRAN y se solucionaron los siguientes ejemplos en un Laptop de 4096 MB en RAM y 800 MHz de velocidad de procesador.

Para llevar a cabo la deformación del dominio se han utilizado las ecuaciones que se describen a continuación (fig. 2 ):

Figura 2. a) Malla de triángulos utilizada en las simulaciones contiene 2.944 nodos y 5.686 elementos, b) Representación de las ecuaciones (15), c) representación de las ecuaciones (16).

Deformación senoidal con una oscilación dado por (ver fig. 2 b)

$h(x,y,t)=Acos(x_{t}\pi )$

( 15)

con un máximo localizado (ver fig. 2 c):

$h(x,y,t)=A\left(cos(2x_{t}\pi )+cos(2y_{t}\pi )\right)$

( 16)

Cada una de estas ecuaciones se puede ver representada en la figura 2 .

3. Ejemplos de aplicación y resultados numéricos

Todas las simulaciones se han hecho con un paso de tiempo Δt  = 0, 05. El tiempo final de la simulación es t  = 80, 0 adimensional. La malla que se utilizo se muestra en la fig. 2 , y contiene 2.944 nodos y 5.686 elementos.

3.1. Primer conjunto de parámetros

Se han hecho 2 tipos de simulaciones, según las ecuaciones (15) y (16). En el primer conjunto de parámetros que se ha utilizado esta dado por a  = 0, 1, b  = 0, 9, γ  = 29, 0 y d  = 10, 0, para el sistema de reacción-difusión y A=0,02 para las ecuaciones (15) y (16).

La deformación se lleva a cabo mediante el movimiento de los nodos según un sistema de ecuaciones como los descritos en (15) y (16) (fig. 2 ). Para la fila A se utiliza la ecuación (15), es decir se utiliza una función senoidal con una oscilación. En la fila B se utiliza la ecuación (16), es decir con un valor máximo en la función de deformación.

Con la tabla 1 , se puede establecer que estos parámetros dan un modo de onda (1,0) [28] . En la figura 3 se observa la formación del patrón de onda (1,0). Se puede notar la formación de un patrón con alta concentración en la parte inferior y baja concentración en la parte superior. Nótese que existe una alta deformación en la parte central del dominio (fig. 2 ).

Figura 3. Bloques 1 y 2 representan los parámetros de las ecuaciones (1) y (3) así: 1: a  = 0, 1, b  = 0, 9, γ  = 29, 0 y d  = 10, 0. Fila A: representan la ecuación (15); fila B: ecuación (16).

En la fila 1B, de la figura 3 , se muestran los resultados según la deformación descrita por la ecuación (16). Nótese que, en t=60, la concentración disminuye en el centro del dominio, donde se presenta la mayor deformación (fig. 2 ). En el siguiente instante de tiempo (t=80) se forma un halo de baja concentración alrededor de aquellas zonas de mayor gradiente de deformación. Se debe advertir que no se llega al estado estable, por tanto se tiene un cambio permanente en los patrones debido a la deformación de la superficie.

3.2. Primer conjunto de parámetros

Se han hecho 2 tipos de simulaciones, según las ecuaciones (15) y (16). En el primer conjunto de parámetros que se ha utilizado esta dado por a  = 0, 1, b  = 0, 9, γ  = 230, 82 y d  = 8, 6676, para el sistema de reacción-difusión y A=0,02 para las ecuaciones (15) y (16).

En la fig. 4 se presentan los resultados para un conjunto de parámetros que forman patrones de Turing (2,2), según la tabla 1[28] . En la fila 2a) [con deformación según la ecuación (15)], en t=60, se observa la formación de un patrón con modo de onda (2,2). Conforme pasa el tiempo se establece un patrón con un mínimo en el centro del dominio. En la parte lateral derecha e izquierda se forman manchas con mayor cercanía, más estrechas. En la figura 2 B [que se deforma según la ecuación (16)] se forman manchas de alta concentración en el centro y un conjunto de anillos que oscilan entre baja y alta concentración.

Figura 4. Bloques 2, se utilizan los siguientes parámetros a  = 0, 1, b  = 0, 9, γ  = 230, 82 y d  = 8, 6676. Las filas representan diferentes ecuaciones de deformación de la superficie: Fila A: representan la ecuación (15); fila B: ecuación (16).

4. Discusión y conclusiones

En este artículo se utilizó el método de los elementos finitos para solucionar un sistema de ecuaciones de reacción difusión de carácter no lineal. Para estudiar la versatilidad del método se desarrollaron varios ejemplos de reacción-difusión, con un sistema reactivo del tipo Schnakenberg. El sistema no lineal se ha resuelto mediante el método de Newton Raphson. Los parámetros del modelo cumplen las restricciones del espacio de Turing, por tanto, los patrones que se obtienen en la solución exhiben inestabilidades espaciales, denominadas: de Turing.

En este trabajo se logró modelar la deformación de superficies en donde se desarrolla el proceso de reacción y difusión, similar a lo que ocurriría en, por ejemplo, la epidermis de una animal, donde un proceso reactivo permite la formación de manchas características de la piel de algunas especies. Para este fin se utilizó la formulación de mecánica de medio continuo que permite aproximar la solución mediante un esquema lagrangiano total. En este método se utilizó la configuración de referencia [o indeformada (1)] para llevar a cabo los cálculos numéricos en cada instante de tiempo. Esta formulación permite generalizar la metodología propuesta por Madzvamuse et al. [38]  and [39] , en donde se desarrolla un modelo de movimiento de malla para problemas bidimensionales planos. En este artículo se logra ampliar dicha formulación para cualquier movimiento espacial de una superficie en deformación.

De los resultados (Figura 3  and Figura 4 ) se puede concluir que debido a la deformación de la superficie se pueden obtener patrones que no corresponden a aquellos predichos por la teoría de espacios de Turing (ver referencia [28] y tabla 1 ). La distribución de concentración en cada ejemplo sigue el patrón del movimiento de la malla. En este sentido se pierde el efecto del sistema reacción difusión en el espacio de Turing y se obtienen patrones que tienen máximos (o mínimos) en aquellas zonas de mayor (o menor) desplazamiento de los puntos del dominio. Este resultado puede explicar porque en algunas zonas de la piel de aquellos animales que exhiben manchas en su dermis puede cambiar el patrón superficial, o, inclusive desaparecer.

Estos resultados pueden compararse con los fenómenos que ocurren en la epidermis de los animales que presentan patrones especiales en su piel. Por tanto se puede especular que en aquellas zonas del cuerpo donde se presenta una gran deformación de la piel se puede perder los patrones superficiales y las manchas características. Por el contrario, en zonas de baja deformación pueden cambiar los patrones presentando manchas localizadas diferentes al resto de la distribución de manchas en la piel.

Bibliografía

1. [1] I. Babuska, F. Ihlenburg, E. Paik, S. Sauter; A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution; Comp Methods Appl Mech Eng, 128 (1995), pp. 325–359
2. [2] D. Garzon-Alvarado, J. García-Aznar, M. Doblare; Appearance and location of secondary ossification centres may be explained by a reaction-diffusion mechanism; Comput Biol Med, 39 (6) (2009), pp. 554–561
3. [3] S. Ferreira, M. Martins, M. Vilela; Reaction-diffusion model for the growth of avascular tumor; Phys Rev, 65 (2) (2002)
4. [4] M. Chaplain, A. Ganesh, I. Graham; Spatio-temporal pattern formation on spherical surfaces: Numerical simulation and application to solid tumor growth; J Math Biol, 42 (2001), pp. 387–423
5. [5] A. Madzvamuse; A numerical approach to the study of spatial pattern formation in the ligaments of arcoid bivalves; Bull Math Biol, 64 (2002), pp. 501–530
6. [6] S. Kondo, R. Asai; A reaction-diffusion wave on the skin of the marine anglefish, Pomacanthus; Nature, 376 (1995), pp. 765–768
7. [7] F. Crauste, M. Lhassan, A. Kacha; A delay reaction-diffusion model of the dynamics of botulinum in fish; Math Biosci, 216 (2008), pp. 17–29
8. [8] F. Rossi, S. Ristori, M. Rustici, N. Marchettini, E. Tiezzi; Dynamics of pattern formation in biomimetic systems; J Theor Biol, 255 (2008), pp. 404–412
9. [9] B. Rothschild, J. Ault; Population-dynamic instability as a cause of patch structure; Ecol Model, 93 (1996), pp. 237–239
10. [10] T. Nozakura, S. Ikeuchi; Formation of dissipative structures in galaxies; Astrophys J, 279 (1984), pp. 40–52
11. [11] R. Smith; Optimal and near-optimal advection-diffusion finite-difference schemes iii. Black-Scholes equation; Proc Math Phys Eng Sci, 456 (2000), pp. 1019–1028
12. [12] O. Richter; Modelling dispersal of populations and genetic information by finite element methods; Environ Model Software, 23 (2008), pp. 206–214
13. [13] L. Ferragut, M. Asensio, S. Monedero; A numerical method for solving convection–reaction–diffusion multivalued equations in fire spread modelling; Adv Eng Software, 38 (2007), pp. 366–371
14. [14] D. Garzón. Simulación de procesos de reacción-difusión: Aplicación a la morfogénesis del tejido óseo, Ph.D. Thesis. Universidad de Zaragoza. 2007.
15. [15] A. Madzvamuse, A. Wathen, P. Maini; A moving grid finite element method applied to a model biological pattern generator; J Comput Phys, 190 (2003), pp. 478–500
16. [16] O. Hirayama, R. Takaki; Thermal convection of a fluid with temperature-dependent viscosity; Fluid Dyn Res, 12–1 (1988), pp. 35–47
17. [17] M. Ardes, F. Busse, J. Wicht; Thermal convection in rotating spherical shells; Phys Earth Planet Int, 99 (1997), pp. 55–67
18. [18] J. Lir, T. Lin; Visualization of roll patterns in Rayleigh–Bénard convection of air in rectangular shallow cavity; Int J Heat Mass Trans, 44 (2001), pp. 2889–2902
19. [19] Y. Balkarei, A. GrigorYants, Y. Rhzanov, M. Elinson; Regenerative oscillations, spatial-temporal single pulses and static inhomogeneous structures in optically bistable semiconductors; Opt. Commun, 66 (1988), pp. 161–166
20. [20] V.I. Krinsky; Self-organisation: Auto-waves and structures far from equilibrium; Ed Springer (1984)
21. [21] L. Zhang, S. Liu; Stability and pattern formation in a coupled arbitrary order of autocatalysis system; Appl Math Model, 33 (2009), pp. 884–896
22. [22] H. Frederik, P. Maini, A. Madzvamuse, A. Wathen, T. Sekimura; Pigmentation pattern formation in butterflies: experiments and models; C R Biol, 326 (2003), pp. 717–727
23. [23] F. Yi, J. Wei, J. Shi; Bifurcation and spatiotemporal patterns in a homogeneous diffusive predator–prey system; J of Diff Eq, 246 (2009), pp. 1944–1977
24. [24] M. Baurmanna, T. Gross, U. Feudel; Instabilities in spatially extended predator–prey systems: Spatio-temporal patterns in the neighborhood of Turing–Hopf bifurcations; J Theor Biol, 245 (2007), pp. 220–229
25. [25] T. Nozakura, S. Ikeuchi; Formation of dissipative structures in galaxies; Astrophys. J, 279 (1984), pp. 40–52
26. [26] J. García-Aznar, J. Kuiper, M. Gómez-Benito, M. Doblaré, J. Richardson; Computational simulation of fracture healing: Influence of interfragmentary movement on the callus growth; J Biomech, 40 (3) (2007), pp. 1467–1476
27. [27] Z. Mei; Numerical Bifurcation Analysis for reaction-diffusion equations; Springer Verlag, Alemania (2000)
28. [28] Madvamuse A. A numerical approach to the study of spatial pattern formation D. Phil. Thesis. Oxford University: UK. 2000.
29. [29] S. Rüdiger, E.M. Nicola, J. Casademunt, L. Kramer; Theory of pattern forming systems under travelling-wave forcing; Phys Rep, 447 (3–6) (2007), pp. 73–111
30. [30] F. Sagués, D.G. Míguez, E.M. Nicola, A.P. Muñuzuri, J. Casademunt, L. Kramer; Travelling-stripe forcing of Turing patterns; Phys D: Nonlin Phenom, 199 (1–2) (2004), pp. 235–242
31. [31] L. Eugene, K.G. Allgower; Numerical path following; Handbook of Num Anal, 5 (1997), pp. 3–207
32. [32] K.J. Painter, H.G. Othmer, P.K. Maini; Stripe formation in juvenile Pomacanthus via chemotactic response to a reaction-diffusion mechanism; Proc Nat Acad Sci USA, 96 (1999), pp. 5549–5554
33. [33] K.J. Painter, P.K. Maini, H.G. Othmer; A chemotactic model for the advance and retreat of the primitive streak in avian development; Bull Math Biol, 62 (2000), pp. 501–525
34. [34] A. Madzvamuse; Turing instability conditions for growing domains with divergence free mesh velocity; Nonlin Anal: Theory Meth Appl, 71 (12) (2009), pp. 2250–2257
35. [35] J. Lefevre, J. A Mangin; Reaction-diffusion model of human brain development; PLOS Comput Biol, 6 (4) (2010), p. 1-10
36. [36] S. Seirin Lee, E.A. Gaffney, R.E. Baker; The dynamics of Turing patterns for morphogen-regulated growing domains with cellular response delays; Bull Math Biol (2011) http://dx.doi.org/10.1007/s11538-011-9634-8
37. [37] A. Turing; The chemical basis of morphogenesis; Phil. Trans. R. Soc. Lond. B, 237 (1952), pp. 37–72
38. [38] A. Madzvamuse, T. Sekimura, R.D.K. Thomas, A.J. Wathen, P.K. Maini; A moving grid finite element method for the study of spatial pattern formation in biological problems; T. Sekimura, S. Noji, N. Nueno, P.K. Maini (Eds.), Morphogenesis and Pattern Formation in Biological Systems – Experiments and Models, Springer-Verlag, Tokyo (2003), p. 59-65
39. [39] A. Madzvamuse, P.K. Maini, A.J. Wathen; A moving grid finite element method applied to a model biological pattern generator; J Comp Phys, 190 (2) (2003), pp. 478–500
40. [40] A. Geirer, H. Meinhardt; A theory of biological pattern formation; Kybernetik, 12 (1972), pp. 30–39
41. [41] A. De Wit; Spatial patterns and spatiotemporal dynamics in chemical systems; Adv Chem Phys, 109 (1999), pp. 435–513
42. [42] P.K. Maini, K.J. Painter, H.N.P. Chau; Spatial pattern formation in chemical and biological systems; Faraday Trans, 93 (20) (1997), pp. 3601–3610
43. [43] R. Kapral, K. Showalter; Chemical Waves and Patterns; Kluwer (1995)
44. [44] J. Lefevre, J. Mangin; A reaction-diffusion model of human brain development; PLOS Comput Biol, 6 (2010), pp. 1–10
45. [45] L. Harrison, S. Wehner, D. Holloway; Complex morphogenesis of surfaces: theory and experiment on coupling of reaction-diffusion patterning to growth; Faraday Dis, 120 (2002), pp. 277–293
46. [46] D. Holloway, L. Harrison; Pattern selection in plants: Coupling chemical dynamics to surface growth in three dimensions; Ann Bot, 101 (2008), p. 361
47. [47] T.J.R. Hughes; The finite element method: linear static and dynamic finite element analysis; Courier Dover Publications, New York (2003)
48. [48] J. Hoffman; Numerical methods for engineers and scientists; Ed McGraw Hill (1992)
49. [49] G.A. Holzapfel; Nonlinear solid mechanics; John Wiley (2000)
50. [50] T. Belytschko, W.K. Liu, B. Moran; Nonlinear finite elements for continua and structures; John Wiley and Sons (2000) Disponible en: http://www.amazon.com/Nonlinear-Finite-Elements-Continua-Structures/dp/0471987743