Año 2 · Módulo 15: Dinámica espacial y sistemas continuos

144.1. Cuando la uniformidad deja de ser la solución estable

Una sustancia que difunde tiende a suavizar sus diferencias de concentración. Si colocamos una pequeña cantidad en una región de un recipiente y dejamos que se distribuya, esperamos ver una mancha cada vez más ancha y con un máximo cada vez menor. Hay sistemas, sin embargo, en que el transporte entre regiones acaba favoreciendo la aparición de diferencias espaciales persistentes. Esta posibilidad motivó el artículo de 1952 de Alan Turing sobre la base química de la morfogénesis, una de sus contribuciones más influyentes a la biología matemática.

Turing imaginó sustancias que reaccionan entre sí y se desplazan por difusión dentro de un tejido. Un estado uniforme puede ser estable si cada posición evoluciona de manera aislada y, pese a ello, volverse inestable cuando se permite intercambiar sustancias entre posiciones vecinas. Una perturbación espacial diminuta puede amplificarse hasta producir regiones con concentraciones sistemáticamente diferentes, cuya geometría dependerá de las ecuaciones y del dominio en que actúen. La idea es una explicación posible de la emergencia de estructura sin necesidad de imponer desde el comienzo un dibujo externo.

No toda irregularidad que persiste es un patrón de Turing. Una fuente fija en una pared puede mantener un gradiente; una condición inicial muy desigual puede tardar largo tiempo en relajarse; un flujo externo puede arrastrar material y acumularlo. En esos ejemplos, las diferencias observadas no requieren una inestabilidad del equilibrio uniforme. La expresión inestabilidad inducida por difusión se reserva, en su formulación clásica, para la situación en que un equilibrio espacialmente homogéneo es estable frente a pequeñas perturbaciones uniformes y pierde estabilidad frente a perturbaciones espaciales de determinadas longitudes de onda al incorporar difusión.

Hay también que separar dos cuestiones. La primera es la aparición inicial de una perturbación creciente, que puede describirse mediante la linealización y una relación de dispersión. La segunda es el desarrollo de un patrón de amplitud finita, ya sean rayas, manchas, laberintos u otras distribuciones que pueden permanecer, desplazarse, oscilar o desintegrarse. Una tasa de crecimiento positiva anticipa que la uniformidad se rompe, mientras que las propiedades del patrón maduro dependen de la dinámica no lineal, la geometría, las fronteras y las perturbaciones iniciales.

Estudiaremos primero el mecanismo mínimo con dos campos químicos, derivaremos las condiciones de inestabilidad y resolveremos un modelo concreto. Después examinaremos cómo compiten distintos modos, qué papel cumple el tamaño del medio y cómo validar una simulación. Finalmente consideraremos la interpretación biológica de estas estructuras y los límites de una explicación que usa únicamente reacción y difusión. El objetivo es disponer de un criterio que permita reconocer cuándo una estructura espacial puede atribuirse con fundamento a ese mecanismo.

144.2. Dos especies que reaccionan localmente y se transportan

Denotemos por u(r,t)u\left( \mathbf{r},t \right) y v(r,t)v\left( \mathbf{r},t \right) las concentraciones, o las variables adimensionales apropiadas, de dos especies. En cada posición se transforman según las funciones f(u,v)f(u,v) y g(u,v)g(u,v). Si ambas difunden de modo aproximadamente isotrópico, con coeficientes constantes y positivos, el sistema tiene la forma

∂u∂t=f(u,v)+DuΔu,∂v∂t=g(u,v)+DvΔv.\frac{\partial u}{\partial t} = f(u,v) + D_{u}\Delta u,\quad\quad\frac{\partial v}{\partial t} = g(u,v) + D_{v}\Delta v.

Las funciones de reacción describen la cinética local, mientras que los laplacianos incorporan el transporte espacial. Si las variables son concentraciones físicas, cada término de reacción tiene unidades de concentración por tiempo y cada difusividad tiene unidades de longitud al cuadrado por tiempo. Si el sistema se ha adimensionalizado, los números de la ecuación ya incorporan las escalas de referencia; para volver a compararlos con un experimento habrá que recuperar esas escalas.

Un equilibrio homogéneo (u∗,v∗)\left( u_{*},v_{*} \right) satisface f(u∗,v∗)=g(u∗,v∗)=0f\left( u_{*},v_{*} \right) = g\left( u_{*},v_{*} \right) = 0. Como los campos son uniformes, sus laplacianos son cero. El equilibrio existe tanto en la dinámica espacial como en el sistema de ecuaciones diferenciales ordinarias u˙=f(u,v)\dot{u} = f(u,v), v˙=g(u,v)\dot{v} = g(u,v). Para conocer la estabilidad local sin transporte, calculamos el jacobiano

J=(fufvgugv)(u∗,v∗).\mathbf{J} = \begin{pmatrix} f_{u} & f_{v} \\ g_{u} & g_{v} \end{pmatrix}_{\left( u_{*},v_{*} \right)}.

Si su traza es negativa y su determinante positivo, los dos autovalores tienen parte real negativa y el equilibrio es linealmente estable frente a perturbaciones homogéneas. Puede haber relajación oscilatoria cuando los autovalores forman un par complejo, pero la amplitud tenderá a desaparecer. Esta estabilidad es el requisito inicial para que una pérdida posterior pueda atribuirse al transporte espacial, porque descarta que la reacción fuera ya inestable sin él.

Para un dominio aislado, una condición de frontera frecuente es n⋅∇u=n⋅∇v=0\mathbf{n} \cdot \nabla u = \mathbf{n} \cdot \nabla v = 0. Significa que no hay flujo difusivo a través de la pared cuando las difusividades son escalares constantes. Esta condición tiene dos consecuencias distintas. El transporte no introduce ni extrae material del dominio, aunque la cantidad de cada especie no tiene por qué conservarse, ya que las reacciones pueden producirla o consumirla. La condición determina además qué modos espaciales satisfacen las fronteras, de lo que depende cuáles podrán crecer.

Las interacciones entre las dos especies admiten diferentes interpretaciones. Una puede favorecer su propia producción y estimular o inhibir a la otra; también puede competir por un sustrato o transformar una especie en otra. La etiqueta «activador–inhibidor» ayuda en ciertas familias de modelos, pero no es una clasificación universal que se deduzca de cualquier matriz de dos variables. Para el análisis lineal, lo decisivo son las derivadas parciales de las funciones de reacción evaluadas en el equilibrio. Esas derivadas cuantifican cómo responde la dinámica local ante pequeñas variaciones y, combinadas con las difusividades, permiten decidir si el transporte puede invertir la estabilidad de alguna perturbación.

144.3. Derivación de las condiciones clásicas de Turing

Escribamos u=u∗+ηu = u_{*} + \eta y v=v∗+ζv = v_{*} + \zeta, donde η\eta y ζ\zeta son perturbaciones pequeñas. Al conservar los términos de primer orden de las reacciones obtenemos, para el vector w=(η,ζ)T\mathbf{w} = (\eta,\zeta)^{T}, el sistema lineal

∂tw=Jw+(Du00Dv)Δw.\partial_{t}\mathbf{w} = \mathbf{J}\mathbf{w} + \begin{pmatrix} D_{u} & 0 \\ 0 & D_{v} \end{pmatrix}\Delta\mathbf{w}.

En un medio homogéneo, busquemos una perturbación proporcional a eσt+ik⋅re^{\sigma t + i\mathbf{k} \cdot \mathbf{r}}. El laplaciano de esa expresión es −k2- k^{2} veces la propia función, donde k=∣k∣k = \left| \mathbf{k} \right|. Así, los autovalores temporales son los de la matriz

A(k2)=J−k2diag(Du,Dv).\mathbf{A}\left( k^{2} \right) = \mathbf{J} - k^{2}diag\left( D_{u},D_{v} \right).

Para abreviar, definamos q=k2q = k^{2}. Si el equilibrio local es estable, trJ<0tr\mathbf{J} < 0 y det⁡J>0\det\mathbf{J} > 0. Como Du,Dv>0D_{u},D_{v} > 0, la traza de A(q)\mathbf{A}(q) permanece negativa para todo q≥0q \geq 0. Con la traza negativa, solo puede haber crecimiento si su determinante llega a ser negativo. Al desarrollarlo aparece el polinomio

det⁡A(q)=DuDvq2−(Dvfu+Dugv)q+det⁡J.\det\mathbf{A}(q) = D_{u}D_{v}q^{2} - \left( D_{v}f_{u} + D_{u}g_{v} \right)q + \det\mathbf{J}.

Sea B=Dvfu+DugvB = D_{v}f_{u} + D_{u}g_{v}. Para que este polinomio, positivo en q=0q = 0 y creciente hacia +∞+ \infty cuando qq es grande, tome valores negativos en alguna región positiva, deben cumplirse

B>0,B2>4DuDvdet⁡J.\boxed{B > 0,\quad\quad B^{2} > 4D_{u}D_{v}\det\mathbf{J}.}

Junto a trJ<0tr\mathbf{J} < 0 y det⁡J>0\det\mathbf{J} > 0, estas son las condiciones clásicas para que exista una banda continua de números de onda linealmente inestables en este modelo de dos especies con difusión diagonal positiva. En un dominio acotado hay que comprobar además que algún autovalor espacial permitido cae dentro de la banda. Sin esta comprobación, el continuo infinito puede predecir crecimiento que la geometría del sistema real no admite.

La banda de cuadrados de número de onda se determina resolviendo det⁡A(q)=0\det\mathbf{A}(q) = 0. Sus límites son

q±=B±B2−4DuDvdet⁡J2DuDv,q−<q<q+.q_{\pm} = \frac{B \pm \sqrt{B^{2} - 4D_{u}D_{v}\det\mathbf{J}}}{2D_{u}D_{v}},\quad\quad q_{-} < q < q_{+}.

El modo homogéneo k=0k = 0 permanece estable. También lo son los modos de número de onda suficientemente grande, porque en ellos la difusión domina la respuesta lineal. La inestabilidad aparece en escalas intermedias; esa selección de un intervalo finito permite al mecanismo producir patrones espaciales definidos. La fórmula no determina directamente si el patrón final será estacionario, aunque en el caso clásico de dos especies la rama que cruza a crecimiento suele hacerlo a través de un autovalor real.

Si Du=Dv=DD_{u} = D_{v} = D, se cumple A(q)=J−DqI\mathbf{A}(q) = \mathbf{J} - Dq\mathbf{I}. Cada autovalor cinético se desplaza hacia la izquierda en DqDq, por lo que un equilibrio local estable sigue siendo linealmente estable para toda perturbación espacial. Esta observación demuestra, dentro de las hipótesis adoptadas, por qué se requieren difusividades distintas. Eso no significa que cualquier sistema con difusión igual carezca de patrones, porque las fronteras forzadas, la cinética con varios atractores, el acoplamiento mecánico o el transporte no descrito por esos laplacianos pueden generar estructuras por mecanismos diferentes.

144.4. Cómo la difusión puede desestabilizar una cinética estable

La paradoja aparente viene de asociar la difusión solo con la homogeneización. En una sola sustancia, el transporte difusivo reduce gradientes sin alterar la estabilidad local por medio de otras variables. Con dos especies que interactúan, un desplazamiento espacial modifica tanto las concentraciones como las velocidades de reacción en cada lugar. Si una pequeña elevación de la primera especie promueve procesos que aumentan su concentración y la segunda responde con un efecto compensador, la rapidez relativa con que se desplazan ambas puede decidir si la compensación alcanza a actuar.

Imaginemos un mecanismo abstracto en que una especie favorece una respuesta local que, de mantenerse sin oposición, aumentaría su propia presencia. Otra limita esa acumulación y puede difundirse sobre distancias mayores durante el mismo intervalo. Si la influencia limitadora se redistribuye rápidamente, el efecto compensador local disminuye en algunos lugares mientras se distribuye por sus vecinos. Una diferencia inicial puede amplificarse cerca de donde surgió y quedar contenida en una escala mayor. Esta es una intuición útil de activación local e inhibición de mayor alcance, aunque el patrón exacto de signos en las reacciones y los tiempos característicos deben verificarse en cada modelo.

El argumento no requiere que una molécula decida su dirección de movimiento ni que exista una plantilla espacial ya dibujada. Las reacciones son locales y los flujos difusivos siguen gradientes de concentración. La estructura se debe a cómo interactúan las tasas de ambos procesos cuando aparece una pequeña diferencia. Matemáticamente, cada escala espacial cambia de forma desigual los elementos diagonales del jacobiano; para ciertos modos, esa alteración puede convertir una respuesta neta restauradora en una amplificadora.

El mecanismo depende de escalas. Si la perturbación varía muy lentamente en el espacio, la difusión apenas distingue las regiones y domina la cinética local estable. Si varía con enorme rapidez, la difusión elimina las diferencias antes de que puedan mantenerse. Entre ambos límites puede haber una ventana donde la combinación de acoplamientos favorezca el crecimiento. La banda inestable indica qué longitudes de onda podrían llegar a dominar la etapa temprana de formación de estructura.

La interpretación intuitiva tiene, sin embargo, sus límites. En tejidos vivos, las proteínas pueden difundirse de manera anisótropa, unirse a receptores, degradarse o ser transportadas activamente. Algunas señales viajan entre células sin desplazarse como un soluto libre. Dos mecanismos con patrones visuales similares pueden tener parámetros cinéticos y transportes radicalmente distintos. Por ello, observar manchas periódicas no prueba que exista una pareja simple de activador e inhibidor. Para que el modelo sea plausible hay que identificar sus variables con magnitudes medibles y someter a prueba las relaciones funcionales propuestas.

Incluso cuando se conoce bien una vía química, la relación entre una concentración y un rasgo visible suele pasar por procesos adicionales (expresión génica, propiedades mecánicas, cambios celulares). El patrón químico puede orientar una respuesta morfológica sin determinar todos sus detalles. Las ecuaciones de reacción–difusión formulan una hipótesis mecanística que puede verificarse mediante perturbaciones de parámetros.

144.5. El sistema de Schnakenberg como ejemplo calculable

Para examinar las condiciones sin depender únicamente de argumentos generales, consideremos el modelo adimensional de Schnakenberg

∂tu=a−u+u2v+DuΔu,∂tv=b−u2v+DvΔv,\partial_{t}u = a - u + u^{2}v + D_{u}\Delta u,\quad\quad\partial_{t}v = b - u^{2}v + D_{v}\Delta v,

con a,b>0a,b > 0 y Du,Dv>0D_{u},D_{v} > 0. Sus funciones de reacción describen una cinética idealizada; los parámetros adimensionales no deben leerse como tasas universales de una reacción biológica particular. Su ventaja es que permite obtener un equilibrio y una relación de dispersión explícitos. En el equilibrio se cumplen a−u+u2v=0a - u + u^{2}v = 0 y b−u2v=0b - u^{2}v = 0, de donde se obtiene

u∗=a+b,v∗=b(a+b)2.\boxed{u_{*} = a + b,\quad\quad v_{*} = \frac{b}{(a + b)^{2}}.}

El jacobiano de las reacciones, evaluado en ese estado, es

J=(−1+2u∗v∗u∗2−2u∗v∗−u∗2).\mathbf{J} = \begin{pmatrix} - 1 + 2u_{*}v_{*} & u_{*}^{2} \\ - 2u_{*}v_{*} & - u_{*}^{2} \end{pmatrix}.

Elijamos a=0.1a = 0.1 y b=0.9b = 0.9, con lo que u∗=1u_{*} = 1 y v∗=0.9v_{*} = 0.9. Resulta

J=(0.81−1.8−1),trJ=−0.2,det⁡J=1.\mathbf{J} = \begin{pmatrix} 0.8 & 1 \\ - 1.8 & - 1 \end{pmatrix},\quad\quad tr\mathbf{J} = - 0.2,\quad\quad\det\mathbf{J} = 1.

El equilibrio es localmente estable. Sus autovalores son −0.1±i0.99- 0.1 \pm i\sqrt{0.99}, de modo que ante pequeñas perturbaciones homogéneas hay una relajación oscilatoria amortiguada. Para investigar si la difusión cambia esa conclusión, asignemos Du=0.01D_{u} = 0.01 y Dv=1D_{v} = 1. En estas unidades adimensionales, la segunda especie difunde cien veces más rápido; la elección es deliberada y facilita exhibir el mecanismo.

El coeficiente de la condición de Turing vale B=Dvfu+Dugv=0.8−0.01=0.79B = D_{v}f_{u} + D_{u}g_{v} = 0.8 - 0.01 = 0.79. Como B2=0.6241B^{2} = 0.6241 supera a 4DuDvdet⁡J=0.044D_{u}D_{v}\det\mathbf{J} = 0.04, se satisfacen las condiciones. El polinomio característico espacial queda

det⁡A(q)=0.01q2−0.79q+1.\det\mathbf{A}(q) = 0.01q^{2} - 0.79q + 1.

Sus ceros están aproximadamente en q−=1.287q_{-} = 1.287 y q+=77.713q_{+} = 77.713. Los modos cuyo número de onda cumple 1.134<∣k∣<8.8161.134 < |k| < 8.816 son linealmente inestables en el continuo ideal. El modo uniforme k=0k = 0 sigue amortiguado, por construcción. También decaen las variaciones espaciales extremadamente finas, ya que los términos de difusión acaban dominando para números de onda grandes.

Si estudiamos un intervalo [0,π]\lbrack 0,\pi\rbrack con condiciones de Neumann, los números de onda permitidos son los enteros no negativos. El modo n=1n = 1 queda por debajo de la banda; los modos n=2,…,8n = 2,\ldots,8 quedan dentro; n=9n = 9 queda fuera. Al evaluar las raíces temporales encontramos para n=3n = 3 una tasa dominante aproximada de 0.5390.539 por unidad temporal, frente a 0.4280.428 para n=2n = 2 y 0.1320.132 para n=8n = 8. La evolución inicial favorece el crecimiento de determinadas escalas, aunque las amplitudes dependerán de cuánto esté presente cada modo en el estado inicial.

Este cálculo es una predicción de la linealización. Para conocer el resultado tardío debemos integrar el sistema completo, comprobar que las soluciones mantienen propiedades físicas apropiadas y distinguir un estado estacionario de una estructura que todavía evoluciona lentamente. Una simulación que termina mostrando tres o cuatro máximos puede ser compatible con la tasa de crecimiento de ciertos modos, pero su número final de máximos también depende de la interacción no lineal, el tamaño del dominio y los criterios con que se define un máximo.

Schnakenberg es estable con difusividades iguales y presenta una banda inestable entre k aproximadamente 1,134 y 8,816 cuando la segunda especie difunde cien veces más rápido.
Figura 144.1. Para los parámetros del manuscrito, la cinética homogénea es estable. Con Dᵥ = Dᵤ = 0,01, todos los autovalores se desplazan hacia tasas más negativas. Con Dᵥ = 1 aparece una banda 1,134 < k < 8,816. En L = π con Neumann, los modos n = 2,…,8 caen dentro. La figura establece una inestabilidad lineal, no la existencia ni la forma de un patrón persistente. Ampliar figura ↗

144.6. Del crecimiento inicial a la saturación no lineal

La linealización muestra por qué algunas diferencias espaciales se amplifican al comienzo. Sin embargo, una perturbación que crece como eσte^{\sigma t} y tiene σ>0\sigma > 0 no puede seguir siendo pequeña indefinidamente. Antes de que su amplitud alcance valores arbitrariamente grandes, los términos descartados en el desarrollo de Taylor vuelven a ser importantes. En una reacción química pueden cambiar las velocidades por agotamiento de reactivos, saturación enzimática o modificación de concentraciones que alteran las propiedades locales. La formación de un patrón persistente requiere estudiar cómo esos efectos limitan o reorganizan el crecimiento.

Cerca de ciertos umbrales, cuando domina una amplitud lentamente variable y existen las simetrías adecuadas, la dinámica puede aproximarse mediante una ecuación de amplitud. Para una representación real y espacialmente uniforme de la amplitud de un modo crítico, una forma normal elemental es

dAdt=μA−gA3.\frac{dA}{dt} = \mu A - gA^{3}.

Si g>0g > 0, el equilibrio A=0A = 0 pierde estabilidad cuando μ\mu cambia de signo y aparecen estados no nulos A=±μ/gA = \pm \sqrt{\mu/g} para μ>0\mu > 0. Estos valores son amplitudes en una aproximación reducida, distintas de las concentraciones químicas completas. Tampoco hay garantía universal de que la bifurcación sea supercrítica, porque el valor y el signo de gg se obtienen de la estructura no lineal del modelo; un término cúbico con signo opuesto puede requerir correcciones de orden superior y permitir histéresis o coexistencia de soluciones.

En un dominio extendido, una sola amplitud puede ser insuficiente. Puede haber varios números de onda inestables, orientaciones equivalentes y efectos de fase. Cuando dos modos crecen con tasas próximas, sus interacciones deciden si terminan coexistiendo o si uno reduce la amplitud del otro. En dos dimensiones, la posibilidad de combinar varias direcciones de onda influye en que aparezcan rayas, manchas o redes. Las simetrías y las resonancias entre vectores de onda condicionan esas selecciones, pero el signo de los coeficientes no lineales sigue siendo imprescindible.

El patrón tampoco tiene por qué ser estable durante tiempos largos. Algunas estructuras se desplazan, cambian de orientación, se fragmentan o alternan entre configuraciones. Una imagen tomada en un instante puede sugerir estacionariedad donde hay una evolución lenta. Para certificar un estado final conviene seguir las amplitudes de Fourier, las velocidades de las interfaces y las tasas de cambio de las concentraciones durante un intervalo suficientemente largo comparado con los tiempos cinéticos y difusivos del sistema.

Hay además un fenómeno de selección dependiente de las condiciones iniciales. Una pequeña perturbación aleatoria puede estimular varios modos y conducir a un patrón; una perturbación localizada de amplitud grande puede llevar al sistema a otro estado de atracción. Incluso fuera de la ventana de inestabilidad lineal, la cinética no lineal puede sostener estructuras generadas por perturbaciones finitas. En tal situación conviene describir el mecanismo como una formación de patrón en un sistema de reacción–difusión, sin atribuirle automáticamente el criterio lineal clásico de Turing.

Krause y colaboradores (2024) describieron sistemas con un equilibrio homogéneo inestable a ciertos modos espaciales que, tras formar estructuras transitorias, terminan en otro equilibrio homogéneo estable. La inestabilidad lineal detecta el crecimiento inicial, pero de ella no se deduce que el patrón persista.

La amplitud del tercer modo crece inicialmente como predice la linealización y después se satura; el perfil numérico final muestra variaciones positivas de u y v.
Figura 144.2. Se integra la EDP de Schnakenberg con a = 0,1, b = 0,9, Dᵤ = 0,01 y Dᵥ = 1. La perturbación inicial tiene amplitud 0,001 en cos(3x), con proporción entre especies dada por el autovector creciente. La tasa temprana coincide con la predicción del operador discretizado; luego actúan las no linealidades. El perfil de t = 60 corresponde a 384 celdas. No se infiere estabilidad indefinida de esta instantánea. Ampliar figura ↗

144.7. Geometría, fronteras y selección de manchas o franjas

La misma cinética puede producir distribuciones diferentes cuando cambia el lugar donde ocurre. En un intervalo unidimensional, los modos de Neumann tienen perfiles cosenoidales; sobre un anillo son naturales los modos periódicos. Una lámina bidimensional añade nuevas orientaciones y permite combinar ondas planas. La curvatura de una superficie, su conectividad y las regiones donde el flujo está restringido modifican el espectro del laplaciano y, con él, las posibilidades de crecimiento.

Para una lámina rectangular aislada de lados LxL_{x} y LyL_{y}, los modos propios espaciales pueden escribirse como cos⁡(nπx/Lx)cos⁡(mπy/Ly)\cos\left( n\pi x/L_{x} \right)\cos\left( m\pi y/L_{y} \right), con n,m≥0n,m \geq 0. Su valor espectral es

qnm=π2(n2Lx2+m2Ly2).q_{nm} = \pi^{2}\left( \frac{n^{2}}{L_{x}^{2}} + \frac{m^{2}}{L_{y}^{2}} \right).

Si la banda inestable ocupa (q−,q+)\left( q_{-},q_{+} \right), solo los pares de índices que satisfagan q−<qnm<q+q_{-} < q_{nm} < q_{+} pueden crecer linealmente. Un rectángulo alargado quizá favorezca variaciones a lo largo de una dirección, mientras una región casi cuadrada puede permitir modos con orientaciones comparables. Si dos modos tienen el mismo valor espectral, sus combinaciones crean configuraciones que el análisis lineal considera equivalentes; las no linealidades y las imperfecciones reales pueden seleccionar una de ellas.

Las fronteras también pueden imponer concentraciones. En una condición de Dirichlet, una sustancia puede mantenerse a un nivel fijo en el borde mediante un reservorio externo. En ese caso, el estado homogéneo que usamos para derivar las condiciones de Turing quizá ya no satisfaga el problema físico, salvo que el valor impuesto coincida exactamente con el equilibrio. Hay que linealizar entonces alrededor de una solución que respete esas fronteras. Una condición de flujo impuesto, por su parte, puede producir gradientes incluso cuando el equilibrio uniforme es estable; el patrón observado contendrá un efecto de forzamiento que debe distinguirse del crecimiento espontáneo.

El tamaño absoluto del medio determina cuántas longitudes de onda caben en él. Si el dominio crece junto con un organismo, sus autovalores espaciales cambian con el tiempo. El patrón puede añadir nuevas regiones, sufrir divisiones o reajustarse. En otros casos, los mecanismos biológicos parecen mantener distancias relativas pese al crecimiento, lo que exige modelos con regulación adicional o parámetros dependientes del tamaño.

Estas observaciones también afectan a la comparación entre simulación y experimento. Un gráfico rectangular con bordes periódicos equivale topológicamente a una región que se conecta consigo misma; usarlo para representar una placa física con paredes impermeables cambia el problema. Las aparentes rayas que atraviesan un borde y reaparecen por el opuesto pueden ser coherentes en una simulación periódica, pero resultar imposibles en un experimento con reservorios o barreras materiales. La geometría debe especificarse al mismo nivel que las tasas de reacción.

144.8. Simular patrones sin confundir crecimiento físico con error numérico

Una simulación de reacción–difusión requiere resolver simultáneamente escalas temporales y espaciales que pueden ser muy diferentes. Para una malla unidimensional de celdas uniformes, una aproximación elemental del laplaciano es (ui+1−2ui+ui−1)/Δx2\left( u_{i + 1} - 2u_{i} + u_{i - 1} \right)/\Delta x^{2}. Si usamos un avance de Euler explícito, la primera especie se actualiza como

uin+1=uin+Δt[f(uin,vin)+Duui+1n−2uin+ui−1nΔx2],u_{i}^{n + 1} = u_{i}^{n} + \Delta t\left\lbrack f\left( u_{i}^{n},v_{i}^{n} \right) + D_{u}\frac{u_{i + 1}^{n} - 2u_{i}^{n} + u_{i - 1}^{n}}{\Delta x^{2}} \right\rbrack,

con una fórmula análoga para vv. Esta discretización es fácil de inspeccionar, pero puede ser inestable si Δt\Delta t es demasiado grande. Para difusión pura en una dimensión, la restricción clásica es DΔt/Δx2≤1/2D\Delta t/\Delta x^{2} \leq 1/2. En el sistema acoplado hay que considerar además la mayor difusividad, las tasas de reacción y las regiones del espacio de estados por las que pasa la solución.

Los métodos implícitos o semimplícitos permiten superar algunas restricciones del tratamiento explícito, a costa de resolver sistemas algebraicos o separar cuidadosamente los términos. En cualquier caso, el algoritmo elegido debe preservar las cantidades y restricciones relevantes. Si las ecuaciones describen concentraciones no negativas, una integración que genera valores negativos requiere un examen serio. A veces se debe a un paso demasiado grande; otras, a un esquema que no preserva positividad. Recortar artificialmente los valores negativos puede alterar la masa y las tasas de reacción, por lo que no debe presentarse como una corrección inocua.

Para observar la desestabilización de un equilibrio, la simulación suele empezar con el estado homogéneo más una perturbación pequeña generada con una semilla reproducible. Si la perturbación solo excita modos fuera de la banda inestable, el sistema puede permanecer cercano al equilibrio durante el intervalo de observación. El ruido numérico podría acabar introduciendo otros modos, pero no es una preparación experimental bien definida. Es mejor registrar la forma, la amplitud y el espectro de la perturbación inicial.

La verificación puede organizarse mediante pruebas independientes. Primero, comprobar que el equilibrio calculado anula las reacciones. Segundo, comparar la tasa inicial de crecimiento de un modo aislado con el autovalor de la matriz J−qD\mathbf{J} - q\mathbf{D}. Tercero, refinar por separado Δx\Delta x y Δt\Delta t y observar si las tasas convergen a las predicciones del continuo. Cuarto, verificar que el tratamiento de las fronteras produzca los flujos prescritos. Estas pruebas detectan errores que pueden pasar inadvertidos en una imagen atractiva de manchas o rayas.

Para decidir si el patrón tardío es confiable hacen falta además pruebas de duración y robustez. Conviene repetir simulaciones con perturbaciones iniciales independientes, cuantificar la distancia entre configuraciones, examinar el espectro de Fourier o del laplaciano y comprobar si el perfil continúa cambiando cuando se duplica el tiempo de integración. Si dos mallas generan dibujos visualmente parecidos pero seleccionan números de onda distintos, el resultado todavía no está numéricamente resuelto.

144.9. Qué puede y qué no puede explicar la morfogénesis por reacción–difusión

La propuesta de Turing abrió una vía para estudiar cómo una estructura espacial puede resultar de reglas dinámicas locales. La investigación posterior extendió estos modelos al desarrollo de tejidos, a reacciones químicas y a otros sistemas alejados del equilibrio. En determinados organismos y procesos hay pruebas experimentales compatibles con mecanismos de tipo reacción–difusión; en otros casos, los modelos son todavía una explicación alternativa a mecanismos mecánicos, genéticos o de señalización espacial. Con aplicaciones tan diversas, hay que separar la evidencia específica de las semejanzas visuales.

Una predicción útil es cómo cambia la longitud de onda al modificar una tasa de transporte. Como kk aparece multiplicado por difusividades en la matriz linealizada, cambiar esas tasas puede desplazar o eliminar la banda inestable. Sin embargo, la relación entre un parámetro molecular y la difusividad efectiva del modelo no siempre es simple. La unión reversible a receptores, la circulación de material por una matriz extracelular o la degradación durante el transporte pueden introducir escalas adicionales. Un experimento que modifica una proteína podría alterar simultáneamente su velocidad de movimiento y las reacciones en que participa.

La interpretación exige medir más que la forma final. Idealmente deberían conocerse las variaciones temporales de las concentraciones, sus velocidades de transporte, las respuestas ante pequeñas perturbaciones y las condiciones de contorno. La detección de una banda inicial de crecimiento espacial sería evidencia más específica que la simple observación de rayas. También importan las pruebas de causalidad. Si se modifica un coeficiente relevante, ¿cambia el umbral o la escala en el sentido previsto? ¿Puede una perturbación controlada desencadenar o suprimir el patrón? ¿Se mantienen estas relaciones al cambiar el tamaño de la región?

La revisión histórica de Turing obliga a ser cuidadosos con las etiquetas. Una reconstrucción moderna publicada por John J. Tyson en 2026 recuerda que ciertos patrones pueden emerger mediante perturbaciones de amplitud finita incluso donde las condiciones de inestabilidad lineal clásica no se cumplen para los parámetros analizados. Por ello conviene distinguir, en la literatura y en las simulaciones, la inestabilidad lineal inducida por difusión de la formación de patrones en un sistema que incluye reacción y difusión. La primera puede iniciar la segunda, pero no es su único mecanismo de inicio.

El modelo mínimo supone, además, un medio continuo y bien mezclado a escalas pequeñas. En tejidos formados por un número finito de células, las fluctuaciones y el transporte entre compartimentos pueden importar. La mecánica celular modifica distancias y superficies; la expresión génica incorpora demoras y regulación no lineal. Algunos sistemas muestran patrones controlados por fuerzas elásticas o por crecimiento diferencial. Los componentes adicionales complementan el análisis de reacción–difusión, porque permiten evaluar qué parte de la estructura puede atribuirse a sus mecanismos y cuál requiere procesos que el sistema de dos ecuaciones omite.

La variedad de imágenes que un modelo puede generar no mide su potencia explicativa. Con suficientes parámetros, muchos modelos producen manchas, rayas y laberintos parecidos. Una explicación científica gana fuerza cuando permite anticipar cómo responde el fenómeno a manipulaciones cuantitativas y cuando se comparan explícitamente alternativas. La formación de un patrón es una observación; identificar su mecanismo es un problema de inferencia y validación que exige más información que la apariencia del resultado.

144.10. Del patrón observado a una hipótesis contrastable

Imaginemos que una lámina de tejido cultivado desarrolla bandas de una señal fluorescente. Después de una preparación inicialmente casi uniforme, aparecen variaciones regulares que aumentan de contraste. Una primera inspección permite medir su separación media, pero no establecer por qué se forman. Si el tejido crece durante la observación, la distancia entre bandas podría cambiar simplemente por deformación geométrica. Si la iluminación es desigual, algunas regiones podrían parecer más concentradas sin que exista una diferencia química real. Antes de ajustar una ecuación conviene definir qué magnitud cuantifica la fluorescencia y cómo se corrige su relación con la concentración.

Supongamos, para construir un análisis conceptual, que podemos medir dos sustancias relevantes, estimar sus velocidades locales de reacción y obtener un intervalo plausible para sus difusividades. Con estos datos hallaríamos un equilibrio homogéneo y calcularíamos el jacobiano cinético. Comprobaríamos primero si la dinámica local es estable. Si lo es, construiríamos A(q)=J−qD\mathbf{A}(q) = \mathbf{J} - q\mathbf{D}, hallaríamos la banda inestable y estudiaríamos qué números de onda admite la geometría real del tejido. Este procedimiento predice, antes de simular un dibujo completo, qué longitudes de onda deberían aumentar al comienzo y con qué tasas aproximadas.

El siguiente paso sería introducir una perturbación de amplitud conocida o usar pequeñas fluctuaciones espontáneas medidas con resolución suficiente. Al estimar las amplitudes espaciales durante la etapa temprana, podríamos contrastar la pendiente temporal de cada modo con la relación de dispersión calculada. Si el tejido está lejos de un estado estacionario o sus parámetros cambian continuamente durante el experimento, la hipótesis de una relación de dispersión fija podría fallar. En ese caso sería necesario tratar el estado de referencia y el operador linealizado como dependientes del tiempo y explicitar el alcance de cualquier aproximación cuasiestacionaria.

Para interpretar el patrón tardío necesitaríamos una segunda comparación. Integraríamos las ecuaciones no lineales con las mismas fronteras y una condición inicial representativa, comprobaríamos convergencia numérica y observaríamos si se reproducen la estabilidad, la escala y la evolución temporal de las bandas. Una coincidencia razonable sería evidencia favorable, aunque seguiría abierta la posibilidad de modelos alternativos. La estrategia más convincente es contrastar predicciones que distintos mecanismos no compartan, como sus respuestas ante cambios controlados de transporte, tiempos de reacción o geometría.

Las ecuaciones de reacción–difusión muestran cómo un equilibrio estable puede perder estabilidad espacial al acoplar procesos locales con transporte. La idea es rigurosa, pero tiene condiciones de aplicación precisas. El estado uniforme debe existir, el espectro espacial tiene que contener modos inestables y el comportamiento tardío requiere la dinámica no lineal. Conocer estas condiciones evita identificar como «autoorganización» cualquier estructura regular y permite convertir una imagen sugerente en una hipótesis cuantitativa que indique qué medir, qué perturbación aplicar y qué resultado obligaría a revisar el modelo.

El interés general del mecanismo está en la aparición de una escala espacial a partir de interacciones distribuidas. La difusión conecta lugares distintos; la cinética decide cómo responde cada lugar a esas influencias. Algunas combinaciones amortiguan todas las diferencias; otras permiten que crezca una banda delimitada del espectro. La regularidad visible es una consecuencia posible de esa selección dinámica. Su estudio exige pasar continuamente entre ecuaciones, observaciones y pruebas de estabilidad sin confundir una posibilidad matemática con una explicación física o biológica validada.

Referencias

  • Turing, A. M. (1952). The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society B, 237(641), 37–72. https://doi.org/10.1098/rstb.1952.0012
  • Murray, J. D. (2003). Mathematical Biology II: Spatial Models and Biomedical Applications (3.ª ed.). Springer.
  • Cross, M. C., y Hohenberg, P. C. (1993). Pattern formation outside of equilibrium. Reviews of Modern Physics, 65(3), 851–1112. https://doi.org/10.1103/RevModPhys.65.851
  • Meinhardt, H. (1982). Models of Biological Pattern Formation. Academic Press.
  • Villar-Sepúlveda, E., y Champneys, A. R. (2023). General conditions for Turing and wave instabilities in reaction–diffusion systems. Journal of Mathematical Biology, 86, 39. https://doi.org/10.1007/s00285-023-01870-3
  • Tyson, J. J. (2026). Revisiting Turing’s chemical basis of morphogenesis. Bulletin of Mathematical Biology, 88, artículo 76. https://doi.org/10.1007/s11538-026-01629-z
  • Krause, A. L., Gaffney, E. A., Jewell, T. J., Klika, V., y Walker, B. J. (2024). Turing instabilities are not enough to ensure pattern formation. Bulletin of Mathematical Biology, 86, artículo 21. https://doi.org/10.1007/s11538-023-01250-4