Morfogenesi & Pattern di Turing

Nell'articolo del 1952, Turing esplora dei sistemi reazione-diffusione, con l'obiettivo di spiegare i meccanismi che determinano la formazione di strutture anatomiche in un organismo, come — ad esempio — la pigmentazione dei tessuti. Tali processi, detti di morfogenesi, possono essere modellizzati da un sistema a reazione-reazione, in cui sostanze chimiche — dette morfogeni — interagiscono tra loro e con l'ambiente. Partendo da uno stato di quasi-omogeneità, è possibile assistere a particolari evoluzioni del sistema, causate dall'instabilità dell'equilibrio.

Turing si occupa non solamente degli aspetti matematici, ma ricerca anche giustificazioni biologiche e chimiche per le osservazioni da lui effettuate. Questa relazione si concentrerà tuttavia unicamente sulla componente matematica del problema.

L'anello cellulare

Per illustrare i principi della morfogenesi, consideriamo - seguendo strettamente l'articolo del 1952 di Turing - una particolare configurazione: un anello di cellule simili. Pur non essendo riconducibile ad un corrispettivo biologico, questa configurazione facilita l'analisi matematica del problema e la teoria ad essa applicata può essere eventualmente estesa ad altri casi più complessi.

Per semplicità, consideriamo per ora un sistema in cui solo due morfogeni, $X$ e $Y$, interagiscono tra loro. Ogni cellula scambia i morfogeni con le cellule adiacenti per diffusione. Inoltre, in ogni cellula i morfogeni reagiscono tra loro ed eventualmente con catalizzatori e altre molecole, la cui disponibilità - per semplicità - si presuppone illimitata . Il prodotto di tale reazione è un aumento o una diminuzione dei morfogeni, che dipende dalla concentrazione degli stessi.

Risulta così ovvio il motivo per cui tale sistema è noto come reazione-diffusione.

Si suppone infine che, all'instante $t = 0$, le concentrazioni di $X$ e $Y$ non siano omogenee , a causa di una perturbazione che verrà tuttavia ignorata per $t > 0$.

Il sistema può pertanto essere descritto dalle seguenti $2N$ equazioni differenziali: \[ \begin{cases} \dfrac{dX_r}{dt} = \underbrace{f(X_r, Y_r)}_{\text{reazione}} + \underbrace{\mu (X_{r + 1} - 2X_r + X_{r-1})}_{\text{diffusione}} \\[1em] \dfrac{dY_r}{dt} = g(X_r, Y_r) + \nu (Y_{r + 1} - 2Y_r + Y_{r-1}) \end{cases} \qquad r = 1, ..., N \] dove:

• $N \in \mathbb{N} \;$ è il numero di cellule nell'anello;
• $X_r, Y_r \;$ indicano rispettivamente la concentrazione dei morfogeni $X$ e $Y$ nella cellula $r$;
• $f(X, Y) \, , \; g(X, Y)$ modellizzano la reazione dei morfogeni.
• $\mu \, \; \nu$ sono i coefficienti di diffusione rispettivamente dei morfogeni $X$ e $Y$.

Si noti che, data la struttura anulare, la cellula $X_N$ corrisponde alla $X_0$ e, più in generale $X_{r_0} \equiv X_{N - r_0}$.

Osservazioni

• Una cellula isolata ha equilibrio $X = h \, , \; Y = k$ quando $f(h, k) = 0 = g(h, k)$.
• Il sistema è pertanto in equilibrio se $X_r = h \, , \; Y_r = k \;\: \forall r = 1, ..., N \,$.
• Supponendo di non essere troppo distanti dall'equilibrio, è possibile scrivere $X_r = h + x_r$, $\, Y_r = k + y_r \,$ e approssimare $f(x + h, y + k) \approx ax + by \, , \; g(x + h, y + k) \approx cx + dy \,$.

Il sistema può quindi essere riscritto nella seguente forma: \[ \begin{cases} \dfrac{dx_r}{dt} = a x_r + b y_r + \mu (x_{r + 1} - 2x_r + x_{r-1}) \\[.5em] \dfrac{dy_r}{dt} = c x_r + d y_r + \nu (y_{r + 1} - 2y_r + y_{r-1}) \end{cases} \]

Con il seguente cambio di variabili \[ \begin{cases} x_r = \displaystyle \sum_{s = 0}^{N - 1} e^{\large \frac{2 \pi i r s}{N}} \xi_s \\[.5em] y_r = \displaystyle \sum_{s = 0}^{N - 1} e^{\large \frac{2 \pi i r s}{N}} \eta_s \\[.5em] \end{cases} \quad \Longrightarrow \quad \begin{cases} \xi_r = \dfrac{1}{N} \displaystyle \sum_{s = 1}^{N} e^{\large \frac{2 \pi i r s}{N}} x_s \\[.5em] \eta_r = \dfrac{1}{N} \displaystyle \sum_{s = 1}^{N} e^{\large \frac{2 \pi i r s}{N}} y_s \\[.5em] \end{cases} \] il sistema può inoltre così ricondotto in una forma lineare con coefficienti costanti: \[ \begin{cases} \dfrac{d\xi_s}{dt} = \Big(a - 4 \mu \sin^2 \dfrac{\pi s}{N} \Big) \xi_s + b \eta_s \\[.5em] \dfrac{d\eta_s}{dt} = c \xi_s + \Big(d - 4 \nu \sin^2 \dfrac{\pi s}{N} \Big) \eta_s \end{cases} \]

Soluzioni del sistema discreto

Detti $p_s$ e $p_s^\prime$ gli autovalori della matrice $ \scriptsize \begin{pmatrix} a - 4 \mu \sin^2 \dfrac{\pi s}{N} & \small b \\ \small c & d - 4 \nu \sin^2 \dfrac{\pi s}{N} \end{pmatrix} $, che pertanto soddisfano l'equazione \begin{equation} \Big(p - a + 4 \mu \sin^2 \dfrac{\pi s}{N} \Big) \Big(p - d + 4 \nu \sin^2 \dfrac{\pi s}{N} \Big) = bc \, , \end{equation} e $A_s, \, B_s, \, C_s, \, D_s$ le componenti degli autovettori associati, che soddisfano \begin{equation} \begin{cases} A_s \Big(p_s - a + 4 \mu \sin^2 \dfrac{\pi s}{N} \Big) = b C_s \\[.5em] B_s \Big(p_s - a + 4 \mu \sin^2 \dfrac{\pi s}{N} \Big) = bD_s \, , \end{cases} \end{equation} le soluzioni delle equazioni differenziali risultano essere \[ \begin{cases} \xi_s = A_s e^{p_s t} + B_s e^{p_s^\prime t} \\[.5em] \eta_s = C_s e^{p_s t} + D_s e^{p_s^\prime t} \, , \end{cases} \] ovvero, in termini di $X$ e $Y$, \begin{equation} \begin{cases}\label{eq:discrete-sol} X_r = h + \displaystyle \sum_{s = 1}^{N} \big(A_s e^{p_s t} + B_s e^{p_s^\prime t} \big) e^{\Large \frac{2 \pi i r s}{N}} \\[.5em] Y_r = k + \displaystyle \sum_{s = 1}^{N} \big(C_s e^{p_s t} + D_s e^{p_s^\prime t} \big) e^{\Large \frac{2 \pi i r s}{N}} \, . \end{cases} \end{equation}

Soluzioni del sistema continuo

Il modello può essere esteso al caso continuo, dove $\rho$ è il raggio dell'anello cellulare e $\theta$ descrive la posizione di un punto sull'anello. In tal caso, i coefficienti di diffusione sono $\mu^\prime \,$ e $\nu^\prime$, con $\mu = \mu^\prime \Big(\dfrac{N}{2 \pi \rho} \Big)^2 \,$ e, analogamente, $\nu = \nu^\prime \Big(\dfrac{N}{2 \pi \rho} \Big)^2$.

Le equazioni differenziali che descrivono il modello diventano \begin{equation} \begin{cases} \dfrac{\partial X_\theta}{\partial t} = a(X_\theta - h) + b(Y_\theta - k) + \dfrac{\mu^\prime}{\rho^2} \dfrac{\partial^2 X_\theta}{\partial \theta^2} \\[.5em] \dfrac{\partial Y_\theta}{\partial t} = c(X_\theta - h) + d(Y_\theta - k) + \dfrac{\nu^\prime}{\rho^2} \dfrac{\partial^2 X_\theta}{\partial \theta^2} \end{cases} \end{equation} con soluzioni date da \begin{equation} \begin{cases} X_\theta = h + \displaystyle \sum_{s = -\infty}^{+\infty} \big(A_s e^{p_s t} + B_s e^{p_s^\prime t} \big) e^{\large i \theta s} \\[.5em] Y_\theta = k + \displaystyle \sum_{s = -\infty}^{+\infty} \big(C_s e^{p_s t} + D_s e^{p_s^\prime t} \big) e^{\large i \theta s} \, . \end{cases} \end{equation} con $p_s$ e $p_s^\prime$ radici di $\Big(p - a + \dfrac{\mu^\prime s^2}{\rho^2} \Big) \Big(p - d + \dfrac{\nu^\prime s^2}{\rho^2} \Big) = bc \,$, mentre $A_s, \; B_s, \; C_s, \; D_s$ soddisfano \begin{equation} \begin{cases} A_s \Big(p_s - a + \dfrac{\mu^\prime s^2}{\rho^2} \Big) = b C_s \\[.5em] B_s \Big(p_s - a + \dfrac{\mu^\prime s^2}{\rho^2} \Big) = b D_s \end{cases} \end{equation}

Comportamento delle soluzioni

Osservazioni

• Data la \[ \small \begin{cases} X_r = h + \displaystyle \sum_{s = 1}^{N} \big(A_s e^{p_s t} + B_s e^{p_s^\prime t} \big) e^{\Large \frac{2 \pi i r s}{N}} \\[.5em] Y_r = k + \displaystyle \sum_{s = 1}^{N} \big(C_s e^{p_s t} + D_s e^{p_s^\prime t} \big) e^{\Large \frac{2 \pi i r s}{N}} \, . \end{cases} \] soluzione del sistema discreto , siamo interessati alla sua parte reale.
• Dopo un certo lasso di tempo, la soluzione risulterà dominata dai termini per cui i corrispondenti $p_s$ hanno la parte reale maggiore.
• Se $p_{s_0}$ è tra questi, allora anche $p_{N - s_0}$ lo è, perché $\: \sin^2 {\large \frac{\pi (N - s_0)}{N}} = \sin^2 {\large \frac{\pi s_0}{N}}$.
• Vale $\, \textrm{Re} (p_{s_0}) = \textrm{Re} (p_{s_0}^\prime) = \textrm{Re} (p_{\small N - s_0}) = \textrm{Re} (p_{\small N - s_0}^\prime) \,$.

Per il sistema discreto, si distinguono così due macro-casi:

• Stazionario: se $\, p_{s_0} = I \,$ è reale, allora $X_r - h \,$ e $\, Y_r - k$ tendono asintoticamente alla forma \begin{equation} \begin{cases} \underbrace{X_r - h}_{ \begin{gathered} \text{distanza} \\[-.3em] \text{da'eq.} \end{gathered} } = 2 \, e^{I t} \, \text{Re}\Big( A_{s_0} \, e^{\large \frac{2 \pi i s_0 r}{N}} \Big) \\[.5em] Y_r - k = 2 \, e^{I t} \, \text{Re}\Big( C_{s_0} \, e^{\large \frac{2 \pi i s_0 r}{N}} \Big) \end{cases} \end{equation}
• Oscillatorio: se $\, p_{s_0} = I + i \omega \,$ è complesso, con $I, \omega \in \mathbb{R}$, allora $X_r - h \,$ e $\, Y_r - k$ tendono alla forma \begin{equation} \begin{cases} X_r - h = 2 \, e^{I t} \, \text{Re}\Big( A_{s_0} \, e^{\large \frac{2 \pi i s_0 r}{N} + i \omega t} + A_{\small N - s_0} \, e^{\large \frac{2 \pi i s_0 r}{N} - i \omega t} \Big) \\[.5em] Y_r - k = 2 \, e^{I t} \, \text{Re}\Big( C_{s_0} \, e^{\large \frac{2 \pi i s_0 r}{N} + i \omega t} + C_{\small N - s_0} \, e^{\large \frac{2 \pi i s_0 r}{N} - i \omega t} \Big) \end{cases} \end{equation}

Nel caso continuo, questi corrispondono ai seguenti macro-casi:

• Stazionario: se $\, p_{s_0} = I \,$ è reale, allora \begin{equation} \begin{cases} x_\theta - h = 2 \, \text{Re}\Big( A_{s_0} \, e^{i \theta s_0 + I t} \Big) \\[.5em] y_\theta - k = 2 \, \text{Re}\Big( C_{s_0} \, e^{i \theta s_0 + I t} \Big) \end{cases} \end{equation}
• Oscillatorio: se $\, p_{s_0} = I + i \omega \,$ è complesso, allora \begin{equation} \begin{cases} x_\theta - h = 2 \, e^{I t} \, \text{Re}\Big( A_{s_0} \, e^{i \theta s_0 + i \omega t} + A_{\small N - s_0} \, e^{i \theta s_0 - i \omega t} \Big) \\[.5em] y_\theta - k = 2 \, e^{I t} \, \text{Re}\Big( C_{s_0} \, e^{i \theta s + i \omega t} + C_{\small N - s_0} \, e^{i \theta s_0 - i \omega t} \Big) \end{cases} \end{equation}

Osservazioni

• Nel caso stazionario, si possono appunto osservare onde stazionarie aventi $s_0$ creste e avvallamenti lungo l'anello. Essendo i coefficienti $A_{s_0}$ e $C_{s_0}$ in

  Caso discreto: \[ \small A_{s_0} \Big(p_{s_0} - a + 4 \mu \sin^2 \dfrac{\pi s_0}{N} \Big) = b C_{s_0} \, \]

  Caso continuo:
  \[ \small A_{s_0} \Big(p_{s_0} - a + \dfrac{\mu^\prime s_0^2}{\rho^2} \Big) = b C_{s_0} \, \]

  relazione tra loro , il pattern di un morfogeno determina quello dell'altro. Inoltre, quando il sistema è instabile, ovvero $I > 0$, le onde divento via via più pronunciate.
• Nel caso oscillatorio, si formano due onde che si muovono in direzioni opposte lungo l'anello, con frequenza $\, \dfrac{\omega}{2 \pi}$.

La lunghezza d'onda dipende non solo dai parametri del sistema, ma - nel caso continuo - anche dalla dimensione dell'anello, poiché ne dovrà essere un sottomultiplo. Proprio sulla base della lunghezza d'onda è quindi possibile differenziare ulteriori sottocasi, per un totale di sei, di cui due ottenibili solamente con tre o più morfogeni. Ci limiteremo tuttavia all'analisi dei sistemi con due morfogeni.

Caso stazionario, lunghezza d'onda $\gg 1$

Il contenuto di tutte le cellule è identico e non si osservano flussi evidenti causati dalla diffusione. Le cellule si comportano pertanto come se fossero isolate.

Questa configurazione si può ottenere quando $\, \mu = \nu = \frac{1}{4}$, $\, a = d$, $\, b = c = 1$.
In tal caso, $\, p_s = a - \sin^2 {\large \frac{\pi s}{N}} + 1 \,$ è reale ed massimo per $s = 0$. Per $\, a > 1$, il sistema è dunque instabile.

Nota: le ondulazioni visibili nel grafico di tempo non sono direttamente prodotti dall'evoluzione del sistema, ma sono causati dal rumore introdotto nella configurazione iniziale, che viene via via amplificato.

Durante le mie esplorazioni, ho notato che, per la configurazione $\, a = d = -1$, il sistema sembra tendere asintoticamente ad un equilibrio. Particolarmente interessante è il fatto che l'equilibrio sembra essere raggiunto indipendentemente dal rumore introdotto iniziale nel sistema.

Introducendo invece un rumore continuo nel tempo, il sistema sembra stabilizzarsi, ma - come ci si aspetterebbe - non raggiunge uno stato di equilibrio definitivo come nel caso precedente.

Studiando il caso continuo, per un anello di raggio $\, \rho \,$ sufficientemente grande da poter accomodare ogni tipo di onda, è possibile ricavare condizioni più precise che caratterizzano i vari sottocasi. Detto $\, U_s = \dfrac{s^2}{\rho^2} = \Big( \dfrac{2 \pi}{\lambda} \Big)^ \,2$, dove $\, \lambda \,$ indica la lunghezza d'onda, si può scrivere $\, p_s \,$ in funzione di $\, U_s$: \[ p_s = \dfrac{a + d}{2} - \dfrac{\mu^\prime + \nu^\prime}{2} U_s \pm \sqrt{\Big(\dfrac{\mu^\prime - \nu^\prime}{2} U_s + \dfrac{d - a}{2} \Big)^2 + bc} \, . \] La parte reale di $\, p \,$ rappresenta l'instabilità delle onde di lunghezza d'onda $\, \lambda = \dfrac{2 \pi} {\sqrt{U_s}}$. La componente dominante della soluzione è data dal massimo di $\, \text{Re}(p_s)$, verificato per $\, U_s \to 0^+$, oppure per $\, \, U_s \to +\infty$, oppure ancora per un punto intermedio, \[ U_s = \Big(a - d + \dfrac{\mu^\prime + \nu^\prime}{\sqrt{\mu^\prime \nu^\prime}} \sqrt{-bc} \Big) \dfrac{1}{\mu^\prime - \nu^\prime} \, , \] da cui si ricava \[ p_s = \big(d \mu^\prime - a \nu^\prime - 2 \sqrt{\mu^\prime \nu^\prime} \big) \dfrac{\sqrt{-bc}}{\mu^\prime - \mu^\prime} \, . \]

Per garantire la stazionarietà dell'onda dominante è necessario innanzitutto imporre che imporre che \[ \small p_{\large s} = \dfrac{a + d}{2} - \dfrac{\mu^\prime + \nu^\prime}{2} U_{\large s} \pm \sqrt{\Big(\dfrac{\mu^\prime - \nu^\prime}{2} U_{\large s} + \dfrac{d - a}{2} \Big)^2 + bc} \] $\, p_s \,$ sia reale, cioè \[ \Big(\dfrac{\mu^\prime - \nu^\prime}{2} U_{\large s} + \dfrac{d - a}{2} \Big)^2 + bc \geq 0 \, . \] Inoltre, affinché $\, \lambda = {\large \frac{2 \pi}{\sqrt{U_s}}} \gg 1 \,$ (o meglio $\, \lambda \to +\infty$) è possibile porre $\, U_s \to 0^+ \,$. Combinando tutte condizioni e il fatto che non debbano esistere $U_s$ intermedi per cui $\, p_s \,$ è massimo, si ottiene che il caso desiderato è dato da $\, bc \gt 0 \,$ oppure da $\, bc \lt 0 \,$ e \[ \dfrac{d - a}{\sqrt{-bc}} \gt \dfrac{\mu^\prime + \nu^\prime}{\sqrt{\mu^\prime \nu^\prime}} \quad \lor \quad \dfrac{d - a}{\sqrt{-bc}} \lt -2 \, . \] L'instabilità della soluzione è verificata per $\, p_s > 0$, ovvero $\, a + d > 0 \,$ oppure $\, bc \gt ad$.

Caso oscillatorio, lunghezza d'onda $\gg 1$

Come per l'analogo stazionario, le cellule si comportano come se fossero isolate. In tal caso, tuttavia, l'allontanamento dall'equilibrio segue un moto oscillatorio.

Questa configurazione si può ottenere quando $\, \mu = \nu = \frac{1}{4}$, $\, a = d$, $\, b = -c = 1$.
In tal caso, $\, p_s = a - \sin^2 \frac{\pi s}{N} \pm i \,$ è complesso e la parte reale è massima per $\, s = 0$.

Come per il caso stazionario, ho individuato la configurazione $\, a = d = -1 \,$ per cui il sistema sembra tendere asintoticamente all'equilibrio, indipendentemente dal rumore introdotto iniziale nel sistema.

Un rumore continuo nel tempo non impedisce al sistema di stabilizzarsi, ma - come previsto - non si verifica l'avvicinamento asintotico all'equilibrio.

Il caso oscillatorio si verifica per \[ \small p_{\large s} = \dfrac{a + d}{2} - \dfrac{\mu^\prime + \nu^\prime}{2} U_{\large s} \pm \sqrt{\Big(\dfrac{\mu^\prime - \nu^\prime}{2} U_{\large s} + \dfrac{d - a}{2} \Big)^2 + bc} \] $\, p_s \,$ complesso, pertanto è necessario imporre \[ \Big(\dfrac{\mu^\prime - \nu^\prime}{2} U_{\large s} + \dfrac{d - a}{2} \Big)^2 + bc \lt 0 \, . \] Inoltre, come per il caso stazionario, $\lambda \to +\infty \,$ è data da $\, U \to 0^+$. Tuttavia, La mia ipotesi è che la condizione derivi dalle disuguaglianze che caratterizzano i casi successivi, in particolare quello delle onde stazionarie con lunghezze d'onda finita. Turing raffina la disuguaglianza risultate , ottenendo \[ bc \lt 0 \quad \; \land \quad -2 \lt \dfrac{d - a}{\sqrt{-bc}} \lt \dfrac{4 \sqrt{\mu^\prime \nu^\prime}}{\mu^\prime + \nu^\prime} \, . \] L'instabilità è data banalmente da $\, a + d \gt 0$.

Caso stazionario, lunghezza d'onda breve

In questo caso, le onde risultanti sono stazionarie, ma la lunghezza d'onda estremamente breve genera pattern frastagliati.
Le configurazioni che seguono sono generate da $\, \mu = 1, \; \nu = 0$, $b = - c = 1$, $a = J_1 - 1$, $d = J_1 \,$ con $\, J_1 = {\large \frac{1}{4}}$. In questo caso, \[ p_s = J_1 - \dfrac{1}{2} - 2 \sin^2 \dfrac{\pi s}{N} + \sqrt{\Big( 2 \sin^2 \dfrac{\pi s}{N} + \dfrac{1}{2} \Big)^2 - 1} \, \] ed è massimo quando $\, \sin^2 \dfrac{\pi s}{N} \,$ è massimo, cioè $\, s = \dfrac{1}{2} N \;\, \lor \;\, s = \dfrac{1}{2} (N - 1) \,$.

In generale, il contenuto di ogni cellula tenderà ad essere quanto più diverso possibile da quello delle cellule adiacenti. Come si può osservare, per un numero di cellule sufficientemente grande, nonostante la generale frastagliatura, esistono punti in cui la differenza tra il contenuto di due cellule adiacenti è contenuta, mentre altri è più marcata. In caso di instabilità, la differenza crescerà con l'avanzare del tempo.

Esistono configurazioni dei parametri per cui il sistema sembra raggiungere un equilibrio.

Ad esempio, il grafico seguente è generato da $\, J_1 = 0.21$.

Per $\, J_1 = 0.20$, le onde sembrano via via appiattirsi, tendendo asintoticamente ad uno stato di equilibrio.

Quando il numero di cellule $N$ è contenuto, si possono osservare due differenti fenomeni. Per $\, N \,$ pari, il contenuto di ogni cellula $r$ tenderà ad essere simile a quello delle cellule $\, r + 2 \,$ e $\, r - 2$, ma quanto più dissimile da quello delle cellule adiacenti $\, r + 1 \,$ e $\, r - 1 \,$.

Per $N$ dispari, invece, ciò non è possibile e - tendenzialmente - la differenza tra cellule adiacenti sarà contenuta in un punto dell'anello e massima nel punto diametralmente opposto.

Il caso di onde stazionarie con lunghezza d'onda $\, \lambda = {\large \frac{2 \pi}{\sqrt{U}}} \to 0^+ \,$ (nel caso continuo) è verificato per \[ \small p_{\large s} = \dfrac{a + d}{2} - \dfrac{\mu^\prime + \nu^\prime}{2} U_{\large s} \pm \sqrt{\Big(\dfrac{\mu^\prime - \nu^\prime}{2} U_{\large s} + \dfrac{d - a}{2} \Big)^2 + bc} \] $\, p_s \,$ reale con \[ \small U_s = \Big(a - d + \dfrac{\mu^\prime + \nu^\prime}{\sqrt{\mu^\prime \nu^\prime}} \sqrt{-bc} \Big) \dfrac{1}{\mu^\prime - \nu^\prime} \] $\, U_s \to +\infty \,$ , ovvero $bc \lt 0 \,$ (per garantire l'esistenza di $U$) e, supponendo $\, \mu^\prime > \nu^\prime$, $\, \nu \to 0^+$. L'instabilità è data da $\, a + d \gt 0$.

Caso stazionario, lunghezza d'onda finita

Il seguente è il caso più rilevante dal punto di vista della biologia, poiché la concentrazione dei morfogeni nelle cellule è determinato da un'onda con lunghezza d'onda finita. Se si immagina che i morfogeni determino il colore di un tessuto, ad esempio, una configurazione di questo determinerebbe la formazione di motivi striati.

Una simile configurazione è generata, ad esempio, da $\, \mu = 1$, $\, \nu = \dfrac{1}{2}$, $\, a = J_2 -2$, $\, b = \dfrac{5}{2}$, $\, c = -\dfrac{5}{4}$, $\, d = J_2 + \dfrac{3}{2}$. Pertanto \[ \small \Big(p_s - a + 4 \mu \sin^2 \dfrac{\pi s}{N} \Big) \Big(p_s - d + 4 \nu \sin^2 \dfrac{\pi s}{N} \Big) = bc \] $\, p_s \,$ risolve \[ (p_s - J_2)^2 + \Big(\dfrac{1}{2} + \dfrac{3}{2} U_s \Big)(p - J_2) + \dfrac{1}{2} \Big( U_s - \dfrac{1}{2} \Big)^2 = 0 \, , \] che ha una soluzione $\, p_s = J_2 \,$ per $\, U_s = \dfrac{1}{2} \,$. Nel grafico che segue $\, p_s = J_2 = 0$.

Per $p_s = J_2 \lt 0$, l'onda si appiattisce, tendendo asintoticamente ad uno stato di equilibrio stabile.

Mentre, per $\, p_s = J_2 \gt 0 \,$ il sistema è instabile.

Il caso continuo, stazionario e con lunghezza d'onda finita è verificato per $\, bc \lt 0$, che garantisce l'esista di \[ \small U_s = \Big(a - d + \dfrac{\mu^\prime + \nu^\prime}{\sqrt{\mu^\prime \nu^\prime}} \sqrt{-bc} \Big) \dfrac{1}{\mu^\prime - \nu^\prime} \] $\, U_s \,$ , e da \[ \dfrac{4 \sqrt(\mu^\prime \nu^\prime)}{\mu^\prime + \nu^\prime} \lt \dfrac{d - a}{\sqrt{-bc}} \lt \dfrac{\mu^\prime + \nu^\prime}{\sqrt{\mu^\prime \nu^\prime}} \, . \] La prima disuguaglianza si ricava imponendo $\, p_{\text{max}} \gt p_0$. La seconda, invece, garantisce \[ \small U_s = \Big(a - d + \dfrac{\mu^\prime + \nu^\prime}{\sqrt{\mu^\prime \nu^\prime}} \sqrt{-bc} \Big) \dfrac{1}{\mu^\prime - \nu^\prime} \gt 0 \] $\, U_s > 0 \,$

Il sistema è instabile se \[ \small p_s = I = \big(d \mu^\prime - a \nu^\prime - 2 \sqrt{\mu^\prime \nu^\prime} \big) \dfrac{\sqrt{-bc}}{\mu^\prime - \mu^\prime} \gt 0 \] $\, p_{\text{max}} = I \gt 0 \,$ , cioè \[ \dfrac{1}{\sqrt{-bc}} \Bigg( d \, \sqrt{\dfrac{\mu^\prime}{\nu^\prime}} - a \, \sqrt{\dfrac{\nu^\prime}{\mu^\prime}} \Bigg) > 2 \, . \]

Altri casi

Due ulteriori casi oscillatori sono ottenibili considerando sistemi con tre morfogeni.

Caso oscillatorio, lunghezza d'onda finita

Caso oscillatorio, lunghezza d'onda breve

Modello generale per due morfogeni

Ulteriori considerazioni

Nel proseguo dell'articolo, Turing esplora più nel dettaglio ulteriori questioni matematiche, tra cui

• l'applicazione di un rumore continuo alla concentrazione dei morfogeni;
• reazioni $\, f \,$ e $\, g \,$ non lineari;
• diverse approssimazioni per il caso stazionario con lunghezza d'onda $\, \lambda \gg 1$.

Questi temi non saranno tuttavia affrontati in questa relazione.

Particolarmente interessante — dato l'anno della pubblicazione, il 1952 — è il fatto che Turing riconosca il ruolo che un computer potrebbe svolgere nella comprensione del problema, soprattutto per quanto riguarda il caso non lineare. Per questa ragione, ho pensato che l'inclusione dei grafici dinamici potesse essere uno strumento in sintonia con il pensiero dell'autore.

Tessuto cellulare

È possibile estendere il modello di reazione-diffusione a sistemi più complessi, come ad esempio una superficie piana composta da cellule. Nel caso che segue, si considera un griglia quadrata. Ogni unità della griglia (pixel) corrisponde ad una cellula. Ogni cellula confina con le otto cellule che la circondano (quattro con un lato in comune, quattro con un vertice comune) e pertanto la diffusione coinvolge queste unità. Dato il costo computazionale, per ottenere l'illusione di un tessuto esteso, alla griglia è stata applicata una topologia toroidale — detto banalmente, si comporta come la griglia del celebre gioco Snake.

Per modellizzare il sistema, è stata implementata la seguente formulazione: \begin{equation} \begin{cases} X_{i,j}^{\text{new}} = X_{i,j} + (\mu \nabla^2 X_{i,j} - X_{i,j} Y_{i,j}^2 + f (1 - X_{i,j})) \Delta t \\[.5em] Y_{i,j}^{\text{new}} = Y_{i,j} + (\nu \nabla^2 Y_{i,j} + X_{i,j} Y_{i,j}^2 - (k + f) Y_{i,j}) \Delta t \end{cases} \end{equation} dove

• $i, \, j \,$ sono le coordinate della cellula;
• $\nabla \,$ indica un prodotto di convoluzione bidimensionale, con pesi $\, - 1 \,$ per la cellula $\, (i, j) \,$, $\, \frac{1}{2} \,$ per quelle lateralmente adiacenti, $\, \frac{1}{4} \,$ per quelle diagonalmente adiacenti;
• $f \,$ è il cosiddetto fattore feed;
• $k \,$ è il fattore kill.

Fonte: karlsims.com/rd.html

Bibliografia

• A. M. Turing, The Chemical Basis of Morphogenesis , Philosophical Transactions of the Royal Society of London, 1952.
• P. Ball, Forging patterns and making waves from biology to geology: a commentary on Turing (1952) ‘The chemical basis of morphogenesis’ , Philosophical Transactions of the Royal Society of London, 2015.

Sitografia

• K. Sims, Reaction-Diffusion Tutorial, via karlsims.com.
• D. Shiffman, Coding Challenge #13: Reaction Diffusion Algorithm in p5.js, via You Tube.
• Ian, Turing Patterns, via Degenerate State.