Anja - Projekt

[1] C.M.: Projekt-Idee

Selbstantriebene Agenten in zufälligen Potential-Landschaften

Das Problem der Diffusion von Teilchen in Potential-Landschaften taucht in vielen Bereichen der Physik auf. Normalerweise werden die Teilchen von thermischen Zufallskräften (weißes Rauschen) getrieben. Der Transport über längere Distanzen erfordert die Überwindung von Barrieren zwischen benachbarten Tälern der Potential-Landschaft. Die Diffusion wird limitiert durch die "mean first passage time" für solche Barrieren-Überschreitungen, welche exponentiell mit der Barrierenhöhe zunimmt.

Vor einem ähnlichen Transport-Problem stehen lebende Zellen, die durch biologisches Gewebe wandern. Handelt es sich bei diesem Gewebe z.B. um ein engmaschiges, elastisches Faser-Netzwerk (Kollagen..), so erfordert die Migration der Zelle auch hier das wiederholte Überwinden von Engstellen. In einfachster Näherung kann wieder jedem Zellort (= Mittelpunkt der kugelförmigen Modell-Zelle) eine elastische potentielle Energie zugeordnet werden. Im Gegensatz zu einem thermisch getriebenen Teilchen in der Potential-Landschaft generiert die Zelle jedoch aktive Kräfte, es handelt sich also um einen selbstgetriebenen Agenten. Die Zelle kann ihren Krafterzeugungs-Apparat ("Aus-Strecken" von Fillopodien, mechanische Verbindung mit der Umwelt durch fokale Adhäsionen, Bildung von Akto-Myosin-Stressfasern, Erzeugung von Zugkräften, Auflösen von Adhäsionen ..) sicherlich zu einem gewissen Grad steuern. Allerdings scheint bei der zeitlichen Dynamik der resultierenden Antriebs-Kräfte auch ein erheblicher, rein stochastischer Anteil vorzuliegen.

Es stellt sich nun die Frage, mit welcher Krafterzeugungs-Strategie die Zelle in einer gegebenen Potential-Landschaft am effektivsten und effizientesten vorankommt.

Im einfachsten Fall ist die Zellkraft als Funktion der Zeit ein reiner Zufalls-Prozeß, völlig unabhängig von der lokalen Eigenschaften der externen Potential-Landschaft ("Blinde Strategien"). Der zeitlich fluktuierende Kraftvektor kann dann z.B. durch Wahrscheinlichkeits-Verteilungen $P(\vec{F})$ und Autokorrelations-Funktionen $C_{\vec{F}\vec{F}}(\tau)$ beschrieben werden. Angenommen, die von der Zelle investierte mittlere Gesamt-Leistung (Integral über die frequenzabhängige Leistungsdichte der Zellkraft) sowie die Potential-Landschaft seien vorgegeben. Welche statistischen Eigenschaften $P(\vec{F})$ und $C_{\vec{F}\vec{F}}(\tau)$ sind dann optimal im Sinne einer möglichst "schnellen" Diffusion ?

Noch interessanter sind "Strategien mit Rückkopplung": Angenommen, die Zelle kann nach jedem Zeitschritt $n$ feststellen, welche Verschiebung $\Delta\vec{R}_n$ die Anwendung der Kraft $\vec{F}_n$ bewirkt hat und auf der Grundlage dieser lokalen Umwelt-Information die Wahrscheinlichkeiten für den folgenden Kraftschritt $\vec{F}_{n+1}$ beeinflussen (Sensor-Motor-Schleife). Welche Strategie ist dann optimal bei vorgegebenen Randbedingungen ?


[2] C.M.: Alpha-Modell

Wir beginnen mit einem sehr einfachen Modell: Ein punktförmiger Agent befindet sich in einer eindimensionalen Potential-Landschaft $U(x)$ und spürt aufgrund einer Hintergrund-Viskosität geschwindigkeits-proportionale Reibungskräfte $F_{frict}=-\;\gamma\;v$ mit Reibungskonstante $\gamma$. Der Agent erzeugt zeitabhängige Antriebskräfte $F(t)$ mit der Verteilungsfunktion $P(F)$ und der Autokorrelation $C_{FF}(\tau)$. Seine Bewegungsgleichung lautet somit in überdämpfter Näherung:

(1)
\begin{align} \dot{x}=\left( F(t) - \frac{\partial}{\partial x}U(x) \right)/\gamma. \end{align}

Dies hat die Struktur eine Langevin-Gleichung $\dot{x}=D(x)+S(t)$ mit deterministischem Term $D(x)=-\frac{\partial}{\partial x}U(x)/\gamma$ und stochastischem Term $S(t)=F(t)/\gamma$. Die statistischen Eigenschaften der Zufallskraft $F(t)$ werden in unserem Fall aber komplexer sein als bei den üblichen Langevin-Gleichungen für thermische Teilchen.

Insbesondere betrachten wir zwei Arten von prototypischen Kraft-Fluktuationen:

  • Exponential-korrelierte, gaußverteilte (EG-) Fluktuationen:
(2)
\begin{align} C_{FF}(\tau) = \sigma^2 e^{-\tau/\tau_c}\;\;\mbox{und}\;\;P(F)\propto e^{-\frac{1}{2}(F/\sigma)^2} \end{align}
  • Langzeit-korrelierte, "random-phase" (LR-) Fluktuationen:
(3)
\begin{align} C_{FF}(\tau) \propto \tau^{-\beta}\;\;\mbox{falls}\;\;\tau\in [\tau_1,\tau_2]. \end{align}

Durch die Beschränkung der Zeitkonstanten auf ein breites, aber endliches Intervall bleibt das Leistungsspektrum $|F(\omega)| ^2$ normierbar.

Die Bezeichung "random-phase" bedeutet hier, daß die Phasen $\varphi(\omega)$ der Kraft-Fourier-Komponenten $\tilde{F}(\omega)=|F(\omega)|\cdot e^{i\varphi(\omega)}$ voneinander statistisch unabhängig und im Intervall $[0,2\pi]$ gleichverteilt gewählt werden.

  • Konstante Kraft: $F(t)=F_0$.
  • Periodische Kraft: $F(t)=f_0\;\sin(\om t)$.

Folgende Potential-Landschaften könnten untersucht werden:

  • Flaches Potential $U(x)=0$
  • Harmonisches Potential $U(x)=\frac{k}{2}x^2$
  • Gaußpotential $U(x)=-h\; e^{-\frac{1}{2}(x/w)^2}$
  • Doppelgaußpotential $U(x)=-h\; \left( e^{-\frac{1}{2}((x+\frac{d}{2})/w)^2} + e^{-\frac{1}{2}((x-\frac{d}{2})/w)^2} \right)$
  • Multigaußpotential $U(x)= \sum_m -h_m\; e^{-\frac{1}{2}((x-x_m)/w_m)^2}$
  • Fourier-Reihen-Potential $U(x)=\sum_m a_m \;\sin(k_mx-\eta_m)$

[2] C.M.: Alpha-Modell: Numerische Simulation

  • EG-Kraftfluktuationen werden mit der Box-Müller-Methode erzeugt.
  • Für die LR-Fluktuationen wird die Fourier-Methode verwendet.
  • Wir machen zusätzlich die Näherung, daß die Zufallskraft $F(t)$ stückweise konstant ist über kleine Simulations-Zeitintervalle $\Delta t_{sim}$. Dieses Zeitintervall muß immer kleiner gewählt werden als die kürzeste Korrelations-Zeitkonstante der Kraft.
  • Die Langevin-Gleichung wird mit einem Runge-Kutta-Verfahren integriert.
  • Die resultierende Teilchen-Trajektorie $x(t)$ wird dann statistisch ausgewertet. Insbesondere wird das Mean Squared Displacement (MSD) bestimmt: $\overline{x^2}(\tau)=<(x(t+\tau)-x(t))^2>_t$.

[3] C.M.: Alpha-Modell: Erste Tests

Zunächst sollen zum Test der numerischen Verfahren einige einfache Fälle simuliert werden, für die analytische Lösungen existieren.

Wir betrachten hierzu EG-Fluktuationen mit Standardabweichung $\sigma$ und Korrelationszeit $\tau_c$.

  • Freie Diffusion:

Für flaches Potential erwartet man zwei Regime in der MSD:

Ballistisches Regime für $\tau<<\tau_c$ mit $\overline{x^2}(\tau)= \left( \frac{\sigma^2}{\gamma^2} \right) t^2$.

Diffusives Regime für $\tau>>\tau_c$ mit $\overline{x^2}(\tau)= 2\left( \frac{\sigma^2 \tau_c}{\gamma^2} \right) t$.

Auf langen Zeitskalen $\tau>>\tau_c$ erscheint die kurzzeit-korrelierte Zufallskraft wie weißes Rauschen

(4)
\begin{align} <F(t_1)F(t_2)>=\Gamma_F \;\delta(t_1-t_2) \end{align}

mit einer effektiven Fluktuations-Stärke

(5)
\begin{align} \Gamma_F = 2 \sigma^2 \tau_c. \end{align}
  • Diffusion im harmonischen Potential:

Wir betrachten hier nur den Fall eines relativ weichen Potentials (kleine Federkonstante $k$), so daß die charakteristische Übergangszeit

(6)
\begin{align} \tau_k = \frac{\gamma}{2k} >> \tau_c \end{align}

weit über der Korrelationszeit der Zufallskraft liegt. In diesem Fall kann man die Zufallskraft als weißes Rauschen nähern. Man erwartet dann wiederum zwei Regime:

Diffusives Regime für $\tau<<\tau_k$ mit $\overline{x^2}(\tau)= \left( \frac{\Gamma_f}{\gamma^2} \right) t = \left( \frac{2 \sigma^2 \tau_c}{\gamma^2} \right) t$.

Plateau-Regime für $\tau>>\tau_k$ mit $\overline{x^2}(\tau)= \frac{\Gamma_f}{2k\gamma} = \frac{\sigma^2 \tau_c}{k\gamma}$.


[4] A.M.: Alpha-Modell:Erste Testergebnisse

In [3] bei Punkt freie Diffusion fehlt in der Formel für das diffusive Regime ein Faktor 2. Außerdem ist wahrscheinlich im zweiten Teil der Formel (2) P(F) gemeint.
Die Wahrscheinlichkeits Verteilung der Zufallskraft entspricht der erwarteten Gaußform (2) mit den entsprechend eingestellten Werten für die Varianz. Bei der Korrelation ergibt sich eine gute Übereinstimmung mit Formel (2) nur für $\tau_c=1$, für zwei größere Werte stimmt zwar die Steigung im logarithmischen Plot, die Kurven sind jedoch verschoben. Das Integral über die Korrelationsfunktion wurde konstant auf N=1 gehalten. (siehe file)
Die freie Diffusion wurde für vier verschiedene Werte von $\tau_c$ mit N=10 simuliert, ab $\tau_c=10$ tritt auch das ballistische Regime auf. (siehe file F_ext=0.dat)
Für das harmonische Potential wurde bei drei verschiedenen Korrelationszeiten die Abhängigkeit von der Federkonstante untersucht (N=10=const.). Hierbei wurde gemäß Gleichung (6) beachtet, dass die Federkonstante nicht zu groß wird. Die Plateauwerte sind in allen Simulationen um einen Faktor 2 zu hoch. Ab $\tau_c=10$ tritt wieder das ballistische Regime auf. (siehe files tau_c=1_3k.dat, tau_c=1e-1_3k.dat, tau_c=10_2k.dat)
Die Wahrscheinlichekeitsverteilungen für x sind gaußförmig, für zu kleine Federkonstanten treten aber Abweichungen auf.


[5] C.M.: Re: Alpha-Modell:Erste Testergebnisse

  • In [3] bei Punkt freie Diffusion fehlt in der Formel für das diffusive Regime ein Faktor 2. Außerdem ist wahrscheinlich im zweiten Teil der Formel (2) P(F) gemeint.

Habe ich verbessert.

  • Die Wahrscheinlichkeits Verteilung der Zufallskraft entspricht der erwarteten Gaußform (2) mit den entsprechend eingestellten Werten für die Varianz.

Ausgezeichnet.

  • Bei der Korrelation ergibt sich eine gute Übereinstimmung mit Formel (2) nur für $\tau_c=1$, für zwei größere Werte stimmt zwar die Steigung im logarithmischen Plot, die Kurven sind jedoch verschoben. Das Integral über die Korrelationsfunktion wurde konstant auf N=1 gehalten. (siehe file)

Das klingt nach einem Normierungsproblem. Die Auto-Korrelationsfunktion der Kraft sollte doch bei $\tau=0$ den Wert $\sigma^2$ annehmen. Tatsächlich schneiden sich die simulierten Kurven ja bei $\tau=0$ im Wert $\sigma^2=1$. Vielleicht war die Kraft-Varianz im Programm ja wirklich immer auf 1 eingestellt ? Vielleicht könntest Du den verwendeten cpp-File auch einmal uploaden. Es könnte aber auch sein, daß beim Abspeichern der numerisch bestimmten Autokorr.-Funktion künstlich auf Maximalwert 1 normiert wird. Man müßte in der entsprechenden Klassen-Methode nachsehen:

class DRecXY {
…..
void corrCoeff(int DTSTP, char* froot, int PAR) {
…..
do {
sum=0.0;
for (N=0;N<=(NMX-DT);N++) {
sum+=(ff[N]*ff[N+DT]);
}
double cor=sum/((double)(NMX-DT+1));
double cc=cor/var;
corco.record(dt*(double)DT,cc);
if (DTSTP>0) {
DT+=DTSTP;
} else {
DT*=2;
}
} while(DT<DTMAX);

In der Tat, die Routine berechnet den Korrelations-Koeffizienten, also die Autokorrelation geteilt durch die Varianz. Du kannst entweder die fettgedruckte Zeile ersetzen durch

double cc=cor*dt;

, oder aber normierte Fits verwenden.

  • Die freie Diffusion wurde für vier verschiedene Werte von $\tau_c$ mit N=10 simuliert, ab $\tau_c=10$ tritt auch das ballistische Regime auf. (siehe file F_ext=0.dat)

Sehr schön.

  • Für das harmonische Potential wurde bei drei verschiedenen Korrelationszeiten die Abhängigkeit von der Federkonstante untersucht (N=10=const.). Hierbei wurde gemäß Gleichung (6) beachtet, dass die Federkonstante nicht zu groß wird. Die Plateauwerte sind in allen Simulationen um einen Faktor 2 zu hoch. Ab $\tau_c=10$ tritt wieder das ballistische Regime auf. (siehe files tau_c=1_3k.dat, tau_c=1e-1_3k.dat, tau_c=10_2k.dat)

Das ist sehr interessant. Deutet auf einen analytischen Rechenfehler hin. Wir sollten mit folgender (nachgeprüften) Herleitung aus einer meiner Publikationen vergleichen:

http://lpmt-theory.wikidot.com/local--files/anja-projekt/msdWell

Die Reibungskonstante ist dort mit $\alpha$ bezeichnet, außerdem sind Ort und Kraft zwei-dimensionale Vektoren. Geht man wieder auf 1D zurück, entfällt der Faktor 2 in Formel (B8). Somit muß man die Formeln (B9) sowie das Endresultat (B10) mit einem Faktor 1/2 multiplizieren. Dies ergibt dann also

(7)
\begin{align} \overline{x^2}(\tau)= \frac{\Gamma_f}{k\gamma}\;\left[ 1-e^{\tau/\tau0} \right], \end{align}

was zu Deinen Simulationen besser passen dürfte.

Der Unterschied dieser Herleitung zu jener aus der TGBP-Vorlesung

http://lpmt-theory.wikidot.com/local--files/anja-projekt/DiffCorrForce

ist lediglich, daß die Zeit-Integrale von $-\infty$ statt von $0$ loslaufen.

Dann sollte aber der gleiche Faktor 2 in den Formeln für die freie Diffusion mit exp-korrelierter Antriebskraft falsch sein. Könntest Du vielleicht diese Herleitung nochmal korrekt durchführen ? Man müßte im TGBP-Skript schon auf Seite 56 unten die beiden Integrationen von $-\infty$ starten lassen.


[6] A.M.: Analytische Herleitung der MSD für harmonisches Potential

  • In der Tat, die Routine berechnet den Korrelations-Koeffizienten, also die Autokorrelation geteilt durch die Varianz.

Dies sollte also keine Auswirkungen auf die Simulationsresultate haben.

  • Man müßte im TGBP-Skript schon auf Seite 56 unten die beiden Integrationen von $-\infty$ starten lassen.

Habe ich versucht, es tritt aber ein divergierendes Integral auf:

(8)
\begin{align} <x^2>= \frac{<F^2>}{\gamma^2}(2\tau_c \int_{-\infty}^{t}dt_1-\tau_c^2\exp(\frac{2t}{\tau_c})) \end{align}

Eigentlich betrifft die Herleitung auf S. 56 ja auch die freie Diffusion, wo die Vorfaktoren noch mit der Simulation verträglich waren. Gleichung (7) sollte sich vielleicht eher mit der Herleitung ab S. 59 im TGBP-Skript ergeben. Werde das gleich mal ausprobieren.
Ergebnis:
Ersetzt man in der Herleitung einfach die untere Integrationsgrenze in den Zeitintegralen durch $-\infty$ ergibt sich $<x^2>=\frac{\Gamma_F}{2k\gamma}$, was (bis auf den Faktor 1/2 aufgrund der unterschiedlichen Dimensionen) mit dem Wert für c(0) aus dem Paper übereinstimmt. Im TGBP_Skript sind die oberen Integrationsgrenzen ja auch identisch, sodass dies konsistent ist. Die Herleitungen unterscheiden sich also außerdem noch in der "Definition" der MSD. Mit Integration ab $-\infty$ und der Formel für die MDS aus dem Paper bekomme ich auch das gleiche raus.


[7] A.M.: Ergebnisse für Gaußpotential:

Bei N=1 wurde für $\tau_c=1$ die MSD bei konstanter Breite w von 1 Variation der Tiefe des Potentialtopfes untersucht. Mit steigender Tiefe nimmt der Sättigungswert ab. Der Exponent im Anfangsbereich beträgt ungefähr 0.5. (depth_fit_n=1,tau_c=1.dat)
Für $\tau_c=10$ ergibt sich für die Tiefenabhängigkeit bei $w=1$ im Sättigungsbereich ein ähnliches Bild. Der Exponent beträgt hier jedoch 1 wie im diffusiven Regime. (depth_fit_n=1,tau_c=1e1.dat)
Die Wahrscheinlichkeitsverteilung für x ist näherungsweise gaußförmig. (p(x)_depth_fit.dat)
Bei Variation der Breite bei konstanter Tiefe d von 10 zeigt sich ein Anstieg des Sättigungswertes mit steigender Breite. Die Exponenten ändern sich hierbei allmählich von 1 (für w=1) auf 2(für w=10). (width_fit_n=1,tau_c=1e1.dat)
Für N=10 tritt für zu kleine Werte für $\tau_c$ aufgrund der höheren Varianz für die Kraft kein Sättigungsregime in der MSD auf. Für $\tau_c=1$ geht der Exponent in der MSD nach einem kurzen Bereich mit höherem Wert in 1 über. Der genaue Anfangswert lässt sich aufgrund der Kürze nicht genau bestimmen. (n=10,tau_c=1_fit.dat)
Die Wahrscheinlichkeitsverteilung ähnelt der für freie Diffusion. (n=10,tau_c=1_p(x).dat)
Bei $\tau_c=2$ ist er zu Beginn ungefähr 2 und es tritt später in relativ guter Übereinstimmung mit der kritischen Zeitkonstanten der Übergang in das diffusive Regime auf. (n=10,tau_c=2_fit.dat)
In der Wahrscheinlichkeitsverteilung tritt schon ein einzelnes deutliches Maximum auf. (n=10,tau_c=2_p(x).dat)


[8] C.M.: Re: Ergebnisse für Gaußpotential:

  • N=1

Ich finde einige Deiner Resultate ziemlich spannend, insbesondere den anscheinend durchstimmbaren, fraktionierten Exponenten in der MSD, der in Abbildung width_fit_n=1,tau_c=1e1.dat vor dem Plateau auftritt.

Ist dieses Verhalten eine besondere Eigenschaft des Gaußpotentials ? Oder würde es auch schon in einem harmonischen Potential auftreten ? Wenn man den Gauß $U(x)=-h\; e^{-\frac{1}{2}(x/w)^2}$ für kleine Abweichungen vom Minimum taylor-entwickelt, erhält man (wenn ich richtig gerechnet habe) $U_{tay}(x)=\frac{h}{2w^2} x^2$, was einer effektiven Steifigkeit $k=\frac{h}{w^2}$ entspricht. Man könnte also einmal eine Gauß-MSD direkt auf die entsprechende Parabel-MSD legen.

Sollte die fraktionale Diffusion auch im harmonischen Topf funktionieren, wäre eventuell sogar eine analytische Theorie möglich. Man könnte wieder die Greensfunktion für das visko-elastische Teilchen nehmen und damit die MSD als Integral über die Kraft-Autokorrelations-Funktion ausdrücken. Nun würde man aber nicht wie bisher eine $\delta-$Korrelation, sondern eben eine exponential-korrelierte $C_{FF}(\tau) = \sigma^2 e^{-\tau/\tau_c}$ einsetzen. Die ganzen exp-Funktionen sollten ja analytisch integrierbar bleiben. Wenn man die analytische MSD für den allgemeinen Fall erstmal hätte, könnte man vielleicht für kleine $\tau$ zeigen, daß sie wie ein Potenzgesetz verläuft.

Es wäre weiterhin interessant, ob man den fraktionierten Übergangsbereich zeitlich weiter ausdehnen kann, z.B. indem man die Korrelationszeit der Kraft noch größer macht.

  • N=10

Es ist auch super, daß man im Multi-Topf ein strukturlos diffusives Resultat erhält. Man kann offenbar die ganze komplexe Potential-Landschaft im Sinne eines "Coarse Grainings" als glattes Medium mit effektiver Viskosität ansehen ! Ein wenig Sorge macht mir der Rand des Systems, aber mit den P(x)-Verteilungen kann man ja sehen, wie oft sich das Teilchen außerhalb aller Töpfe aufhält. Die Tiefe aller Töpfe war ja noch identisch gewählt, oder ? Vielleicht könnte man stattdessen ein Einzel-Sinuspotential $U(x)=a \;\sin(k_mx)$ verwenden.


[9] A.M.: Verwechslung der Bezeichnung N

Mit N habe ich nicht die Anzahl der Töpfe, sondern $\tau_c \sigma^2$ bezeichnet. Deshalb ist für $tau_c=1$ die Varianz der Kraft im Vergleich zum Fall $N=1$ höher, sodass das Teilchen leichter aus dem Potentialtopf entkommen kann.


[10] C.M.: Überdämpftes Teilchen im harmonischen Potential mit exp-korrelierten Zufallskräften

Anja hat die MSD für den allgemeinen Fall berechnet. Dabei hat sich folgende Vorgehensweise bewährt:

(9)
\begin{align} \overline{x^2}(t) = 2\left( <x(0)x(0)> -<x(t)x(0)> \right), \end{align}

wobei

(10)
\begin{align} x(t)=\int_{-\infty}^t G(t-t^{\prime}) F(t^{\prime}) dt^{\prime}, \end{align}

so daß

(11)
\begin{align} <x(t)x(0)> = \int_{-\infty}^t dt_1 \int_{-\infty}^0 dt_2 \;\;G(t-t_1)G(t-t_2) <F(t_1)F(t_2)>. \end{align}

In unserem Fall lautet die Greensfunktion

(12)
\begin{align} G(t-t^{\prime}) = \;\Theta(t-t^{\prime})\; \frac{1}{\gamma} e^{-(k/\gamma)(t-t^{\prime})} \end{align}

und der Kraft-Korrelator

(13)
\begin{align} <F(t_1)F(t_2)> = \sigma^2 e^{-|t_1-t_2|/\tau_c}. \end{align}

Wir definieren eine Dekorrelationsrate $\nu_c=\frac{1}{\tau_c}$ und eine viskoelastische Relaxationsrate $\nu_k=\frac{k}{\gamma}$. Die resultierende MSD im allgemeinen Fall lautet dann

.

(14)
\begin{align} \overline{x^2}(t) = \left[ \frac{2\sigma^2}{\gamma^2 (\nu_k^2-\nu_c^2)} \right]\;\left[ (1-e^{-\nu_ct})-\frac{\nu_c}{\nu_k}(1-e^{-\nu_kt})\right]. \end{align}

.

Das Bild unten zeigt den Verlauf für $\tau_c=0.01$ und $\tau_k=10$:

msd3Regimes.jpg

Im Grenzfall kurzer Zeiten (verglichen sowohl mit $\tau_c$ als auch mit $\tau_k=1/\nu_k$) kann man die exp-Funktionen bis zur zweiten Ordnung taylorentwickeln und erhält eine ballistische MSD

(15)
\begin{align} \overline{x^2}(t \rightarrow 0) = \frac{\sigma^2\nu_c}{\gamma^2 (\nu_c+\nu_k)}\;t^2. \end{align}

Im Grenzfall veschwindenden Potentials $k \rightarrow 0$ und bei nicht zu großen Zeiten $t << \tau_k$ kann man nur die $\nu_k$-abhängige exp-Funktion entwickeln und erhält eine zunächst ballistische, später diffusive MSD

(16)
\begin{align} \overline{x^2}(k \rightarrow 0) = \frac{2\sigma^2\tau_c^2}{\gamma^2}\;\left[ (t/\tau_c) - (1-e^{-t/\tau_c}) \right]. \end{align}

Dies stimmt überein mit dem Resultat aus TGBP Seite 58.

Der Grenzfall kurzer Korrelationszeit $\tau_c << t , \tau_k$ ist ebenfalls von Interesse. Man kann dann den Term $e^{-\nu_ct}$ vernachlässigen und weiterhin ausnutzen, daß $\nu_c>>\nu_k$. Es ergibt sich

(17)
\begin{align} \overline{x^2}(\tau_c \rightarrow 0) = \left( \frac{(2\sigma^2\tau_c)}{k\gamma} \right)\;\left[ 1-e^{t/\tau_k} \right]. \end{align}

Mit der effektiven weißen Rausch-Stärke $\Gamma_F = 2\sigma^2\tau_c$ stimmt dies mit Formel (7) oben überein, nicht aber mit der Formel aus TGBP Seite 60.


[11] A.M.: Harmonische Näherung des Gaußpotentials

Für den Fall eines Gaußpotentials mit $w=4$, $d=10$ wurde ein harmonisches Potential mit $k=0.625$ in den Plot eingefügt. Die beiden Kurven stimmen sehr gut überein, insbesondere ergibt sich auch ein fraktionierter Exponent. Lediglich die Sättigungswerte weichen geringfügig ab.
Zum Vergleich wurde auch die analytische Kurve gemäß Gleichung (14) geplottet. Diese stimmt im Sättigungsbereich wie erwartet besser mit der harmonischen Näherung überein. Der Anfangswert ist im Vergleich zu den simulierten Kurven niedriger, sodass sich im Anfangsbereich auch ein höherer Exponent von ungefähr 2 ergibt. (vgl_theo_gauss_harmon.dat)


[12] C.M.: Re: Harmonische Näherung des Gaußpotentials

Nachdem die Gauß-Resultate so gut zu den harmonischen, und diese zu der analytischen MSD passen, muß man nach meiner Meinung eigentlich davon ausgehen, daß die MSD für sehr kleine Zeiten ballistisch ist. Der fraktionale Exponent kann dann wohl nur näherungweise in einem Übergangsregime auftreten. Ich habe aber im Augenblick keine Idee, wie man so etwas analytisch aus der allgemeinen MSD-Formel herausdestillieren könnte.

Wenn man alle Parameter gleich läßt, aber die Varianz $\sigma^2$ der Rauschkraft erhöht, sollte das Teilchen doch irgendwann die Ausläufer des Gaußpotentials erkunden und es müßten sich Abweichungen zum harmonischen Fall zeigen.


[13] C.M.: Logarithmische Ableitung

Um die Frage des fraktionalen MSD-Exponenten im Rahmen der analytischen Formel zu klären, habe ich den lokalen Exponenten als logarithmische Ableitung der MSD bestimmt:

(18)
\begin{align} e(t)= \frac{\partial \log_{10}(\overline{x^2}(t))}{\partial \log_{10}(t)} \end{align}

Für die gleichen Paramter ($\tau_c=0.01$ und $\tau_k=10$ wie bei der Abbildung in Post [10] findet man:

msdExp

Offensichtlich ergibt sich bei diesen speziellen Parametern kein klares Plateau bei einem gebrochenzahligen Exponenten.

Man könnte noch Hoffnungen hegen, wenn sich die beiden charakteristischen Zeitkonstanten sehr ähnlich werden (Raten: $\nu_k=1$, $\nu_c=1.01$). Aber auch dann:

msdExp2

[14] C.M. Statistik des Well-to-well Transfers

Nachdem bereits die Hälfte der Projektwoche um ist (-: sollten wir uns vielleicht langsam von den (allerdings sehr wichtigen) Testfällen lösen und unser eigentliches Problem angehen: Wie muß der selbstgetriebene Agent die Korrelationen seiner Kraft wählen, damit er möglichst schnell von einem Potentialtopf in den nächsten kommt ?

In Anbetracht der Probleme mit der "FindThreshCross"-Funktion schlage ich folgendes Computer-Experiment vor: Überdämpftes Teilchen mit exp-korrelierten Kräften im Doppelgauß-Potential.

Wenn die zentrale Barriere bei $x=0$ liegt, läßt sich eine Barrierenüberschreitung sehr leicht durch den Vorzeichenwechsel von $x(t)$ detektieren. Die Zeitintervalle $\Delta t_i = t_{i+1}-t_i$ zwischen zwei aufeinanderfolgenden Überschreitungs-Zeitpunkten $t_k$ wird histogramiert und $P(\Delta t)$ bestimmt.

Für die Zufallskraft würde man die effektive weiße Rausch-Stärke $\Gamma_F = 2\sigma^2\tau_c$ konstant halten und die Korrelationszeit $\tau_c$ verändern. Welchen Einfluß hat dies auf $P(\Delta t)$ ?

.

Man muß bei den Simulationen u.a. auf folgende Punkte achten:

  • Die beiden äußeren Barrieren des Doppelgaußpotentials müssen so hoch sein, daß das Teilchen niemals daraus entkommen kann.
  • Es müssen natürlich hinreichend viele Well-to-well-Transfers stattfinden können, damit ein Histogramm überhaupt Sinn macht. Dies erfordert ggf. lange Laufzeiten.

.

Die neue Funktion befndet sich im File potdif1.cpp.

.

Das Programm enthält folgende zusätzliche Neuerungen:

  • Bei NWELL=2 werden die beiden Töpfe symmetrisch um 0 plaziert.
  • Der Abstand der Töpfe steht in der Variablen wDist.
  • Die effektive Rauschstärke ist mit rGamma=2.0 so gewählt, daß sich bei einer Korrelationszeit rTau=1.0 genug Transfers ergeben.
  • Die Standardabweichung rSigma wird dann automatisch gesetzt.
  • Die neue Routine SignChangeStatistics wird in ComputeTrajectory aufgerufen und speichert $P(\Delta t)$ ab. Dabei druckt sie noch die Anzahl der Transfers KMX aus.
  • Einige unnötige Dinge wurden ausdokumentiert.

Der erste Test zeigt, daß $P(\Delta t)$ sehr empfindlich auf rTau reagiert. Schon bei 10 Prozent Änderung der Korrelationszeit sieht man Effekte in $P(\Delta t)$.


[14] A.M.: Re: Harmonische Näherung des Gaußpotentials

Für $\tau_c=10$ und $w=1$, $d=10$ wurde nun $\sigma^2$ variiert. Bis zu einem Schwellwert von 1.36 treten kaum Abweichungen zwischen harmonischem Potential und Gaußpotential auf. Für höhere Varianzen ändert sich der qualitative Verlauf der MSD für das Gaußpotential stark. Bei weiterer Erhöhung verschwindet der Sättigungsbereich schließlich ganz. (gauss_harmon_sigma.dat)

Super !


[15] A.M.: Re: Statistik des Well-to-well Transfers

Sollen dann neben verschiedenen Werten für $\tau_c$ auch unterschiedliche $\Gamma_F$ getestet werden?

Ist sicher auch interessant. Muß jetzt leider ins Labmeeting )-:

$P(\Delta t)$ scheint exponentiell abzufallen.

Wenn $P(\Delta t)$ exp. abfällt, dann wäre der Transfer ein Poisson-Prozeß mit einer festen mittleren Transferzeit $\tau_T$ , was schön ist. Das spannende ist aber, daß der Zusammenhang von $\tau_T$ mit den System-Parametern, insbesondere mit $\tau_c$, auf den ersten Blick nicht leicht vorhersagbar ist…

Für $\Gamma=2.0$ steigt $\tau_T$ bei Erhöhung von $\tau_c$ in nichtlinearer Weise an. (6tau_fit_.dat)

Sehr interessant ! Soweit wir jetzt sehen, ist die beste Strategie für den Agenten, kurze Korrelationszeiten zu verwenden, weil dann die mittlere Transferzeit kürzer wird. Wegen der Nichtlinearität wird das zunehmend schwieriger.

Bei Erhöhung von $\Gamma$ fällt die mittlere Transferzeit $\tau_T$ ab. (3gamma_fit.dat)

OK, das macht Sinn: Höhere Rausch-Amplitude —> Schnellerer Transfer. Ich würde fast erwarten, daß $\tau_T$ exponentiell mit dem Verhältnis aus Barrierenhöhe und $\Gamma$ zusammenhängt… Vielleicht …

In der Tat steigt $\tau_T$ exponentiell mit $\frac{B}{\Gamma}$an. (Tau_T_fit.dat)
Auch die Abhängigkeit von $\tau_T$ und $\tau_c$ kann man gut mit einem exponentiellen Anstieg fitten. (tau_t_vs_tau_c_fit.dat)

Großartig !


[16] A.M.: Doppelgaußpotential ohne Barriere

Bei einem Abstand von 2.0 der symmetrisch um den Ursprung liegenden Gaußpotentiale mit $w=1$, $d=10$ schien das resultierende Gesamtpotential in der Mitte keine Barriere zu haben, was scheinbar im Widerspruch zu $P(x)$ stand, wo zwei Maxima erkennbar waren. Das Potential hat in der Tat bei $x=0$ ein Minimum und sonst keine weiteren Extrema. Es ist allerdings sehr flach, sodass zum Vergleich ein Potential bestehend aus einem konstanten Teil in der Mitte und zwei linear ansteigenden Ästen als Begrenzung simuliert wurde. Auch bei diesem zeigten sich zwei Maxima in $P(x)$, deren Position gut mit den Unstetigkeitsstellen des Potentials übereinstimmt. Beim Doppelgaußpotential ist der Übergang kontinuierlich, sodass die Maxima nach innen verschoben sind. (landscape_p(x)_vgl_eckig.dat)

OK, dann ging das schon mit rechten Dingen zu !


[17] C.M. Brainstorming

Im folgenden einige unstrukturierte Bemerkungen und unfertige Ideen ("schriftliches Denken") zu unserem Projekt:

  • Ich denke, für Deinen Bericht hast Du eigentlich jetzt schon genug Material. Wir brauchen also von nun an nicht mehr möglichst viele Resultate produzieren, sondern können uns wieder mehr vom Fernziel der Zellmigration leiten lassen. Die Modelle brauchen dann nicht mehr so generisch allgemein und analytisch zugänglich sein wie bisher.
  • Es wäre aber wünschenswert, daß die neuen Modelle auf den Erkenntnissen aufbauen, die wir aus den alten gewonnen haben. Insbesondere sollten wir den Transfer im Doppeltopf im Sinne der Zellmigration betrachten.
  • Den exponentiellen Anstieg der Transfer-Wartezeiten mit der Kraft-Korrelationszeit habe ich noch nicht wirklich verstanden. Bei einer korrelierten Kraft können ja mehrere ähnlich große - und insbesondere mit dem gleichen Vorzeichen behaftete - Kraftwerte unmittelbar aufeinander folgen, während bei einem weißen Rauchen in jedem Zeitschritt ganz neu gewürfelt wird. Meine ursprüngliche naive Intuition war, daß sich die hiermit verbundenen räumlichen Verschiebungen des Teilchens im Potential irgendwie kohärent aufaddieren würden und das Teilchen so leichter die Grenze der Barriere zum Nachbartopf erreicht. Wäre dem so, würde die Transfer-Wartezeit mit der Kraft-Korrelationszeit sinken, nicht steigen.
  • Kommt diese Diskrepanz einfach von der überdämpften Dynamik unseres Teilchens ? Die Auslenkung des Teilchens von seinem (linken oder rechten) Topfminimum ist ja in jedem Augenblick nur durch den Momentanwert der Kraft bestimmt. Es gibt da also nicht wirklich ein Gedächtnis, durch das sich aufeinanderfolgende Kräfte aufaddieren könnten (Nur immer nach jedem Topf-zu-Top-Transfer führt die gleiche Momentankraft plötzlich zu einer anderen Position). Was die stationäre Verteilung $P(F)$ der Momentankräfte angeht, liefert unser Box-Müller-Verfahren ja immer brav seine Gaußverteilungen. Wenn wir bei konstantem $\Gamma_F$ die Korrelationszeit $\tau_c$ erhöhen, geht aber $\sigma^2$ in die Knie, d.h. $P(F)$ wird schmäler und erreicht die nötige Grenze zum Nachbartopf nicht mehr so leicht ! Hinzu kommt noch, daß bei korrelierten Kräften (für Barrienüberwindungen) zu kleine $F(t)$ unnötigerweise mehrfach hintereinander auftreten - was Zeit kostet. Könnte das die Erklärung sein ?

Wäre durchaus möglich. Sollte man zum Test einmal Simulationen mit konstantem $\sigma^2$ bei gleichzeitiger Änderung von $\tau_c$ durchführen?

Sehr gute Idee !

Leider scheint sich die Hypothese nicht zu bestätigen. $\tau_T$ steigt auch bei konstantem $\sigma^2$ weiter mit $\tau_c$ an. (sigma=1,tau.dat)

  • Ursprünglich ging ich eigentlich auch von einem 2D-Modell mit Spinnen-Geometrie aus: Die Modell-Zelle erzeugt ihre Antriebskraft durch M Stressfasern, die von ihrem Zentrum aus radial in verschiedene Richtungen $\vec{e}_m$ zeigen. Die Gesamtkraft ist die Vektor-Summe der gerichteten Einzelkräfte $f_m(t)$ jeder Faser, $\vec{F}(t)=\sum_{m=1}^M \vec{e}_m f_m(t)$. Ich stelle mir für die Einzelfaserkräfte $f_m(t)$ Ketten von positiv definiten Kraftpulsen verschiedener Länge vor, z.B. Ketten von Sinus-Halbbögen mit jedesmal anderer zeitlicher Länge $T_k$. Was ist die optimale Verteilung $P(T)$ dieser Pulslängen im Hinblick auf effiziente Bewegung in einer 2D Potentiallandschaft ?
  • Im obigen 2D-Spinnen-Modell besteht ja offensichtlich die Möglichkeit, daß sich die Einzelfaserkräfte kohärent überlagern und so kurzzeitig sehr hohe Gesamtkräfte auftreten können. Ein Teil der Fasern kann die Zelle schon halb den Potentialberg hinaufziehen und sie dort eine Zeit lang halten. Dann braucht nur noch eine weitere Faser hinzukommen und die Zelle über die Barriere bringen. Macht das nach Deiner Meinung irgendeinen Sinn ?

Klingt für mich plausibel.

  • Eine andere Intuition, die verschiedene Leute zu haben scheinen, ist die, daß eine superbreite Powerlaw-Verteilung von Zeitkonstanten $P(T)$ besonders günstig ist, weil sie damit "für alle denkbaren Fälle gerüstet" ist.

[18] C.M. Auf zu explorativen, rückgekoppelten Migrationsstrategien im 2D

Ich bin nach sorgfältigem Abwägen zur Überzeugung gelangt, daß wir die blinden Migrations-Strategien hinter uns lassen und gleich einen Sprung machen sollten in Richtung auf explorative, rückgekoppelte Modelle, die in Hinblick auf lebenden Zellen wirklich interessant sind. Nur die kämen dann ja auch für eine Publikation infrage. Mir schweben da einige sehr spannende Modelle vor …

Zwar ist die Projektwoche ist nun schon fast vorbei und das Schreiben der neuen Simulations-Programme wird einen signifikanten Teil des heutigen Tages kosten. Aber nachdem inklusive Vor-und Nachbereitung ja anscheinend 3 Wochen zur Verfügung stehen, könntest Du ja z.B. schon mit dem Schreiben des Berichtes beginnen (bitte auch noch die paar Arbeiten von gestern dokumentieren), während ich heute eine Rohversion des neuen Programmes zusammenstricke. Schlimmstenfalls müßten wir auch noch nächste Woche ein paar Simulationen laufen lassen. Wäre das akzeptabel ?


[19] C.M.: Lambda-Modell

Wir behalten den äußeren Rahmen des Alpha-Modells bei, also punktförmiger Agent, Hintergrundviskosität, externe Potential- und intern generierte Kräfte. Wieder seien die internen Kräfte für einen Zeitraum $\Delta t_{sim}$ konstant. Allerdings gehen wir nun zu zwei Raumdimensionen über, so daß die Bewegungsgleichung lautet

(19)
\begin{align} \frac{d}{dt}\;\vec{R} = \left( \vec{F}(t) - \frac{\partial}{\partial \vec{R}}U(\vec{R}) \right)/\gamma. \end{align}

Als Potentiallandschaft kommt z.B. ein 2D-Multigauß infrage.

Der wesentliche Unterschied des Lambda-Modells zum Alpha-Modell besteht in der explorativen, rückgekoppelten Dynamik der internen Kräfte $\vec{F(t)}$. Der Name "Lambda-Modell" soll diesen großen qualitativen Sprung ausdrücken.

Der Agent besitzte $S_{mx}$ "Stressfasern", die innerhalb des Zeitraums $\Delta t_{sim}$ sowohl hinsichtlich ihrer Zugrichtung $\vec{e}_s$ (Einheitsvektor), als auch hinsichtlich ihrer Zugkraft $f_s=f_0$ konstant sind (Also kein kontinuierliches An- und Abklingen der Kraft mehr). Die auf den Agenten wirkende Gesamtkraft ist wieder die Summe der Einzelfaserkräfte:

(20)
\begin{align} \vec{F}(t_n) = f_0 \sum_{s=1}^{S_{mx}} \vec{e}_s(t_n). \end{align}

Die Simulation zerfällt in kontinuierliche Abschnitte $t \in [t_n,t_{n+1}[$, in denen die $\vec{e}_s$ konstant sind und die Bewegung des Agenten durch numerisches Lösen der Langevin-Gleichung für festes $\vec{F}(t_n)$ im 2D-Potential $U(\vec{R})$ bestimmt wird. Zu den diskreten Zeitpunkten $t_n = n \Delta t_{sim}$ aber werden die Zugrichtungen $\vec{e}_s$ synchron aktualisiert. Dies geschieht nach folgender Regel:

Für jede Faser $s$ wird die Komponente $\Delta L_s$ bestimmt, um den sich der Agent im letzten Zeitintervall $\Delta t_{sim}$ relativ zur alten Zugrichtung $\vec{e}_s(t_{n-1})$ verschoben hat, also das Skalarprodukt

(21)
\begin{align} \Delta L_s = (\vec{R}(t_n)-\vec{R}(t_{n-1})) \cdot \vec{e}_s(t_{n-1}). \end{align}

Falls die Verschiebungskomponente groß genug war, $\Delta L_s > \Delta L_s^{(min)}$, bleibt die Zugrichtung der Faser unverändert, $\vec{e}_s(t_{n})=\vec{e}_s(t_{n-1})$. Falls nicht, wird die Faser aufgelöst und eine völlig zufällige neue Richtung gewählt, $\vec{e}_s(t_{n})=\mbox{random}$. Auf die gleiche Weise wird mit allen $S_{mx}$ Fasern verfahren. Die nützlichen Fasern (die viel zur Gesamtverschiebung der Agenten beigetragen haben) behalten ihre Zugrichtung bei, die nicht nützlichen probieren eine neue Zugrichtung aus.


[20] C.M.: Erste Vorstufe: Lambda1

Ich habe als Vorübung das Alpha-Programm auf 2D erweitert (lambda1.zip).

Dabei waren einige Veränderungen notwendig:

  • Diverse Orts-und Kraft-Variablen sind jetzt vom Typ vec2D statt double.
  • Da die Teilchen bei stärkerem Rauschen oft aus dem Simulationsbereich entkommen, wurden ein äußeres Confinement-Potential eingeführt. Dieses steigt jenseits der Grenzen jeweils rampenförmig an mit einer wählbaren Steifheit kConf.
  • Die Ausgabe traj.dat enthält jetzt die x- und y-Koordinaten des Teilchens nach jedem Zeitschritt dT. Die Zeitpunkte selbst werden nicht mehr gespeichert.
  • Die graphische Darstellung ist nun natürlich umständlicher (Prozedur StoreData). Die 2D Potentiallandschaft wird zunächst mit NMX*NMX Punkten abgerastert und als 2D-Matrix land2D.dat gespeichert. Manche Grafik-Programme (sogar Easyplot) können daraus einen Konturplot anfertigen. Allerdenings möchte man die Teilchentrajektorie gerne darüberlegen, was sich als nicht-trivial erweist.
  • Daher folgender Ausweg: Es wird auch eine Datei obstacles.dat gespeichert. Dabei handelt es sich um zufällig verteilte Punkte (Hindernisse für die Zelle), deren Dichte mit der Höhe des lokalen Potentials zunimmt. Man muß diese Punkte mit Easyplot erst einmal im Style "Mark Pts." (ohne "Connect Pts.") speichern. Die weißen Bereiche sind die Potential-Minima.
  • Dann kann man die Trajektorie traj.dat im Style "Connect Pts." (ohne "Mark Pts.") darüber zeichnen.
  • Man darf sich nicht von dem Querformat irritieren lassen, mit dem Easyplot ggf. auch quadratische Simulationsbereiche zeichnet.

[21] A.M.: Erste Ergebnisse für Lamda1-Modell

In einer festen Zufallslandschaft wurde bei festem $\Gamma=20$ die Korrelationszeit variiert. Wie im 1D-Fall kommt das Teilchen bei höheren Korrelationszeiten schlechter voran. Für $\tau_c=1$ und $\tau_c=10$ nimmt das Teilchen zu Beginn noch den gleichen Weg, bei $\tau_c=100$ landet es in einem anderen Potentialtopf, aus dem es auch nicht mehr herauskommt.
Nur bei $\tau_c=1$ führt eine Verlängerung der Simulationszeit zu einer Fortsetzung der Trajektorie in noch nicht besuchte Gebiete. Für die längeren Korrelationszeiten bleibt das Teilchen im zuerst erreichten Potentialtopf gefangen, weswegen auf ein Hochladen dieser files Speicherplatzgründen verzichtet werden kann. (tmx=1e4,tau=1.dat, tmx=1e4,tau=1e1.dat, tmx=1e4,tau=1e2.dat, tmx=1e5,tau=1.dat)
Nun bietet sich die Untersuchung von variablen Korrelationszeiten bei festem $\sigma^2$ an.


[22] C.M.: Zweite Vorstufe: Lambda2

Die nächste Version lambda2 ist fertig:

  • Der Agent hat bis jetzt nur eine einzelne Stressfaser, deren Zugrichtung nach jedem Intervall dT gemäß der lambda-Regeln aus Post [19] aktualisiert wird.
  • Die entsprechende minimale Verschiebungskomponente $\Delta L_s^{(min)}$ heißt im Programm dLMin. Wenn man sie sehr groß setzt, werden alle Verschiebungen als unbefriedigend gewertet, was zu einem Random Walk führen sollte (ungetestet).
  • Der Parameter TMX bezeichnet die Anzahl der kontinuierlichen Drift-Abschnitte (konstante Zug-Kraft), während der neue Parameter KMX die Auflösung angibt, mit dem der Runge-Kutta-Algor. die DGL innerhalb jedes Drift-Abschnittes berechnet. KMX beeinflußt nicht die Rechen-Genauigkeit.
  • Es stellt sich heraus, daß es Potentiale gibt, aus deren Minima der Agent auch mit der besten Strategie nicht entkommen kann, da der Gradient überall größer ist als seine Maximalkraft. Umgekehrt gibt es zu schwache Potentiale, über die der Agent einfach drüber-rauscht.
  • Daher setzte ich nun die Zufalls-Potentiale statt aus Gauß-Minima aus Gauß-Maxima zusammen. Die (lösbare) Aufgabe des Agenten besteht dann darin, die Hindernisse geschickt zu umschiffen.
  • Der eingestellte Testfall (orange) zeigt dieses Verhalten für dLMin=1.0. Im "Random-Walk-Fall" (hellgrün) dLMin=100.0 scheint er schlechter voranzukommen.
avoider1
  • Alles noch mit Vorbehalt. Programm könnte leicht noch Fehler enthalten.
  • Vielleicht sollten wir gar nicht mehrere Stressfasern einführen, bis dieses Modell einigermaßen verstanden ist.
  • Guter Tip: In Easyplot die Option "Connect Pts." überhaupt nicht mehr verwenden - spart Nerven.

[23] C.M.: Bemerkungen zum Lambda2-Modell

  • Im Hinblick auf Zellmigration ist das exakte Beibehalten der Zugrichtung der einen Stressfaser über viele Schritte hinweg unrealistisch. Die Stressfaser wird ja in Wahrheit spätestens dann durch eine neue ersetzt, sobald sich die Zelle an ihr bis zum fokalen Adhäsionspunkt herangezogen hat.
  • Im Lambda-Modell mit mehreren Fasern könnte man jede einzelne Faser nach einer gewissen Anzahl von Schritten getrost sterben lassen.
  • Wichtiger Hinweis: Es spart viel Zeit, wenn man in CMLib/SampledVecFunctions.h in der Methode SampledVecFu::solveDE die Zeile "printf("solveDE…\n");" ausdokumentiert !
  • Lustigerweise erinnern die Lambda-Modelle an den Quanten-Hall-Effekt: Im semiklassischen Limes bewegen sich da ja die Elektronen so wie unser Agent in einer Zufalls-Potentiallandschaft. Insbesondere gibt es dort auch die "edge currents", in welchen die Elektronen um den Systemrand herumlaufen.

[24] A.M.: Ergebnisse für Lambda2-Modell

  • Für kleine Werte von dLMin kommt der Agent (wie oben ja auch schon zu sehen) schneller voran, da er eine Richtung lange beibehält, wenn keine Hindernisse im Weg sind. Dieses Verhalten ändert sich qualitativ zwischen $dLMin=4$ und $dLMin=5$, wo er nur einen Teil der Potentiallandschaft überhaupt besucht ($TMX=100$). Ab einem Wert von $10$ ändert sich die Trajektorie bei weiterer Vergrößerung von dLMin überhaupt nicht mehr. (5dlmin.dat, random_walk.dat)
  • Bei $dLMin=3$ ist mindestens eine Kraft von $5$ nötig, damit der Agent sich schnell durch die Potentiallandschaft bewegen kann. Für kleinere Kräfte ist die Schrittweite sehr gering, was zudem noch bewirkt, dass dLMin seltener überschritten wird. Für große Kräfte wir auch die Schrittweite sehr hoch, sodass sich der Agent oft am Rand der Potentiallandschaft auf, weil er dort nicht mehr weiterkommt durch das confinement aufgehalten wird. (4force.dat)

Es wäre jetzt natürlich schön, wenn wir diese Unterschiede (anstatt durch einzelne konkrete Trajektorien) durch gemittelte Größen, wie z.B. die MSD, belegen könnten. Prinzipiell hätte ich eine MSD-Funktion für 2D-Trajektorien, aber der Systemrand führt sicher zu Artefakten. Vielleicht müssen wir doch periodische Randbedingungen einführen. Das erfordert wieder einen gewissen Programmier-Aufwand. Als einfache Alternative könnten wir als Potential eine periodische (aber mit vielen Obertönen versehene) 2D-Sinus-Reihe mit zufälligen Amplituden $A_{mn}$ und Phasen $\varphi_{mn}\in\left[ 0\ldots 2\pi \right]$ verwenden:

(22)
\begin{align} U(x,y)=\sum_{m=1,2,\ldots}\;\sum_{n=1,2,\ldots} A_{mn}\;\sin\left( (mx+ny)(2\pi/L)+\varphi_{mn} \right) \end{align}

Zur grafischen Darstellung könnte man im Nachhinein bestimmen, welche Fläche der Agent besucht hat und die Ausgabe dementsprechend dimensionieren.

  • Frage: Wie kann man in Easyplot die Landschaften quadratisch darstellen?

Man kann natürlich das Fenster von Easyplot per Augenmaß quadratisch ziehen und dann einen Screenshot davon machen. Für eine schönere Präsentation sollte man Origin oder MatLab verwenden. Dann könnte man ja mit land2D.dat die Landschaft als Konturplot oder sogar perspektivisch als Berg-und-Tal-Bild plotten und die Trajektorie drüberlegen.


[25] C.M.: Konzept 2D-Sinus Reihe

Ein erster Entwurf zur Implementierung der 2D-Sinus Reihe im Programm lambda2 findet sich hier: sinReihe.h.

Zusätzlich braucht man noch eine 2D-MSD-Routine. Eine solche wurde als Methode der DRecXY-Klasse hinzugefügt (siehe aktuelle Bibliothek CMLib).

Man könnte dann im lambda2-Programm in der Funktion ComputeTrajectory nach der Zeile traj.store("data/traj.dat"); die zusätzliche Zeile traj.MSD2D(0,"data/msd2d",0,dT); einfügen.


[26] A.M.: Re: Konzept 2D-Sinusreihe

Die 2D-Sinusreihe wurde nun implementiert (lamda3.cpp), auch die MSD-Funktion wurde eingebunden und scheint zu funktionieren.

Hey, super ! Werde ich mir gleich morgen früh mit Genuß ansehen !


[27] A.M.: Frage zu 3D-Plots

Wie kann man in Origin die Landschaft als Konturplot oder perspektivisch darstellen? Ich habe beim Einlesen der Textdatei eine Menge Spalten bekommen, mit denen ich bis jetzt keine sinnvollen Landschaften plotten konnte. Muss man die Textdatei vorher verändern und welche Funktionen von Origin nutzt man zum dann Erstellen der Plots?

In meiner alten Origin-Version 7.0 geht man folgenmaßen vor: Zuerst die 2D-Daten reinladen mit Datei/Import/ASCII. Dann die Datentabelle in der linken oberen Ecke anklicken, um sie als ausgewählt zu markieren. Dann Bearbeiten/In Matrix konvertieren/Direkt. Dann Diagramme/3D-Farbabbildung. Das sieht meistens sehr häßlich aus, da Schnellmodus mit geringer Auflösung eingestellt. Daher Format/Layer/Größe und Performance: Entwurfsmodus-Kästchen beide deaktivieren. Sicher kann man auch irgendwie die Ansichts-Perspektive ändern.


[28] A.M.: Ergebnisse für 2D-Sinusreihe

  • Bei Vergrößerung der Amplitude der Sinusschwingungen zeigt sich ein auf den ersten Blick merkwürdiger Sprung in der MSD. Zunächst sinken die Sättigungswerte von $A=1$ bis $A=3$, für $A=5$ und höhere Werte liegen die Sättigungswerte dann wieder über dem für $A=2$ und nehmen mit steigender Amplitude ab. Vielleicht liefert ein Vergleich der über die Landschaft gelegten Trajektorien eine Erklärung. (JMX=10,lPer=100,6msd_fit.dat)

Sehr spannend ! Die Kurve für amp=1 sieht ja besonders interessant aus !

Ich denke auch man muß die Trajektorien zusammen mit der Landschaft (obstacles) ansehen. Die Plateaus in der MSD sprechen ja dafür, daß der Agent in Potentialtöpfen gefangen wird. In welchen Topf er fällt, könnte von der Amplitude abhängen. Die Höhe des jeweiligen Plateaus könnte mit der Breite des jeweiligen Topfes zusammenhängen. Aber das sind einstweilen natürlich nur wilde Spekulationen.

Betrachtet man die Trajektorie zusammen mit den obstacles scheint das Teilchen immer im selben Potentialtopf gefangen zu sein. Je höher die Amplitude desto mehr ist es auch innerhalb des Potentialtopfes beschränkt. Bei $amp=1$ überschreitet es nach Passieren von zwei Potentialtöpfen die Systemgrenze und kann dort wieder frei diffundieren, was zum Exponent von $1$ am Ende der MSD-Kurve passt. (traj_amp=1.dat, traj_amp=3.dat, traj_amp=5.dat, traj_amp=10.dat)
Nach Vergleich mit der 2D-Landschaft (leider im Moment noch per Augenmaß) ist der Sprung im Sättigungswert der MSD eventuell darauf zurückzuführen, dass ab $amp=5$ das Teilchen in einem kleineren Potentialtopf im Potentialtopf gefangen ist. (2D_land.dat)

  • Bei Erhöhung von JMX bei lPer=100 und A=10 werden die Berge und Täler höher bzw. tiefer und gleichzeitig schmäler. Die Landschaften sehen für höheres JMX „wilder“ aus. (JMX=10.dat, JMX=20.dat, JMX=30.dat)

Ich bin ziemlich unzufrieden mit den 3D-Grafiken sowohl von Easyplot, als auch Origin. Vielleicht lohnt es sich doch noch, auf MatLab oder Maple umzusteigen. Ich werde mal sehen, wie einfach es ist, dort unsere 2D-Daten zu importieren.


[29] C.M. 3D-Grafik mit Maple (Version 7)

  • m:=ImportMatrix("c:\\mirror\\actual\\projects\\ANJA\\lambda\\data\\land2D.dat", format=rectangular, delimiter=" ");
  • with(plots):
  • matrixplot(m,heights=histogram,axes=boxed);
3DLands1.jpg

Irgendetwas scheint mit den Sinus-Reihen-Potentialen nicht zu stimmen. Das obige Beispiel scheint nur aus Wellen zu bestehen, die in Querrichtung $\vec{k} \propto (1,1)$ verlaufen. Das wäre ein seltsamer Zufall, wenn die Koeffizienten der Reihe gerade so ausgewürfelt wurden. Dem sollte man nachgehen …..

…. OK, ein zweiter Versuch mit höheren Frequenz-Anteilen sieht schon glaubhafter aus:

3DLands2.jpg

Jetzt müßte man nur noch die Trajektorie drüberlegen können …


[30] Influence of the Noise Spectrum on Stochastic Acceleration (K. Mallick)

Folgendes neue Paper enthält einige für unser Projekt ganz interessante Aspekte:

  • Die Gleichung (1) im Paper unterscheidet sich von unserem Problem allerdings in vielen Punkten: Es gibt einen Trägheitsterm. Das Potential ist ein einfacher, nicht-parabolischer Potentialtopf, der bei großen Auslenkungen steifer wird (umgekehrt wie beim Gauß). Das Rauschen koppelt multiplikativ an die Variable x(t).
  • Auch Mallik behandelt die Fälle von weißem Rauchen, von exponentiell-korreliertem Rauschen (Ornstein-Uhlenbeck-Prozeß), sowie von Rauschen mit Powerlaw-Korrelationen.
  • Interessant ist die Idee der "Parametrischen Resonanz": Die Energie des Teilchens nimmt besonders stark zu, wenn die im Rauschen vorkommenden Frequenzen mit charakteristischen Frequenzen des Systems übereinstimmen. So eine ähnliches Konzept steckte ja auch hinter der Idee, daß die breite Verteilung von Powerlaw-Rauschen "für alle Fälle gewappnet" ist.
  • Nützlich auch die Idee, sich die Dynamik der Energie des Teilchens anzusehen. Welche Formel würde (18) bei uns entsprechen ?
  • Fraglich ist, ob sich entweder die "Stochastic Averaging"- oder die "Carmeli-Nitzan"-Methode auf unseren Fall anwenden läßt.

[31] C.M. Nachtrag zur Transfer-Statistik

Für ein abschließendes Urteil zur Transfer-Statistik wäre es vielleicht noch interessant, aus den Daten von sigma=1,tau.dat die mittlere Transferzeit $\tau_T$ als Funktion von $\tau_c$ zu plotten, auch wenn das nur 3 Punkte ergibt. Könnte ein nicht-exponentieller, vielleicht eher linearer Zusammenhang sein. Das könnte man dann vielleicht interpretieren.

Ich habe noch ein paar mehr Punkte berechnet und in der Tat einen linearen Zusammenhang gefunden (tau_t_vs_tau_c_fit_lin)


[32] C.M. Nochmals Brainstorming

Alte Modelle:

  • Aus den Untersuchungen zur Transfer-Statistik zwischen Doppeltöpfen kann man wohl folgende Lehre ziehen: Wenn man bei gegebener Gesamt-Rauschleistung optimalen Transfer will, muß man die Korrelationszeit möglichst kurz machen (schnelles Ausprobieren neuer Kräfte hilft) und die Varianz breit (höhere Wahrscheinlichkeit, die Barrierenhöhe zu erreichen). Das ist doch eigentlich ein logisches, klares Resultat.
  • Längere Korrelationszeiten würden sich höchstens auszahlen, wenn die Möglichkeit zur konstruktiven Überlagerung verschiedener Kraftfasern besteht.

Lambda-Modelle:

  • Man könnte zwei extreme Arten von Potential-Landschaften unterscheiden: (a) Binäre Potentiale, dessen Werte entweder Null oder Unendlich sein können. Ich stelle mir z.B. eine Ebene mit U=0 vor, in die undurchdringliche, kreisförmige Inseln mit $U=\infty$ verteilt werden. (b) Glatte Potentiale, in denen der lokale Betrag der Steigung=Kraft alle kontinuierlichen Werte zwischen $-F_{mx}$ und $+F_{mx}$ annehmen kann. Also wie beim Sinus-Reihen-Potential.
  • In einer binären Potentiallandschaft würde der lambda2-Agent so lange geradeaus laufen, bis er gegen ein Hindernis läuft. Sobald dies passiert, würde er neue Kraftrichtungen ausprobieren, bis er wieder freie Bahn hat.
  • In einer glatten Potentiallandschaft ist die Verschiebung des lambda2-Agenten in jedem Zeitschritt $\Delta t_{sim}$ proportional zur Vektor-Differenz aus Eigenkraft und lokaler Potentialkraft. So lange die Verschiebungskomponente in Eigenrichtung groß genug ist, bleibt der Eigenkraftvektor unverändert. Da sich aber i.a. die Potentialkraft stetig ändert, beschreibt die Bahn des Agenten einen Bogen. Irgendwann, z.B. an einem steiler werdenden Hang, wird der Eigenkraftvektor aber doch geändert.
  • Änderungen des Eigenkraftvektors erfolgen in rein zufällige neue Richtungen. Trotz der Komplikation mit den gekrümmten Bahnen dazwischen würde man erwarten, daß auf Zeitskalen größer als die mittlere Zeit zwischen zwei Änderungen des Eigenkraftvektors bestenfalls eine diffusive Bewegung mit MSD-Exponent 1 herauskommt. Schlimmstenfalls läuft der Agent in einem Potentialtopf hin-und her, ohne je zu entkommen. Dann ergäbe sich ein Plateau in der MSD.

[33] C.M. Modulierte Diffusionskonstanten

In F. Schweitzer's Buch "Brownian Agents" scheint sich die Aktivität der Agenten darauf zu beschränken, den Reibungs-Koeffizienten in der Langevin-Gleichung zeitlich zu variieren. Erscheint mir etwas unnatürlich. Aber das führt auf ein neues Primitiv-Modell: Ein Agent könnte einfach eine weiße, gaußverteilte Rauschkraft produzieren und deren Stärke auf "intelligente" Weise modulieren. Z.B. könnte er, wenn er schon lange nicht mehr weit vom Fleck gekommen ist (gefangen in einem Topf), seine Rauschamplitude kurzzeitig erhöhen. Dies wäre ähnlich zur Chemotaxis von Bakterien. Vielleicht hilft es sogar schon, wenn er die Rauschamplitude zufällig schwanken läßt, etwa in Form einer exp-korrelierten Zufalls-Amplitude mit relativ langer Korrelationszeit.

Die stationäre PDF der Momentankraft wäre dann so etwas wie eine gaußförmig gewichtete Überlagerung von verschieden breiten Gaußkurven. Dies ergibt (gefühlsmäßig) eine nicht-gaußförmige PDF mit breiten Tails. Letztere wären nützlich zur Überwindung von Barrieren. Bei Vergleichen müßte man wieder die Gesamtleistung konstant halten.

Das wäre allerdings wieder ein neues Sub-Projekt: "Modulated Noise Model"…..


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