14.4: Análisis de Estabilidad Lineal de Sistemas de Reacción-Difusión

Es posible que haya encontrado que el análisis de estabilidad lineal de los modelos de campo continuo no es tan fácil como el de los modelos no espaciales. Para esto último, tenemos una herramienta muy conveniente llamada las matrices jacobinas, y el análisis de estabilidad es simplemente calcular una matriz jacobiana y luego investigar sus valores propios. Todo es tan mecanicista y automático, comparado con lo que pasamos en el apartado anterior. Quizás te preguntes, ¿no hay atajos más fáciles para analizar la estabilidad de los modelos de campo continuos?

Bueno, si te sientes así, te convertirás en un gran fan de los sistemas de reacción-difusión que discutimos en la Sección 13.6. Su análisis de estabilidad lineal es mucho más fácil, debido a la clara separación de la dinámica de reacción local y la dinámica de difusión espacial. Para ser más específicos, ¡puedes devolver la matriz jacobiana al análisis! Aquí está cómo y por qué funciona.

Considere la posibilidad de realizar un análisis de estabilidad lineal con el siguiente sistema estándar de reacción-difusión:

\[\frac{\partial{f_{1}}}{\partial{t}} =R_{1} (f_{1},f_{2}, \cdots, f_{n})+D_{1}\nabla^{2} f_{1} \label{(14.80)} \]

\[\frac{\partial{f_{2}}}{\partial{t}} =R_{2}(f_{1}, f_{2}, \cdots, f_{n})+D_{2}\nabla^{2}f_{2} \label{(14.81)} \]

\[\vdots \nonumber \]

\[\frac{\partial{f_{n}}}{\partial{t}} =R_{n}(f_{1}, f_{2}, \cdots, f_{n})+D_{n}\nabla^{2}f_{n} \label{(14.82)} \]

El estado de equilibrio homogéneo de este sistema\((f_{1eq},f_{2eq},...,f_{neq})\),, es una solución de las siguientes ecuaciones:

\[0=R_{1}(f_{1eq}, f_{2eq}, \cdots. f_{neq}) \label{(14.83)} \]

\[0=R_{2}(f_{1eq}, f_{2eq}, \cdots. f_{neq}) \label{(14.84)} \]

\[\vdots \nonumber \]

\[0=R_{n}(f_{1eq}, f_{2eq}, \cdots. f_{neq}) \label{(14.85)} \]

Para realizar un análisis de estabilidad lineal, reemplazamos las variables de estado originales de la siguiente manera:

\[f_{i}(x,t) \Rightarrow f_{ieq} +\Delta{f_{i}(x,t)} =f_{ieq} +\sin{(\omega{x} +\phi)}\Delta{f_{1}(t)} \text{for all i}\label{(14.86)} \]

Este reemplazo convierte las ecuaciones dinámicas en la siguiente forma:

\[ S\dfrac{\partial{\Delta}f_{1}}{\partial{t}} =R_{1}(f_{1eq} +S\Delta{f_{1}}. f_{2eq} +S\Delta{f_{2}}, \cdots f_{neq} +S\Delta{f_{n}}) -D_{1}\omega^{2}S\Delta{f_{1}} \label{(14.87)} \]

\[ S\dfrac{\partial{\Delta}f_{2}}{\partial{t}} =R_{2}(f_{1eq} +S\Delta{f_{1}}. f_{2eq} +S\Delta{f_{2}}, \cdots f_{neq} +S\Delta{f_{n}}) -D_{2}\omega^{2}S\Delta{f_{2}} \label{(14.88)} \]

\[\vdots \nonumber \]

\[ S\frac{\partial{\Delta}f_{n}}{\partial{t}} =R_{n}(f_{1eq} +S\Delta{f_{1}}, f_{2eq} +S\Delta{f_{2}}, \cdots f_{neq} +S\Delta{f_{n}}) -D_{n}\omega^{2}S\Delta{f_{n}} \label{(14.89)} \]

Aquí utilicé\(S = sin(ωx + φ)\) sólo en las expresiones anteriores para acortarlas. Estas ecuaciones se pueden resumir en forma de un solo vector sobre\(∆f\),

\[\sin{(\omega{x} +\phi)}\frac{\partial{\Delta{f}}}{\partial{t}} =R(f_{eq} +\sin{(\omega{x} +\phi)}\Delta{f}) -D\omega^{2}\sin{\omega{x} +\phi)}\Delta{f}, \label{(14.90)} \]
donde\(R\) es una función vectorial que representa todos los términos de reacción, y\(D\) es una matriz diagonal cuyos componentes diagonales son\(D_i\) para la posición i-ésima. Ahora que se han simplificado todos los términos de difusión, si también podemos linealizar los términos de reacción, podemos completar la tarea de linealización. Y aquí es donde la matriz jacobiana vuelve a ser puesta en el centro de atención. Los términos de reacción son todos locales sin ningún operador espacial involucrado, y por lo tanto, a partir de la discusión en la Sección 5.7, sabemos que la función vectorial\(R(f_{eq} + \sin{(ωx + φ)}∆f)\) puede aproximarse linealmente de la siguiente manera:

\[R(f_{eq} +\sin{(\omega{x} +\phi)\Delta{f}})\approx R(f_{eq}) + \begin{pmatrix} \dfrac{\partial{R_{1}}}{\partial{f_{1}}} & \dfrac{\partial{R_{1}}}{\partial{f_{2}}} & \cdots & \dfrac{\partial{R_{1}}}{\partial{f_{n}}} \\ \dfrac{\partial{R_{2}}}{\partial{f_{1}}} & \dfrac{\partial{R_{2}}}{\partial{f_{2}}} & \cdots & \dfrac{\partial{R_{2}}}{\partial{f_{n}}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial{R_{n}}}{\partial{f_{1}}} & \dfrac{\partial{R_{n}}}{\partial{f_{2}}} & \cdots & \dfrac{\partial{R_{n}}}{\partial{f_{n}}} \end{pmatrix}| _{f=f_{eq}} \sin{(\omega{x} +\phi)}\Delta{f} \label{(14.91)} \]

\[=\sin{(\omega{x} +\phi)} \begin{pmatrix} \dfrac{\partial{R_{1}}}{\partial{f_{1}}} &\dfrac{\partial{R_{1}}}{\partial{f_{2}}} &\cdots & \dfrac{\partial{R_{1}}}{\partial{f_{n}}} \\ \dfrac{\partial{R_{2}}}{\partial{f_{1}}} &\dfrac{\partial{R_{2}}}{\partial{f_{2}}} &\cdots & \dfrac{\partial{R_{2}}}{\partial{f_{n}}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial{R_{n}}}{\partial{f_{1}}} &\dfrac{\partial{R_{n}}}{\partial{f_{2}}} &\cdots & \dfrac{\partial{R_{n}}}{\partial{f_{n}}} \end{pmatrix} |_{f=f_{eq}} \Delta{f}\label{(14.92)} \]

Tenga en cuenta que podemos eliminar\(R(f_eq)\) debido a las Eqs. \ ref {(14.83)},\ ref {(14.84)}\ ref {(14.85)}. Al conectar este resultado a la Ecuación\ ref {(14.90)}, obtenemos
\[\sin{(\omega{x} +\phi)} \frac{\partial{\Delta{f}}}{\partial{t}} =\sin{(\omega{x} +\phi)}J|_{f=f_{eq}} \Delta{f} =D{\omega^{2}} \sin{(\omega{x} +\phi)}\Delta{f}, \label{(14.93)} \]

\[\frac{\partial{\Delta{f}}}{\partial{t}} - (J-D{\omega^{2}})|_{f=f_{eq}}\Delta{f}, \label{(14.94)} \]

donde\(J\) está la matriz jacobiana de los términos de reacción (\(R\)). ¡Muy sencillo! Ahora solo necesitamos calcular los valores propios de esta matriz de coeficientes para estudiar la estabilidad del sistema.

La estabilidad de un sistema de reacción-difusión en su estado de equilibrio homogéneo se\(f_{eq}\) puede estudiar calculando los valores propios de

\[(J-D\omega^{2})|_{f=f_{eq}}, \label{(14.95)} \]

donde\(J\) está la matriz jacobiana de los términos de reacción,\(D\) es la matriz diagonal hecha de constantes de difusión, y\(w\) es un parámetro que determina la frecuencia espacial de las perturbaciones.

Este atajo en el análisis de estabilidad lineal es posible gracias a la clara separación de términos de reacción y difusión en los sistemas de reacción-difusión. Espero que entiendas ahora parte de las razones por las que tantos investigadores son aficionados a este marco de modelado.

Apliquemos este nuevo conocimiento a algún ejemplo. Aquí está el modelo de Turing que discutimos antes:

\[\frac{\partial{u}}{\partial{t}} = a(u-h) +b(v-k)+D_{u}\nabla^{2}{u} \label{(14.96)} \]

\[\frac{\partial{v}}{\partial{t}} =c(v-h) +d(v-k) +D_{v} \nabla^{2{v} \label{(14.97)}} \]

Usando la ecuación\ ref {(14.95)}, podemos calcular inmediatamente su matriz de coeficientes:

\[(\begin{pmatrix} a & \\ c &d \end{pmatrix} - \begin{pmatrix} D_{u} & 0 \\ 0 & D_{v} \end{pmatrix} \omega^{2})|_{(u, v)=(h, k)} = \begin{pmatrix} a-D_{u}\omega^{2} & b \\ c & d-D_{v} \omega^{2} \end{pmatrix} \label{(14.98)} \]

De la discusión en la Sección 7.4, ya sabemos que, para que esta matriz muestre estabilidad, su determinante debe ser positivo y su traza debe ser negativa (ver Fig. 7.4.1). Por lo tanto, la condición para que el estado de equilibrio homogéneo de este sistema sea estable es que ambas de las dos desigualdades siguientes deben ser ciertas para todos los valores reales de\(ω\):

\[0 < (a- D_{u} \omega^{2})(d-D_{v}\omega^{2}) -bc \label{(14.99)} \]

\[0 > a- D_{u}\omega^{2} +d-D_{v} \omega^{2} \label{(14.100)} \]

Estas desigualdades se pueden reescribir usando\(det(A)\) y\(Tr(A)\) de\( A =\begin{pmatrix}a & b \\ c & d \end{pmatrix}\), de la siguiente manera:

\[aD_{v}\omega^{2} +dD_{u}\omega^{2}-D_{u}D_{v}\omega^{4} <\det{(A)} \label{(14.101)} \]

\[D_{u}\omega^{2} +D_{v}\omega^{2} > Tr{(A)} \label{(14.102)} \]

Ahora, imagina que el modelo no espacial original sin términos de difusión ya era estable, es decir,\(det(A) > 0 \) y\(Tr(A) < 0\). ¿Existe alguna posibilidad de que la introducción de la difusión al modelo pueda desestabilizar el sistema por sí mismo? La segunda desigualdad siempre es cierta para negativo\ (Tr (A)), porque su lado izquierdo no puede ser negativo. Pero la primera desigualdad puede ser violada, si

\[g(z) =-D_{u}D_{v}z^{2} +(aD_{v} +dD_{u})z -\det{(A)} \qquad {(with \ z = \omega^{2})} \label{(14.103)} \]

puede tomar un valor positivo para algunos\(z > 0\). \(g(z)\)se puede reescribir como

\[g(z) =-D_{u}D_{v} (z -\dfrac{aD_{v} +dD_{u}}{2D_{uD_{v}}})^{2} +\frac{(aD_{v} +dD_{u})^{2}}{4D_{u}D_{v}} -\det{(A)}. \label{(14.104)} \]

Existen dos escenarios potenciales en los que este polinomio puede ser positivo para algunos\(z > 0\), como se muestra en la Fig. 14.4.1.

Figura\(\PageIndex{1}\): Dos escenarios posibles en los que\(g(z)\) puede tomar un valor positivo para algunos\(z > 0\). Izquierda: Cuando el pico existe en el lado positivo de\(z\). Derecha: Cuando el pico existe en el lado negativo de\(z\).

Si el pico existe en el lado positivo de\(z (aD_{v} + dD_{u} > 0\); Fig. 14.4.1 izquierda), la única condición es que el pico sobresalga por encima del\(z\) eje -eje, i.e.

\[\dfrac{(aD_{v} +dD_{u})^{2}}{4D_{u}D_{v}} -\det{(A)} >0. \label{(14.105)} \]

O bien, si el pico existe en el lado negativo\(z (aD_{v} +dD_{u} < 0\); Fig. 14.4.1 derecha), la condición es que la intercepción de\(g(z)\) sea positiva, es decir,

\[g(0) =-\det{(A)} >0, \label{(14.106)} \]

pero esto no puede ser cierto si el modelo no espacial original es estable. Por lo tanto, la única posibilidad de difusión para desestabilizar el sistema por lo demás estable es el primer caso, cuya condición puede simplificarse para

\[aD_{v} +dD_{u} >2 \sqrt{D_{u}D_{v} \det{(A).}} \label{(14.107)} \]

Probemos qué tan bien se aplica esta predicción a la dinámica real de los modelos de Turing. En el capítulo anterior, se utilizó\((a,b,c,d) = (1,−1,2,−1.5)\) y\((D_{u},D_{v}) = (10−^{4},6 × 10^{−4)}\) para generar el resultado de simulación mostrado en la Fig. 13.17. Con estos ajustes de parámetros, det (A) = −1.5− (−2) = 0.5 > 0 y Tr (A) = −0.5 < 0, por lo que el sistema sería estable si no hubiera términos de difusión. Sin embargo,

\[aD_{v} +dD_{u} = 6×10^{−4 }−1.5×106{−4} = 4.5×10^{−4}, \label{(14.108)} \]

\[2\sqrt{D_{u}D_{v}\det{(A)}} =2\sqrt{10^{−4} ×6×10^{−4} ×0.5 } = 2×10^{−4} \sqrt{3} \approx 3.464 ×10^{−4} \label{(14.109)} \]

por lo tanto la desigualdad\ ref {(14.107)} se mantiene. Esto indica que el estado de equilibrio homogéneo debe ser inestable y deben surgir patrones espaciales no homogéneos que realmente se pueden ver en la Fig. 13.17. Como Briefly se menciona en la Sección 13.6, esto se denomina inestabilidad inducida por difusión. Se trata de un fenómeno bastante contrario a la intuición, porque la difusión suele considerarse un proceso donde una estructura no homogénea está siendo destruida por el movimiento aleatorio. Pero aquí, el sistema es estable en su estado homogéneo con difusión, pero puede crear espontáneamente estructuras no homogéneas con difusión. Este es un muy buen ejemplo de lo alucinante que puede ser a veces el comportamiento de sistemas complejos.

Ejercicio\(\PageIndex{1}\)

Realizar un análisis de estabilidad lineal del modelo de depredador-presa espacialmente extendido, alrededor de su estado de equilibrio homogéneo distinto de cero, y discutir los resultados:

\[\frac{\partial{r}}{\partial{t}} =ar-brf +D_{r}\nabla^{2}r \label{(14.110)} \]

\[\frac{\partial{f}}{\partial{t}} =-cf -drf +D_{f} \nabla^{2}f \label{(14.111)} \]

Supongamos que todos los parámetros del modelo son positivos.

Ejercicio\(\PageIndex{2}\)

Realizar un análisis de estabilidad lineal del modelo de Gray-Scott, alrededor de su estado de equilibrio homogéneo\((u_{eq},v_{eq}) = (1,0)\), y discutir los resultados:

\[\frac{\partial{u}}{\partial{t}} =F(1-u) uv^{2} +D_{u} \nabla^{2}u \label{(14.112)} \]

\ [\ frac {\ parcial {v}} {\ parcial {t}} =- (F+k) v +uv^ {2} +D_ {v}\ nabla^ {2} v\ etiqueta {(14.113)}]

Nuevamente, supongamos que todos los parámetros del modelo son positivos.

Hay algunas predicciones más útiles que podemos hacer sobre la formación espontánea de patrones en sistemas de reacción-difusión. Sigamos usando como ejemplo el modelo de Turing discutido anteriormente. Podemos calcular los valores propios reales de la matriz de coeficientes, de la siguiente manera:

\[\begin{align} \begin{vmatrix} 1−10−4ω^{2} −λ & -1 \\ 2 & −1.5−6×10−4ω^{2} −λ \end {vmatrix}|=0 \label{(14.114)} \\[4pt] (1 −10−4ω^{2} −λ)(−1.5−6×10−4ω^{2} −λ) -(-2) =0 \label{(14.115)} \\[4pt] \lambda^{2} + (0.5 + 7×10−4ω^{2})λ + (1−10^{−4}ω^{2})(−1.5−6×10^{−4}ω^{2}) + 2 = 0 \label{(14.116)} \\[4pt] \lambda =\frac{1}{2} (-(0.5 + 7×10^{−4}ω^{2} ) \pm \sqrt{(0.5 + 7×10^{−4}ω^{2})2 −4(1−10^{−4}ω^{2})(−1.5−6×10^{−4}ω^{2})−8} ) \label{(14.117)} \\[4pt] =\frac{1}{2} −(0.5 + 7×10^{−4}ω^{2})±\sqrt{2.5×10^{−7}w^{4} + 2.5×10^{−3}w^{2} −1.75}) \label{(14.118)} \end{align} \]

De estos dos valores propios, el que podría tener una parte real positiva es el que tiene el signo “\(+\)” (llamémoslo\ (λ_ {+}).

Aquí, lo que vamos a hacer es calcular el valor de\(ω\) que alcance la mayor parte real de\(λ_+\). Esta es una pregunta significativa, porque la mayor parte real de los valores propios corresponde a la función propia dominante (\(\sin{(ωx + φ)}\)) que crece más rápido, que debería ser el patrón espacial más visible que surge en el estado del sistema. Si descubrimos tal valor de\(ω\), entonces nos\(2π/ω\) da la escala de longitud de la función propia dominante.

Podemos obtener una respuesta a esta pregunta analizando dónde\(λ_+\) ocurre el extremo de. Para simplificar el análisis, dejamos de\(z = ω^2\) nuevo y usamos\(z\) como variable independiente, de la siguiente manera:

\[\begin{align} \frac{d\lambda_{+}}{dz} =\frac{1}{2}( -7×10^{-4}+\frac{5 ×10^{-7}z +2.5 × 10^{-3}}{2\sqrt{2.5×10^{-7} +2.5×10^{-3}z -1.75}}) =0 \label{(14.119)} \\[4pt] 7× 10^{-4}(2 \sqrt{2.5×10^{-7}z^{2} 2.5×10^{-3}z-1.75)} =5×10^{-7}z +2.5×10^{-3} \label{(14.120)} \\[4pt] 1.96 ×10^{-6}(2.5×10^{-7}z^{2} +2.5 ×10^{-3}z-1.75) =2.5 ×10^{-13}z^{2}+2.5 × 10^{-9}z +6.25 × 10^{-6} \label{(14.121)} \\[4pt] (... blah \ blah \ blah ...) \\[4pt] 2.4 × 10^{-13} z^{2} +2.4× 10^{-9}z -9.68× 10^{-6} =0 \label{(14.122)} \\[4pt] z = 3082.9, −13082.9 \label{(14.123)}\end{align} \]

Uf. Agotado. De todos modos\(z = ω^2\), ya que, el valor de\(ω\) eso corresponde a la función propia dominante es

\[\omega = \sqrt{3082.9} =55.5239. \label{(14.124)} \]

Ahora podemos calcular la escala de longitud de la función propia dominante correspondiente, que es
\[\ =\dfrac{2\pi}{\omega} \approx 0.113162. \label{(14.125)} \]

Este número da la distancia característica entre las crestas (o entre los valles) en la función propia dominante, la cual se mide en longitud unitaria. 0.11 significa que hay alrededor de nueve crestas y valles en una unidad de longitud. De hecho, la simulación mostrada en la Fig. 13.17 se realizó en un dominio cuadrado\([0,1]×[0,1]\) unitario, así que continúe y cuente cuántas ondas se alinean a lo largo de su\(x\) eje o\(y\) eje en su configuración final.

¿Has terminado de contarlos? Sí, ¡el resultado de la simulación de hecho mostró alrededor de nueve ondas en cada eje! Este es un ejemplo de cómo podemos predecir no sólo la estabilidad del estado de equilibrio homogéneo, sino también la escala de longitud característica de los patrones que se forman espontáneamente si el estado de equilibrio resulta inestable.

Ejercicio\(\PageIndex{3}\)

Estimar la escala de longitud de los patrones que se forman en el modelo de Turing anterior con\((a,b,c,d) = (0.5,−1,0.75,−1)\) y\((D_u,D_v) = (10^{−4},10^{−3})\). Después confirme su predicción con simulaciones numéricas.

Ejercicio\(\PageIndex{4}\)

Si ambas constantes de difusión se multiplican por el mismo factor\( ψ\), ¿cómo afecta eso a la escala de longitud de los patrones?

Bien, hagamos solo una predicción más, y ya terminaremos. Aquí predecimos la relación crítica de las dos constantes de difusión, en la que el sistema se encuentra justo en el umbral entre la homogeneización y la formación del patrón. Usando un nuevo parámetro\(ρ = D_{v}/D_{u}\), la condición de inestabilidad (desigualdad\ ref {(14.107)}) se puede simplificar aún más de la siguiente manera:

\[a\rho{D_{u}} +dD_{u} > 2\sqrt{\rho{D_{u}^{2}} \det{(A)}} \label{(14.126)} \]

\[a \rho +d > 2\sqrt{\rho\det{(A)}} \label{(14.127)} \]

Para los valores de los parámetros que usamos anteriormente, esta desigualdad se resuelve de la siguiente manera:

\[\rho -1.5 > 2 \sqrt{0.5\rho}\label{(14.128)} \]

\[\rho^{2} -5\rho +2.25 >0 \text{(with ρ−1.5 > 0)} \label{(14.129)} \]

\[\rho > 4.5 \label{(14.130)} \]

Esto significa que la difusión de\(v\) debe ser al menos 4.5 veces más rápida que\(u\) para provocar la inestabilidad de difusión. En otras palabras, u actúa más localmente, mientras que los efectos del\(v\) alcance sobre rangos espaciales más largos. Si miras hacia atrás en la matriz de coeficientes original\(\begin{pmatrix} 1 & -1 \\ 2 & -1.5 \end{pmatrix}\), te darás cuenta de que\(u\) tiende a aumentar tanto\(u\) y\(v\), mientras\(v\) tiende a suprimir\(u\) y\(v\). Por lo tanto, esto representa dinámicas típicas de “activación de corto alcance e inhibición de largo alcance” que discutimos en la Sección 11.5, que es esencial en muchos procesos de formación de patrones.

La Figura 14.4.2 muestra los resultados de simulación numérica con la relación de las constantes de difusión sistemáticamente variada. De hecho, ¡en realidad se observa una transición brusca de los resultados a través\(ρ = 4.5\)! Este tipo de transición del comportamiento de un sistema de reacción-difusión entre la homogeneización y la formación de patrones se llama bifurcación de Turing, que el propio Turing mostró en su monumental trabajo en la década de 1950 [44].

Ejercicio\(\PageIndex{5}\)

A continuación se muestra una variante del modelo de formación del patrón de Turing:

\[\frac{\partial{u}}{\partial{t}} =u(v-1) -\alpha +D_{u} \nabla^{2}u \label{(14.131)} \]

\[\frac{\partial{v}}{\partial{t}} =\beta -uv +D_{v}\nabla^{2}v\label{(14.132)} \]

Aquí\(α\) y\(β\) son parámetros positivos. Vamos a\((α,β) = (12,16)\) lo largo de este ejercicio. Haga lo siguiente:

1. Encuentra su estado de equilibrio homogéneo.
2. Examinar la estabilidad del estado de equilibrio homogéneo sin términos de difusión.
3. Con\((D_{u},D_{v}) = (10^{−4},10^{−3})\), realizar un análisis de estabilidad lineal de este modelo alrededor del estado de equilibrio homogéneo para determinar si los patrones no homogéneos se forman espontáneamente. Si lo hacen, estime la escala de longitud de los patrones.
4. Determinar la relación crítica de las dos constantes de difusión.
5. Confirma tus predicciones con simulaciones numéricas.

Figura\(\PageIndex{2}\): Resultados de simulación numérica del modelo de formación del patrón de Turing con\((a,b,c,d) = (1,−1,2,−1.5)\),\(D_{u} = 10^{−4},\) y\(D_{v} = ρD_-{u}\). Las densidades de\(u\) se trazan en escala de grises (más oscuro = mayor). El valor de\(ρ\) se da debajo de cada resultado.