Sebastian - Projekt

Collective Cell Invasion


[1] C.M. Material zum Projekt-Start


Erste Ideen und Simulationen:

Biased diffusion model
First experiments
Effect of cell proliferation
Numerical solution of drift-diffusion equation
Cooperation via Chemotaxis
First test of chemotactic model

Programme:

CellInvasion0.zip

Nach dem Entpacken müssen der Bibliotheken-Ordner "CMLib" und der Order mit dem Simulationsprogramm "invade3" nebeneinander im gleichen Über-Verzeichnis liegen. Wenn man invade3/invade3.dsw anklickt, sollte VisualC++6.0 aufgehen. Man kann das Programm dann direkt starten, bzw. editieren und neu kompilieren. Als Resultat werden diverse Konzentrations-Profile in den Ordner data geschrieben. Diese kann man dann z.B. bequem mit Easyplot (siehe unten) darstellen.

Easyplot

Easyplot ist einfach ein exe-File, es muß nicht installiert werden. Ein schlankes, schnelles Grafik-Programm, das aber alles wesentliche kann.


[2] C.M. Experimentelle Daten von Ben

Ben hat ein automatisches Auswertungsprogramm entwickelt. Auf der Grundlage vieler Mikroskop-Aufnahmen aus verschiedenen Tiefen-Ebenen des Gels können damit (kummulierte) Einwanderungs-Profile der Zellen generiert werden. Hier die ersten beiden Datensätze:

Hauszellen1
Darwinzellen1

Die Messung erfolgte t=3 Tage nach dem Start des Einwanderungs-Prozesses (in beiden Datensätzen ?)


[3] C.M. Fit der Daten an nicht-kollektive Theorie

Sebastian hat die Daten für die "Hauszellen1" und die "Darwinzellen1" mit der Lösung einer einfachen, nicht-kollektiven Fokker-Planck-Gleichung verglichen, welche lediglich Diffusion und ortsunabhängige Rückdrift annimmt. Im einfachsten Fall besitzt diese Theorie also nur 2 Fitparameter (Diffusionskonstante D und Rückdrift-Geschwindigkeit v). Der Fit liefert folgende Größenordnungen:

$D = 5\cdot 10^{-3} \mu m^2/s = 5\cdot 10^3 nm^2 / s$
$v = 0.1 nm/s$

Um das veränderte Einwanderungsverhalten unmittelbar unterhalb der Oberfläche zu berücksichtigen wurde weiterhin versucht, dort eine andere lokale Diffusionskonstante anzunehmen.

Die Resultate sind hier dokumentiert:

Resultate


[4] S.P. Neue Idee: Gemischte Population

Es gibt Grund zu der Annahme, dass es sich bei den Zellen nicht um eine Zellart mit fester Diffusionskonstante handelt, sondern eher um mehrere Populationen mit jeweils unterschiedlichen Diffusionskonstanten. Im Folgenden habe ich analytisch die Diffusionsgleichung in 1D gelöst mit einem festen Rand links. Man erhält als Verteilung eine Gaußfunktion. Diese wird dann noch integriert, für zwei verschiedene Diffusionskonstanten geplottet, und schließlich die Summe geplottet und mit logarithmischer y-Achse dargestellt:
$\frac{\partial c}{\partial t}=D \cdot \frac{\partial^2 c}{\partial x^2}$
Setzt man eine deltaförmige Anfangsverteilung $\delta(x) \cdot \gamma$ bei x=0 voraus folgt:
$\rightarrow c(x,t)=\frac{\gamma}{\sqrt{D \pi t}} e^{-\frac{x^2}{4 D t}}$ ($x \epsilon [0,\infty), t \epsilon [0,\infty)$
Wir integrieren:
$\int_x^\infty c(x',t) dx' = \int_x^\infty \frac{\gamma}{\sqrt{D \pi t}} e^{-\frac{x'^2}{4 D t}}dx'=\gamma(\text{erf}(0.5 \sqrt{\frac{1}{Dt}} \cdot x)-1)$
Die Verteilung zu einem bestimmten Zeitpunkt mit zwei unterschiedlichen Diffusionskonstanten sieht dann so aus:
ZweiZelltypen.PNG
Die gelb-braune Kurve ist die Summe der zwei Verteilungen und sieht den experimentellen Daten recht ähnlich.
Ein Fit an die Daten der Hauszellen ergibt:
population2.png
(höhere Auflösung)


[5] S.P.

Zwischenstand #1

experimentelle Erkenntnisse:

  1. Verhalten in Abhängigkeit von der Zelldichte
      • Zelldichte proportional zur Invasionstiefe/Invasionsgeschwindigkeit
      • Zelldichte proportional zum Anteil, der an der Oberfläche verweilenden Zellen (bzw. proportional zur Länge des Knicks am Anfang im log-Plot)
      • Länge des konstanten Teils im Bereich z=0 scheint unkorreliert mit Zelldichte (Messartefakt?)
      • => es existieren Zell-Zell-Wechselwirkungen (chemisch/mechanisch?)
  2. Zeitliches Verhalten
      • scheinbares Konvergieren gegen e-Funktion
      • langsames Verschwinden des Knicks am Anfang
  3. Sonstige Beobachtungen
      • fette Zellen (mit gefressenen Eisenkugeln) invadieren schlechter
      • Eindringen von Beads bis $70 \mu m$ in das Gel (Bead-Durchmesser: $4,5\mu m$)
      • Zellteilung ca. 1x pro tag
      • Sterberate ca. 0
      • Maschenweite des Kollagens ca. $1,3 \mu m$

Modell-Ansätze:

  1. reine Diffusion
      • einfachste Möglichkeit
      • erklärt Daten unzureichend
  2. viele Populationen
      • Modellannahme: Jede Population zeigt reines Diffusionsverhalten, aber es existieren 2 stark unterschiedliche Diffusionskonstanten. (siehe [4])
      • Stärke: erklärt starken Knick in der Verteilung bei kleinen z
      • Probleme: Abweichungen für große z
  3. Diffusion mit Rückdrift
      • Modellannahme: Die Zellen sitzen in einem näherungsweise Dreieckspotential (eig. exp), spüren also eine konstante, rücktreibende Kraft. (Ursache soll eine schnell diffundierende, exponentiell zerfallende Chemikalie sein, auf dessen Konzentrations-Gradienten sie reagiert.) (siehe [3] und First test of chemotactic model)
      • Stärke: erklärt exponentielles Verhalten im Mittelteil, konvergiert gegen exp-Verteilung
      • Probleme: Verhalten für kleine z nicht erklärt.
  4. Perkolation
      • Laufen Zellen nur entlang von Tunnel? Können sie Tunnel graben?
      • Wie ist die Tunnellängenverteilung? => Randeffekte
  5. weiter wilde Spekulationen
      • Einfluss der Schwerkraft (beim Bead Experiment), Tunnel?
      • Dynamik des Kollagens? (thermische Schlängelbewegung, ATP, Restpolymerisation?)

Fazit

Bisher zwei Theorien, die brauchbar erscheinen und einmal das Verhalten für kleine z (2 Diffusionskonstanten) und einmal für große z (Diffusion mit Rückdrift) erklären.

Ausblick

Als nächstes sollte ein Versuch gestartet werden, beide Theorien zu kombinieren:

  • die Zelle produziert eine schnelldiffundierende, mit der Zeit exp-zerfallende Chemikalie
  • sie selber reagiert auf den Gradienten und versucht so zu größeren Konzentrationen zu kommen (sprich zu den anderen Zellen)
  • Wenn eine bestimmte absolute Konzentration erreicht ist, dann schaltet die Zelle ihre Motoren ab, und somit sinkt ihre Diffusionskonstante auf ein Minimum ab. Nach dem Fit aus [4] sinkt die Diffusionskonstante um 2 Größenordnungen. (sprich, wenn in ihrer Umgebung genug Zellen sind, dann bleibt sie dort) (Das erklärt übrigens auch die experimentell gefundenen Abhängigkeiten im Bezug auf die Zelldichte.)

[6a] S.P.

Drift-Diffusions Modell mit Chemikalie und mit variabler Diffusionskonstante
Die Motivation, die oben in [5] "Ausblick" gegeben ist, soll hier konkret mathematisch formuliert werden (Es handelt sich um eine Modifikation der Gleichungen von Cooperation via Chemotaxis):

(I) Diffusion der Zellen:
$\frac{\partial}{\partial t}c(z,t) = - \frac{\partial}{\partial z}\left[v(z,t)\;c(z,t)\right] + \frac{\partial^2}{\partial z^2}[D_c(z,t) c(z,t)]$

(II) Driftgeschwindigkeit:
$v(z,t) = s\; \frac{\partial}{\partial z} g(z,t)$

(III) Sättigungsverhalten:
$D_c(z,t)=D_0+D_{c_0} [1-\Theta(g(z,t)-g_{krit})]$

(IV) Produktion und Zerfall der Locksubstanz:
$\frac{\partial}{\partial t}g(z,t) = p\!\cdot\!c(z,t) + D_g \frac{\partial^2}{\partial z^2}g(z,t) - d\!\cdot\!g(z,t)$

$D_0$ stellt die minimale Diffusion dar, die sich nicht vermeiden lässt. $D_{c_0}$ stellt die zusätzliche Diffusion dar, wenn die Zelle ihre Motoren einschaltet.
Sobald eine kritische Lockstoff Konzentration $g_{krit}$erreicht ist, dann greift die Theta Funktion zu und die Diffusionskonstante reduziert sich auf ein Minimum.
Für Bedeutung der restlichen Variablen, siehe Cooperation via Chemotaxis.
Bemerkung(1)
Gleichung (I) wurde umgeschrieben zu:
(I) Diffusion der Zellen:
$\frac{\partial}{\partial t}c(z,t) = - \frac{\partial}{\partial z}\left[v(z,t)\;c(z,t)\right] + D_c(z,t)\frac{\partial^2}{\partial z^2}[c(z,t)]$
Die Diffusionskonstante wurde aus der Differentiation herausgezogen, da sie wegen der Theta-Funktion fast überall konstant ist. Das Problem der unendlichen Steigung am Übergang entfällt dadurch.
Bemerkung(2a)
Im Bereich der kleinen Diffusionskonstante $D_c=D_0$ dominiert die Rückdrift. Der Knick am Anfang des log Plots der Verteilung ist dann sehr steil, denn wenn die Diffusionskonstante viel kleiner ist als die Rückdrift, vereinfacht sich die Gleichung zu
(I') $\frac{\partial}{\partial t}c(z,t) = - \frac{\partial}{\partial z}\left[v(z,t)\;c(z,t)\right] ]$
Dann wird im Laufe der Zeit der Knick immer steiler. Davon ist bisher nichts bekannt…
Bemerkung(2b)
Man kann nun folgendermaßen argumentieren:
Fall niedriger Konzentration: Die Zelle nimmt den Gradienten der Chemikalie war und verhält sich dementsprechend.
Fall kritischer Konzentration und höher: Die Sinne, die die Chemikaliengradienten feststellen können übersteuern. Die Zelle kann nun keinen Gradienten mehr messen und folglich gibt es keine Rückdrift mehr. Bei dieser Konzentration fühlt sich die Zelle wohl und sie Diffundiert mit einer minimalen Diffusion.
Nun schrumpft der Knick am Anfang nicht mehr, sonder fließt sehr (sehr) langsam auseinander.
Gleichung (II) modifiziert sich dann zu:
(II') $v(z,t) = s\; \frac{\partial}{\partial z} g(z,t) [1-\Theta(g(z,t)-g_{krit})]$
Anm: Fits ergeben, dass ein Gaußprofil sehr gut die Anfangsverteilung widerspiegelt.


[6b] S.P.

Diffusions Modell mit variabler Diffusionskonstante ohne Chemikalie
Das Umschalten der Diffusionskonstante, wie in [6a] Gl. (III) beschrieben, könnte auch durch die Zellkonzentration in unmittelbarer Nähe der Zelle ausgelöst werden.
Eine biologische Erklärung dafür wäre, dass die Zellen sich aneinander festhalten, sobald sie sich nahe gekommen sind. Innerhalb dieses Verbunds ist schwerer sich zu bewegen. Deswegen sinkt die Diffusionskonstante drastisch ab einer kritischen Konzentration, wo die Zellen zu nahe beieinander sind.
Dieses Modell beinhaltet keine Rückdrift.
Weiter unten wird der Prozess auch mit Monte Carlo Methoden nachgestellt.


[7] S.P.

wichtige Erkenntnis am Rande
Im log y Plot sieht ein verschobener Gauß einer Exp-Funktion sehr ähnlich. Also sind Modelle, die auf reiner Diffusion basieren immer noch wichtig!
$exp(-(x+c)^2)=exp(-x^2-2xc-c^2)$
Plot


[8] S.P.

(noch nicht fertig) Das zwei Populationen Modell
Ansatz: 2 Populationen diffundieren unabhängig voneinander.
$P(z)=A_1 exp(-z^2/\sigma_1)+A_2 exp(-z^2/\sigma_2)$
Im Folgenden Fits an experimentelle Daten. Welche Daten lassen sich dadurch beschreiben?


[9] S.P.

(noch nicht fertig) Das Doppel-Gauß-Modell
Ansatz: Die Zelle schaltet ab einem bestimmten Punkt, abhängig von der Konzentration, ihre Diffusionskonstante.
$z<z_krit: ~ P_1(z)=A_1 exp(-z^2/\sigma_1)$
$z>z_krit: ~ P_2(z)=A_2 exp(-(z+c)^2/\sigma_2)$
wobei P(z) immer stetig, vor allem bei P(z_krit)
Im Folgenden Fits an experimentelle Daten. Welche Daten passen dazu?


[10] S.P.

Nicht kumulierte Verteilungen aus Experimenten

Aus den Excel Tabellen von den Experimenten lassen sich die ursprünglichen Verteilungen zurückrechnen. Das stellt die experimentellen Daten in ein neues Licht.
Verteilung 50.000 Zellen
Verteilung (log y) 50.000 Zellen
Ganz besonders interessant ist der Bereich für kleine z: z=0 bis z=10
Hier sieht man, dass es zu einer Klumpenbildung kommt. Die einzelnen Zellhaufen lassen sich mit einem Gauß fitten (Daten ebenfalls 50.000 Zellen).
Gauß Fit
Diese Struktur sieht man in der kumulativen Verteilung gar nicht! Zum Vergleich hier die kumulieren Messdaten: kumuliert 50.000 z=0-10
kumuliert 50.000

Vergleichen wir mit anderen Datensätzen:

  • Auch hier sind diese Strukturen zu erkennen, teilweise auch sehr unregelmäßig. => statistisches Problem. Das tritt auch in MC-Simulation (siehe weiter unten) auf und liegt daran, dass nur um die 1000 Agenten an den realen Messungen teilgenommen haben. => Das Programm smooth.exe glättet die Daten.
  • Der Anfang der (geglätteten) Verteilung lässt sich sehr gut mit einem verschobenen Gauß fitten. => Eingentlich müsste die Verteilung aber ein halber Gauß sein… => random walks scheinen aber auch einen verschobenen Gauß zu produzieren (dem Problem muss ich nochmal auf den Grund gehen!=>liegt wohl an zu geringer Agentenzahl =>Rauschen)

[11] S.P.

Monte-Carlo Simulationen

  • Simulation von Agenten in einer 2D Box mit Rändern in x Richtung und perodischen Randbedingungen in y Richtung.
      1. Simulation mit einer Diffusionskonstante: Plot
      2. Simulation mit zwei Diffusionskonstanten, wobei der Wechsel bei einem festen x passiert: Plot
      3. Simulation mit zwei Diffusionskonstanten, wobei der Wechsel individuell bei jeder Zelle von der Anzahl der Nachbarn abhängt: Plot
  • hier ein Vergleich der obigen Simulationen mit einer Messung: Plot logPlot Es scheint, als ob das Umschalt-Modell abhängig von den Nachbarn am nächsten an die Daten herankommt. Eine neue Simulation mit leicht veränderten Diffusionskonstanten müsste dann noch besser passen…
  • Zum besseren Vergleich kann man mit dem Programm cumulate_data.exe die Daten aus der Simulation wieder kumulative Verteilungen umrechnen. Die Ähnlichkeit zu den experimentellen Daten ist verblüffend!
  • firstPcumPlot.PNG

[11] S.P.

Bemerkung zur Anfangsverteilung
Vermutung: In den Messdaten ist die wellige Oberfläche schon rausgerechnet. Was aber, wenn dieser Ausgleich zu grob ist und kleine Löcher, in die die Zellen aber reinpassen, übersieht? Angenommen die Löcher der ohnehin ausgefransten Collagenoberfläche sind bzgl. ihrer Tiefe gaußverteilt, dann ist auch die Anfangsverteilung der Zellen eine Gaußverteilung. (Und nicht wie ursprünglich angenommen eine Kastenverteilung oder ein halber Gauß.)
=> Wir brauchen eine Messung direkt nachdem die Zellen auf das Collagen aufgebracht wurden!


[12] S.P.

Ergebnisse der Monte-Carlo Simulationen
Verhalten in Abhängigkeit von der Dichte:
==> Simulation nach 3 Tagen:
Density_Dependance.PNG
==> Experiment nach 3 Tagen:
Density_Dependance_exp.PNG
Achtung: Diesen Plot müsste man eigentlich noch überarbeiten. Die absolute Zahl der Zellen gibt an, wie viele Zellen insgesamt verwendet wurden, nicht aber die Zellen pro Gesichtsfeld!!! In einem Gesichtsfeld befinden sich deutlich weniger Zellen (z.B.1000). Die genaue Zahl bekommt man aus dem Excel Sheet.


[13] S.P.

Analytische Lösung des Slow Layer Modells
Reduziert man den "Clustering Effect" auf das "Slow Layer Modell" und macht ein paar zusätzliche Annahmen, dann kann man das Problem auch analytisch Lösen und erhält diesen Plot (nichtkummulierte Verteilung!):
testplot2.gif


[14] S.P.

Analytische Lösung des Slow Layer Modells mehrdimensional
Berücksicht man die Inhomogenitäten senkrecht zur Invasionsrichtung erhält man diesen Plot (nichtkummulierte Verteilung!). Leider muss man das entstehende Integral numerisch lösen:
testplot3.gif


[15] S.P.

Das Drift-Diffusions Modell wurde unterschätzt!
Es ist sehr wohl in der Lage den anfänglichen steilen Abfall im Plot zu produzieren, wie er in den Experimente beobachtet wird:
ddm-fit.PNG

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License