BIO-LGCA

Dieser Artikel wurde auf der Qualitätssicherungsseite des Portals Mathematik eingetragen. Dies geschieht, um die Qualität der Artikel aus dem Themengebiet Mathematik auf ein akzeptables Niveau zu bringen. Bitte hilf mit, die Mängel dieses Artikels zu beseitigen, und beteilige dich bitte an der Diskussion! (Artikel eintragen)

Als BIO-LGCA (englisch biological lattice-gas cellular automata, biologische Gittergasautomaten) bezeichnet man in der mathematischen Biologie eine Art diskreter mathematischer Modelle zur Darstellung beweglicher, interagierender biologischer Agenten. Sie wurden aus den Gittergasautomaten (lattice-gas cellular automata, [LGCA]) der Fluiddynamik abgeleitet und stellen somit eine Unterart der zellulären Automaten (cellular automata, CA) dar. Ein BIO-LGCA-Modell beschreibt Zellen und andere frei bewegliche biologische Agenten als Punktpartikel, die sich auf einem räumlich diskreten Gitter fortbewegen und dabei mit anderen Partikeln in ihrer Umgebung interagieren. Im Gegensatz zu klassischen zellulären Automaten definieren BIO-LGCA die Partikel nicht nur durch ihre Position, sondern auch durch ihre Geschwindigkeit. Dadurch lassen sich Prozesse modellieren und analysieren, die hauptsächlich durch Änderungen des Impulses statt durch Änderungen der Dichte getrieben werden, wie aktive Fluide oder kollektive (Zell-)Migration. Anwendungen dieser Modellklasse finden sich z. B. in der Entschlüsselung der Dynamik von Krebsgeschwülsten.

Modelldefinition

Wie alle zellulären Automaten werden BIO-LGCA durch ein räumliches Gitter ${\mathcal {L}}$ (engl. lattice), einen Zustandsraum ${\mathcal {E}}$, eine Nachbarschaft ${\mathcal {N}}$ und eine Zustandsübergangsregel ${\mathcal {R}}$ definiert:

• Das Gitter (${\mathcal {L}}$) bestimmt die Menge aller möglichen (endlichen) Raumpositionen, die Partikel einnehmen können. Diese Punkte werden üblicherweise aus einer regulären und periodischen Pflasterung des Raums gewonnen. Mathematisch betrachtet ist ${\mathcal {L}}\subset \mathbb {R} ^{d}$ eine diskrete Teilmenge des $d$-dimensionalen Raums. Ein Gitterplatz $\mathbf {r} \in {\mathcal {L}}$ wird auch Knoten genannt.
• Der Zustandsraum (${\mathcal {E}}$) beschreibt die mögliche Konfiguration von Partikeln innerhalb jedes Gitterplatzes $\mathbf {r} \in {\mathcal {L}}$. Während sich in klassischen zellulären Automaten typischerweise nur ein Partikel gleichzeitig in einem Knoten befinden darf, erlauben BIO-LGCA mehrere Partikel mit unterschiedlichen Geschwindigkeiten an einem Gitterplatz. Dadurch ist der Zustandsraum etwas komplexer als in klassischen CA (siehe unten).
• Die Nachbarschaft (${\mathcal {N}}$) umfasst die Gitterplätze in der Umgebung eines Knotens $\mathbf {r} \in {\mathcal {L}}$, deren Zustand die Dynamik dieses Knotens bestimmt. Partikel können nur mit anderen Partikeln innerhalb dieser Nachbarschaft interagieren. Für die Nachbarschaft von Knoten am Rand endlicher Gitter müssen Randbedingungen definiert werden. Mögliche Nachbarschaften und Randbedingungen entsprechen denen für klassische zelluläre Automaten.
• Die Zustandsübergangsregel (${\mathcal {R}}$) legt das Verhalten der Partikel fest, d. h. wie sie sich über den Zeitverlauf fortbewegen, sich fortpflanzen oder sterben. Die Zeit in BIO-LGCA ist diskret. Um die Systemdynamik zu simulieren, wird die Regel in jedem Zeitschritt synchron auf alle Gitterplätze angewandt. Durch die Anwendung der Regel ändert sich der Originalzustand eines Gitterknotens zu einem neuen Zustand. Dieser hängt vom Zustand der Gitterplätze in der Nachbarschaft des betrachteten Knotens ab. Die Regel ist in BIO-LGCA in zwei Schritte aufgeteilt: probabilistische Interaktion und deterministischer Transport. Der Interaktionsschritt kann Neuausrichtung, Fortpflanzung und Tod beschreiben und wird für den modellierten biologischen Prozess spezifisch entworfen. Der Transportschritt befördert Partikel entsprechend ihrer Geschwindigkeit zu benachbarten Gitterplätzen. Dies wird unten näher ausgeführt.

Zustandsraum

Struktur eines BIO-LGCA-Gitterknotens mit sechs Geschwindigkeitskanälen (entsprechend den sechs Nachbarn in einem zweidimensionalen hexagonalen Gitter) und einem Ruhekanal. In diesem Fall sind $b=6$, $a=1$ und die Kapazität $K=7$. Die Kanäle 2, 3, 6 und 7 sind besetzt. Dadurch ergibt sich die Knotenkonfiguration zu $\mathbf {s} =(0,1,1,0,0,1,1)$ und die Anzahl der Partikel ist $n\left(\mathbf {s} \right)=\sum _{i=1}^{K}s_{i}=4$.

Um die Geschwindigkeit von Partikeln explizit zu modellieren, haben Gitterplätze eine besondere Struktur. Jeder Knoten $\mathbf {r} \in {\mathcal {L}}$ ist durch Vektoren mit den benachbarten Knoten verbunden, sog. „Geschwindigkeitskanäle“ $\mathbf {c} _{i}$, $i\in \{1,2,\ldots ,b\}$. Die Anzahl der Geschwindigkeitskanäle $b$ entspricht der Anzahl der benachbarten Gitterplätze und hängt daher von der Gittergeometrie ab ($b=2$ für ein eindimensionales Gitter, $b=6$ für ein zweidimensionales hexagonales Gitter etc.). Im zweidimensionalen Fall werden die Geschwindigkeitskanäle als $\mathbf {c} _{i}=\left(\cos {\frac {2\pi i}{b}},\sin {\frac {2\pi i}{b}}\right)$ definiert. Zusätzlich können beliebig viele Ruhekanäle $a$ festgelegt werden, die die Geschwindigkeit $\mathbf {c} _{i}=(0,0)$ haben, wobei $i\in \{b+1,b+2,\ldots ,b+a\}$. Ein Kanal gilt als besetzt, wenn die Geschwindigkeit eines Partikels im entsprechenden Gitterplatz der Geschwindigkeit dieses Kanals entspricht. Die Besetzung eines Kanals $\mathbf {c} _{i}$ wird durch die Besetzungszahl $s_{i}$ angezeigt. Typischerweise wahren die Partikel dabei das Ausschließungsprinzip: Nur ein Partikel kann sich zu einem gegebenen Zeitpunkt in einem Kanal $\mathbf {c} _{i}$ an einer Gitterposition befinden, $i\in \{1,2,\ldots ,b+a\}$. Ist das Prinzip gewahrt, sind die Besetzungszahlen Boolesche Variablen, d. h. $s_{i}\in {\mathcal {S}}=\{0,1\}$, und jeder Knoten hat eine maximale Kapazität $K=a+b$. Da die Besetzungszahlen aller Kanäle die Anzahl der im Knoten befindlichen Partikel und ihre Geschwindigkeiten erfassen, beschreibt der Vektor $\mathbf {s} =\left(s_{1},s_{2},\ldots ,s_{K}\right)$ vollständig den Zustand eines Gitterplatzes. Damit ergibt sich der Zustandsraum für einen Gitterplatz zu ${\mathcal {E}}={\mathcal {S}}^{K}$.

Regel und Dynamik des Modells

Die Zustände der Gitterplätze werden in jedem diskreten Zeitschritt entsprechend der Regel synchron aktualisiert, um die Dynamik des Modells zu simulieren. Die Zustandsübergangsregel besteht aus zwei Teilen: Im probabilistischen Interaktionsschritt treten die Partikel miteinander in Wechselwirkung, während der deterministische Transportschritt für ihre Fortbewegung verantwortlich ist.

Interaktionsschritt

Abhängig von der Anwendung kann der Interaktionsschritt aus Reaktions- und/oder Neuausrichtungsoperatoren bestehen.

Der Reaktionsoperator ${\mathcal {A}}$ überführt den Zustand $\mathbf {s}$ eines Gitterknotens gemäß einer Übergangswahrscheinlichkeit $P\left(\left.\mathbf {s} \rightarrow \mathbf {s} ^{\mathcal {A}}\right|\mathbf {s} _{\mathcal {N}}\right)$ in einen neuen Zustand $\mathbf {s} ^{\mathcal {A}}$. Die Wahrscheinlichkeit hängt vom Zustand $\mathbf {s} _{\mathcal {N}}$ der benachbarten Gitterplätze ab, um den Einfluss umgebender Partikel auf den Reaktionsprozess einzubeziehen. Mit dem Reaktionsoperator lassen sich Geburt und Tod von Agenten darstellen, da die Anzahl der Partikel im Knoten nicht erhalten werden muss. Die Formel für die Übergangswahrscheinlichkeit wird üblicherweise ad hoc aus phänomenologischen Betrachtungen heraus definiert.

Der Operator für die Neuausrichtung ${\mathcal {O}}$ aktualisiert einen Zustand $\mathbf {s}$ ebenfalls anhand einer Wahrscheinlichkeit $P\left(\left.\mathbf {s} \rightarrow \mathbf {s} ^{\mathcal {O}}\right|\mathbf {s} _{\mathcal {N}}\right)$ zu einem neuen Zustand $\mathbf {s} ^{\mathcal {O}}$. Hier besteht allerdings die Einschränkung, die Zahl der Partikel konstant zu halten. Daher modelliert dieser Operator durch die Umverteilung existierender Partikel in den Geschwindigkeitskanälen nur Änderungen in der Geschwindigkeit. Dies entspricht einer Richtungsänderung der einzelnen Agenten. Die Übergangswahrscheinlichkeit für den Operator kann aus statistischen Beobachtungen (durch Anwendung des Prinzips des maximalen Kalibers) oder aus der Dynamik einzelner Teilchen gewonnen werden. In letzterem Fall wird die diskretisierte stationäre Wahrscheinlichkeitsverteilung der Winkel genutzt, die durch die Fokker-Planck-Gleichung entsprechend der Langevin-Gleichung für die Dynamik der Neuorientierung gegeben ist. Die Übergangswahrscheinlichkeit nimmt typischerweise folgende Form an:

$$P\left(\left.\mathbf {s} \rightarrow \mathbf {s} ^{\mathcal {O}}\right|\mathbf {s} _{\mathcal {N}}\right)={\frac {1}{Z}}e^{-\beta H\left(\mathbf {s} _{\mathcal {N}}\right)}\delta _{n\left(\mathbf {s} \right),n\left(\mathbf {s} ^{\mathcal {O}}\right)}$$

$Z$ ist eine Normierungskonstante, in der statistischen Physik auch kanonische Zustandssumme. $\beta$ ist ein freier Parameter, der analog zur inversen Temperatur in der statistischen Thermodynamik umgekehrt proportional zur Zufälligkeit der Partikelneuorientierung ist. Die Funktion $H\left(\mathbf {s} _{\mathcal {N}}\right)$ nimmt eine Rolle ähnlich einer Energie ein: Durch ihre Neuausrichtung werden Partikel $H$ tendenziell minimieren. Das Kronecker-Delta $\delta _{n\left(\mathbf {s} \right),n\left(\mathbf {s} ^{\mathcal {O}}\right)}$ schließt alle Zustandsübergänge aus, bei denen die Anzahl der Partikel $n$ vor und nach der Anwendung des Operators ($n\left(\mathbf {s} \right)$, $n\left(\mathbf {s} ^{\mathcal {O}}\right)$) nicht gleich ist.

Der Zustand eines Gitterknotens nach Anwendung des Reaktions- und des Neuausrichtungsoperators $\mathbf {s} ^{{\mathcal {A}}\circ {\mathcal {O}}}$ wird Postinteraktions-Konfiguration genannt und als $\mathbf {s} ^{\mathcal {I}}:=\mathbf {s} ^{{\mathcal {A}}\circ {\mathcal {O}}}$ notiert.

Dynamik eines BIO-LGCA-Modells. Jeder Zeitschritt besteht aus den zwei Prozessen Interaktion (interaction) und Transport (migration), die an allen Gitterplätzen gleichzeitig erfolgen. Im Interaktionsschritt werden die Besetzungszahlen, d. h. die Zustände der Gitterplätze, durch den Reaktions- und/oder Neuausrichtungsoperator stochastisch aktualisiert. Im Beispiel wird eine partikelkonservierende Regel angenommen (nur Richtungsänderung, kein Reaktionsoperator). Anschließend, im Transportschritt, werden die Partikel deterministisch zum Geschwindigkeitskanal der gleichen Richtung im entsprechenden benachbarten Knoten bewegt. Die Farben in der Skizze zeigen, in welchem Knoten sich die Partikel am Anfang des Zeitschritts befunden haben, um den Transport zu verdeutlichen.

Transportschritt

Nach dem Interaktionsschritt wird der deterministische Transportschritt synchron auf alle Gitterplätze angewandt. Er simuliert die Bewegung von Agenten entsprechend ihrer Geschwindigkeit und repräsentiert damit den Eigenantrieb lebender Organismen oder Zellen. Formal bestimmt die Besetzungskonfiguration der Geschwindigkeitskanäle nach dem Interaktionsschritt die neue Besetzung der Geschwindigkeitskanäle mit der gleichen Richtung in den benachbarten Gitterplätzen, d. h. $s_{i}(\mathbf {r} +\mathbf {c} _{i})=s_{i}^{\mathcal {I}}(\mathbf {r} )$, vgl. auch die Abbildung.

Wenn Interaktion und Transport abgeschlossen sind, beginnt ein neuer Zeitschritt. Zusammengefasst kann die Dynamik des BIO-LGCA als stochastische mikrodynamische Rekursionsbeziehung formuliert werden:

$$s_{i}(\mathbf {r} +\mathbf {c} _{i},k+1)=s_{i}^{\mathcal {I}}(\mathbf {r} ,k)$$

Beispiele für Interaktionen

BIO-LGCA-Modell polarer Schwarmbildung. In diesem Modell richten Zellen ihre Geschwindigkeit vorzugsweise parallel zum Impuls in ihrer Nachbarschaft aus. Die Farbe der Gitterknoten wird durch die Orientierung der enthaltenen Partikel und deren Übersetzung in den Farbkreis bestimmt. Leere Knoten sind weiß. Die Simulation wurde auf einem hexagonalen Gitter mit periodischen Randbedingungen durchgeführt.

Die Übergangswahrscheinlichkeit für den Reaktions- und/oder Neuausrichtungsoperator ist der Schlüssel, um das Verhalten eines modellierten Systems angemessen abzubilden. Einige elementare Interaktionen und ihre entsprechenden Übergangswahrscheinlichkeiten sind hier aufgeführt.

Random Walk

Wenn keine externen oder internen Stimuli vorliegen, können sich Zellen zufällig ohne eine bevorzugte Richtung bewegen (Random Walk). In diesem Fall kann der Neuausrichtungsoperator mit der Übergangswahrscheinlichkeit

$$P\left(\left.\mathbf {s} \rightarrow \mathbf {s} ^{\mathcal {O}}\right|\mathbf {s} _{\mathcal {N}}\right)={\frac {\delta _{n(\mathbf {s} ),n\left(\mathbf {s} ^{\mathcal {O}}\right)}}{Z}}$$

realisiert werden, wobei $Z=\sum _{\mathbf {s} ^{\mathcal {O}}}\delta _{n\left(\mathbf {s} \right),n\left(\mathbf {s} ^{\mathcal {O}}\right)}$. Diese Übergangsregel erlaubt als Ergebnis der Richtungsänderung mit gleicher Wahrscheinlichkeit jede Konfiguration $\mathbf {s} ^{\mathcal {O}}$, die die gleiche Anzahl von Partikeln aufweist wie die Konfiguration $\mathbf {s}$ vor der Anwendung des Operators.

BIO-LGCA-Modell für erregbare Medien. In diesem Modell fördert der Reaktionsoperator eine schnelle Reproduktion von Partikeln in Geschwindigkeitskanälen, während sie in Ruhekanälen langsam sterben. Partikel in Ruhekanälen hemmen dabei die Fortpflanzung in den Geschwindigkeitskanälen. Die Neuausrichtung erfolgt nach dem Random Walk-Operator im Text. Gitterzellen sind nach der Anzahl der beweglichen Partikel farbcodiert: Je mehr Partikel sich in Geschwindigkeitskanälen aufhalten, desto heller ist die Stelle gefärbt. Ruhende Partikel werden nicht gezeigt. Die Simulation wurde auf einem hexagonalen Gitter mit periodischen Randbedingungen durchgeführt.

Einfacher Geburts- und Todesprozess

Wenn Organismen sich unabhängig von anderen Individuen fortpflanzen und sterben (mit Ausnahme der begrenzten Kapazität ihrer Umgebung), kann mit der Übergangswahrscheinlichkeit

$$P\left(\left.\mathbf {s} \rightarrow \mathbf {s} ^{\mathcal {A}}\right|\mathbf {s} _{\mathcal {N}}\right)=\left[r_{b}\delta _{n\left(\mathbf {s} ^{\mathcal {A}}\right),n\left(\mathbf {s} \right)+1}+r_{d}\delta _{n\left(\mathbf {s} ^{\mathcal {A}}\right),n\left(\mathbf {s} \right)-1}\right]\Theta \left[n\left(\mathbf {s} ^{\mathcal {A}}\right)\right]\Theta \left[n\left(K-\mathbf {s} ^{\mathcal {A}}\right)\right]$$

ein einfacher Geburts- und Todesprozess modelliert werden. $r_{b},r_{d}\in [0,1]$, $r_{b}+r_{d}\leq 1$ sind konstante Geburts- und Sterbewahrscheinlichkeiten. Das Kronecker-Delta $\delta _{i,j}$ stellt sicher, dass in jedem Zeitschritt nur ein Geburts- bzw. Todesprozess stattfindet, während die Heaviside-Funktion $\Theta (x)$ die Partikelanzahl auf positive Werte kleiner oder gleich der Kapazität $K$ begrenzt.

Adhäsive Wechselwirkungen

BIO-LGCA-Modell für Adhäsion. Zellen bewegen sich bevorzugt in Richtung des Zelldichtegradienten. Gitterplätze im Video sind mit zunehmender Zelldichte in intensiverem Blau dargestellt. Leere Knoten sind weiß. Die Simulation wurde auf einem quadratischen Gitter mit periodischen Randbedingungen durchgeführt.

Zellen können durch Cadherin-Moleküle an der Zelloberfläche aneinander haften. Dieser Mechanismus erlaubt ihnen die Bildung von Aggregaten. Die Herausbildung von Zellaggregaten durch adhäsive Biomoleküle lässt sich durch einen Neuausrichtungsoperator mit den Übergangswahrscheinlichkeiten

$$P\left(\left.\mathbf {s} \rightarrow \mathbf {s} ^{\mathcal {O}}\right|\mathbf {s} _{\mathcal {N}}\right)={\frac {1}{Z}}\exp \left[\beta \mathbf {G} \left(\mathbf {s} _{\mathcal {N}}\right)\cdot \mathbf {J} \left(\mathbf {s} ^{\mathcal {O}}\right)\right]$$

modellieren, wobei der Vektor $\mathbf {G} \left(\mathbf {s} _{\mathcal {N}}\right)$ in die Richtung der höchsten Zelldichte zeigt. Er ist definiert als $\mathbf {G} \left(\mathbf {s} _{\mathcal {N}}\right)=\sum _{\mathbf {r} '\in {\mathcal {N}}}\left(\mathbf {r} '-\mathbf {r} \right)n\left(\mathbf {s} _{\mathcal {N}}^{\mathbf {r} '}\right)$. $\mathbf {s} _{\mathcal {N}}^{\mathbf {r} '}$ bezeichnet die Konfiguration des Gitterplatzes $\mathbf {r} '$ innerhalb der Nachbarschaft ${\mathcal {N}}$ eines Knotens $\mathbf {r}$ (vor dem Interaktionsschritt). $\mathbf {J} \left(\mathbf {s} ^{\mathcal {O}}\right)=\sum _{j=1}^{b}s_{j}^{\mathcal {O}}\mathbf {c} _{j}$ ist der Impuls des Knotens nach dem Interaktionsschritt. Diese Übergangsregel begünstigt die Wahl von Konfigurationen, bei denen sich die Zellen entlang des Gradienten der Zelldichte bewegen.

BIO-LGCA-Modell von Zellen, die indirekt durch Chemotaxis interagieren. Individuen produzieren einen diffundierenden Lockstoff mit einer bestimmten biologischen Halbwertszeit und bewegen sich bevorzugt in Richtung des Gradienten dieses Stoffes. Die Farbe eines Knotens im Video ergibt sich durch eine zunehmend dunklere Blaufärbung nach der Zelldichte und eine zunehmend intensivere Gelbfärbung nach der Konzentration des Lockstoffes. Leere Gitterplätze sind weiß. Die Simulation wurde auf einem quadratischen Gitter mit periodischen Randbedingungen durchgeführt.

Mathematische Analyse

Eine exakte Betrachtung stochastischer agentenbasierter Modelle ist aufgrund von Korrelationen höherer Ordnung zwischen allen Agenten oft schwierig. Eine approximative Analysemethode für BIO-LGCA-Modelle besteht darin, sie in eine geeignete deterministische Differenzengleichung für die mittlere Populationsdynamik zu überführen, diese Gleichung mathematisch zu analysieren und aus den Ergebnissen Rückschlüsse auf das Verhalten des ursprünglichen BIO-LGCA-Modells zu ziehen.

Zuerst wird der Erwartungswert der mikrodynamischen Gleichung $s_{m}(\mathbf {r} +\mathbf {c} _{i},k+1)=s_{m}^{\mathcal {I}}(\mathbf {r} ,k)$ bestimmt:

$$f_{m}\left(\mathbf {r} +\mathbf {c} _{m},k+1\right)=\left\langle s_{m}^{\mathcal {I}}\left(\mathbf {r} ,k\right)\right\rangle$$

$\langle \cdot \rangle$ steht für den Erwartungswert. $f_{m}\left(\mathbf {r} ,k\right):=\left\langle s_{m}\left(\mathbf {r} ,k\right)\right\rangle$ bezeichnet den Erwartungswert für die Besetzungszahl des $m$-ten Kanals $\mathbf {c} _{m}$ am Gitterplatz $\mathbf {r}$ zum Zeitpunkt $k$. Der Term auf der rechten Seite der Gleichung, $\left\langle s_{m}^{\mathcal {I}}\left(\mathbf {r} ,k\right)\right\rangle$, ist hochgradig nichtlinear im Bezug auf die Besetzungszahlen des Knotens $\mathbf {r}$ und der Knoten in der Interaktionsnachbarschaft ${\mathcal {N}}$. Dies ergibt sich aus der Form der Übergangswahrscheinlichkeit $P\left(\left.\mathbf {s} \rightarrow \mathbf {s} ^{\mathcal {I}}\right|\mathbf {s} _{\mathcal {N}}\right)$ und den statistischen Eigenschaften der Verteilung von Partikeln auf die Geschwindigkeitskanäle, z. B. durch die Forderung, dass nur ein Partikel einen Kanal besetzen darf (Ausschließungsprinzip, s. o.). Die Nichtlinearität der rechten Seite würde zu den oben erwähnten Korrelationen und Momenten höherer Ordnung zwischen allen beteiligten Kanalbesetzungszahlen führen. Daher wird meist eine Molekularfeldnäherung genutzt, die alle Korrelationen und Momente höherer Ordnung vernachlässigt. Direkte Partikel-Partikel-Interaktionen werden durch Produkte ihrer Erwartungswerte approximiert, anders ausgedrückt: Es seien $X_{1},X_{2},\ldots ,X_{n}$ Zufallsvariablen und $F:\mathbb {R} ^{n}\mapsto \mathbb {R}$ eine Funktion. Dann gilt unter der Molekularfeldnäherung: $\left\langle F\left(X_{1},X_{2},\ldots ,X_{n}\right)\right\rangle \approx F\left(\left\langle X_{1}\right\rangle ,\left\langle X_{2}\right\rangle ,\ldots ,\left\langle X_{n}\right\rangle \right)$. Dadurch kann die Gleichung für den Erwartungswert vereinfacht werden zu:

$$f_{m}\left(\mathbf {r} +\mathbf {c} _{m},k+1\right)={\mathcal {C}}\left(\mathbf {f} \left(\mathbf {r} ,k\right),\mathbf {f} _{\mathcal {N}}\left(\mathbf {r} ,k\right)\right)$$

Dabei ist ${\mathcal {C}}\left(\mathbf {f} \left(\mathbf {r} ,k\right),\mathbf {f} _{\mathcal {N}}\left(\mathbf {r} ,k\right)\right)$ eine nichtlineare Funktion der erwarteten Knotenkonfiguration $\mathbf {f} \left(\mathbf {r} ,k\right)$ und der erwarteten Konfiguration der Interaktionsnachbarschaft $\mathbf {f} _{\mathcal {N}}\left(\mathbf {r} ,k\right)$. Sie hängt von den Übergangswahrscheinlichkeiten und der statistischen Verteilung der Partikel innerhalb der Knoten ab.

Für diese nichtlineare Differenzengleichung lassen sich typischerweise mehrere homogene stationäre Zustände oder Konstanten ${\bar {f}}_{m}$ finden, welche die Differenzengleichung lösen und unabhängig von $\mathbf {r}$ und $k$ sind. Zur Untersuchung der Bedingungen für ihre Stabilität und des Potenzials eines Modells zur Musterbildung kann eine lineare Stabilitätsanalyse vorgenommen werden. Dazu wird die nichtlineare Rekursionsbeziehung linearisiert:

$$f_{m}\left(\mathbf {r} +\mathbf {c} _{m},k+1\right)=\sum _{j=1}^{K}\left.{\frac {\partial {\mathcal {C}}}{\partial f_{j}\left(\mathbf {r} ,k\right)}}\right|_{\mathrm {ss} }f_{j}\left(\mathbf {r} ,k\right)+\sum _{j=1}^{K}\sum _{p=1}^{K}\left.{\frac {\partial {\mathcal {C}}}{\partial f_{j}\left(\mathbf {r} +\mathbf {c} _{p},k\right)}}\right|_{\mathrm {ss} }f_{j}\left(\mathbf {r} +\mathbf {c} _{p},k\right)$$

In der Herleitung wurde eine Von-Neumann-Nachbarschaft angenommen. Der Index $\mathrm {ss}$ steht für den homogenen stationären Zustand (engl. steady state) $f_{m}\left(\mathbf {r} ,k\right)={\bar {f}}_{m},m\in \{1,\ldots ,K\}$. Um eine Gleichung mit ausschließlich zeitlichen Inkrementen zu erhalten, kann anschließend für den Raum eine diskrete Fourier-Transformation (DFT) auf beide Seiten angewandt werden. Durch die Nutzung der Verschiebungseigenschaft der DFT lässt sich die räumliche Verschiebung von vormals $\mathbf {r} +\mathbf {c} _{p}$ in eine Multiplikation mit dem Term $e^{{\frac {2\pi i}{L}}\mathbf {q} \cdot \mathbf {c} _{p}}$ umformen. Nach der Isolation des zeitlichen Inkrements auf der linken Seite ergibt sich die Lattice-Boltzmann-Gleichung:

$${\hat {f}}_{m}\left(\mathbf {q} ,k+1\right)=e^{-{\frac {2\pi i}{L}}\mathbf {q} \cdot \mathbf {c} _{m}}\left\{\sum _{j=1}^{K}\left[\left.{\frac {\partial {\mathcal {C}}}{\partial f_{j}\left(\mathbf {r} ,k\right)}}\right|_{\mathrm {ss} }+\sum _{p=1}^{K}\left.{\frac {\partial {\mathcal {C}}}{\partial f_{j}\left(\mathbf {r} +\mathbf {c} _{p},k\right)}}\right|_{\mathrm {ss} }e^{{\frac {2\pi i}{L}}\mathbf {q} \cdot \mathbf {c} _{p}}\right]{\hat {f}}_{j}\left(\mathbf {q} ,k\right)\right\}$$

$i={\sqrt {-1}}$ bezeichnet die imaginäre Zahl, $L$ ist die Größe des Gitters entlang einer Dimension, $\mathbf {q} \in \{1,2,\ldots ,L\}^{d}$ steht für die Fourier-Wellenzahl und ${\hat {\cdot }}={\mathcal {F}}\{\cdot \}$ zeigt die diskrete Fourier-Transformation an. In Matrixschreibweise vereinfacht sich die Gleichung zu ${\hat {\mathbf {f} }}\left(\mathbf {q} ,k+1\right)=\Gamma {\hat {\mathbf {f} }}\left(\mathbf {q} ,k\right)$, wobei die Matrix $\Gamma$ Boltzmann-Propagator genannt wird und als

$$\Gamma _{m,j}=e^{-{\frac {2\pi i}{L}}\mathbf {q} \cdot \mathbf {c} _{m}}\left[\left.{\frac {\partial {\mathcal {C}}}{\partial f_{j}\left(\mathbf {r} ,k\right)}}\right|_{\mathrm {ss} }+\sum _{p=1}^{K}\left.{\frac {\partial {\mathcal {C}}}{\partial f_{j}\left(\mathbf {r} +\mathbf {c} _{p},k\right)}}\right|_{\mathrm {ss} }e^{{\frac {2\pi i}{L}}\mathbf {q} \cdot \mathbf {c} _{p}}\right]$$

definiert ist. Die Eigenwerte $\lambda \left(\mathbf {q} \right)$ des Boltzmann-Propagators bestimmen die Stabilitätseigenschaften des stationären Zustands:

• Wenn $\left|\lambda \left(\mathbf {q} \right)\right|>1$ ist ($|\cdot |$ bezeichnet den Betrag), wachsen Störungen mit der Wellenzahl $\mathbf {q}$ mit der Zeit an. Wenn $\left|\lambda \left(\mathbf {q} _{\mathrm {max} }\right)\right|>1$ und $\left|\lambda \left(\mathbf {q} _{\mathrm {max} }\right)\right|\geq \left|\lambda \left(\mathbf {q} \right)\right|\forall \mathbf {q} \in \{1,2,\ldots ,L\}^{d}$, dominieren Perturbationen mit der Wellenzahl $\mathbf {q} _{\mathrm {max} }$ und es treten räumliche Muster mit einer deutlich erkennbaren Wellenlänge auf. In allen anderen Fällen ist der stationäre Zustand stabil, Störungen klingen ab.
• Wenn $\mathrm {arg} \left[\lambda \left(q\right)\right]\neq 0$ ist ($\mathrm {arg} (\cdot )$ steht für das Argument), pflanzen sich Störungen in der Domäne fort und führen zu nichtstationärem Populationsverhalten. Anderenfalls scheint die Population auf einer makroskopischen Ebene betrachtet stillzustehen.

Anwendungen

Das Hauptaugenmerk beim Entwurf eines BIO-LGCA-Modells zur Untersuchung eines biologischen Phänomens liegt auf der Definition geeigneter Übergangswahrscheinlichkeiten für den Interaktionsschritt. Daneben können für spezielle Anwendungsfälle auch weitere Anpassungen von Bedeutung sein:

• Durch Änderungen des Zustandsraums können z. B. verschiedene Phänotypen von Zellen eingebunden werden.
• Für die Modellierung von Vorgängen in räumlich begrenzten Bereichen ist die Wahl der Randbedingungen entscheidend.
• Um z. B. experimentelle Interaktionsreichweiten quantitativ nachzubilden, wird die Interaktionsnachbarschaft entsprechend verändert.
• Die Kapazität eines Knotens kann in Verbindung mit der Skalierung des BIO-LGCA zum realen Raum beispielsweise genutzt werden, um Crowding-Effekte für eine bestimmte Zellgröße darzustellen.

Während sich die Wahrscheinlichkeitsverteilung für den Neuausrichtungsoperator mit den oben genannten statistischen und biophysischen Methoden herleiten lässt, kann die Übergangswahrscheinlichkeit für den Reaktionsoperator auch durch die statistische Auswertung von Experimenten in vitro gewonnen werden.

In der Praxis werden BIO-LGCA-Modelle genutzt, um zelluläre, biophysikalische und medizinische Phänomene zu untersuchen. Einige Beispiele umfassen:

• Angiogenese: Dieser Aufsatz vergleicht Beobachtungswerte aus einem in vitro-Experiment mit Endothelzellen und aus BIO-LGCA-Simulationen, um die an der Blutgefäßbildung beteiligten Prozesse und ihre jeweilige Bedeutung zu ermitteln. Im Ergebnis werden die Adhäsion zwischen Zellen, Alignment, Kontaktführung und ein Umbau der extrazellulären Matrix (EZM) als relevant identifiziert. Wechselwirkungen mit größerer Reichweite sind nicht entscheidend für die Angiogenese.
• Aktive Fluide: Mithilfe eines BIO-LGCA-Modells werden die makroskopischen physikalischen Eigenschaften einer Population von Partikeln untersucht, die über polares Alignment miteinander interagieren. Eine Erhöhung der anfänglichen Partikeldichte oder der Interaktionsstärke führt zu einem Phasenübergang zweiter Ordnung: Ein homogener, ungeordneter Anfangszustand geht in einen geordneten, strukturierten, aber bewegten Zustand über.
• Epidemiologie: In diesem Aufsatz wird zur Analyse verschiedener Impfstrategien ein räumliches SIR-Modell im BIO-LGCA-Framework implementiert und zusätzlich genutzt, um den Effekt der Näherung einer räumlich aufgelösten Epidemie mit einem nicht-räumlichen Modell zu untersuchen. Im Modell zeigen sich Impfstrategien in Form von Barrieren wesentlich effektiver als überall einheitlich vorangetriebene Impfungen. Im Vergleich mit räumlichen Modellen überschätzen nicht-räumliche Modelle die Infektionsraten enorm.
• Jamming von Zellen: In vitro- und BIO-LGCA-Modelle werden hier für die Studie der Ausbreitung und Metastasierung von Brustkrebs genutzt. Das Simulationsmodell zeigt auf, dass Metastasen abhängig von der Stärke der Zelladhäsion, der EZM-Dichte und den Interaktionen zwischen Zellen und EZM verschiedene Verhaltensweisen annehmen können. Dazu gehören eine zufällige gasähnliche, eine korrelierte flüssigkeits-ähnliche und eine blockierte Festzustands-ähnliche Beschaffenheit.

Weblinks

• Bio-LGCA-Simulator – Ein Online-Tool für die Simulation von BIO-LGCA im Browser mit elementaren Interaktionen und einstellbaren Parameterwerten.
• BIO-LGCA Python Package – Ein Open-Source-Softwarepaket zur Implementation von BIO-LGCA-Modellsimulationen. Es enthält einige elementare Interaktionen und ist für eigene Regeln erweiterbar.