Nachdem es hier in den letzten Wochen recht ruhig war, moechte ich jetzt eine kleine Serie von Beitraegen starten, wo ich einen Aspekt meiner Arbeit beschreiben will, naemlich stochastische Modelle und was man mit ihnen ueber biologische Systeme lernen kann. Als allererstes, moechte ich die generellen Unterschiede zwischen einem deterministischen System und einem stochastischen System beschreiben. Deterministisch (von lateinisch determinare - ‘abgrenzen’, ‘festlegen’) bedeutet im naturwissenschaftlichen Sinne, dass ein System, wenn der Zustand zu einem beliebigen Zeitpunkt bekannt ist, man den Zustand fuer spaetere Zeitpunkte exakt berechnen kann. Einfachstes Beispiel aus der Physik ist wohl der beruehmte Apfel, der vom Baum faellt. Wenn ich den Ort und die Geschwindigkeit des Apfels zu Anfang seines Falles kenne, kann ich mittels Newtons Gleichung ausrechnen, wann der Apfel am Boden aufschlaegt und wie schnell er dabei ist (freier Fall).
Eine Anwendung aus der Biologie waere zum Beispiel ein Gen, das zu einem bestimmten Zeitpunkt eingeschaltet wird. Das Produkt Y (mRNA) wird nun mit einer festen Rate (Produkt pro Zeit) produziert und zerfaellt mit einer Rate
. Dabei wird der Zerfall des Produktes natuerlich von der aktuellen Menge von Produkt abhaengen, d.h. wenn mehr Produkt da ist, wird auch mehr zerfallen. Zu Anfang ist kein Produkt vorhanden.
Ein solcher Prozess wird sich dabei in ein Gleichgewicht bewegen in dem Produktion und Zerfall sich die Waage halten. Mathematisch wird das durch eine einfache Differentialgleichung beschrieben:

Dabei ist die erste zeitliche Ableitung, also die zeitliche Aenderung der Produktion. Im Gleichgewicht sind Produktion und Zerfall gleich gross, also die zeitliche Ableitung Null. Daraus folgt der Gleichgewichtswert
. Die Loesung der Differentialgleichung steht in der dritten Zeile, sie stellt den eigentlichen zeitlichen Verlauf von Y dar. Man sieht, dass der Exponentialterm fuer grosse t verschwindet und Y mit der Zeit gegen seinen Gleichgewichtswert strebt. Wie schnell oder langsam das geht, haengt einzig vom Zerfallsparameter
ab mit der Halbwertszeit
.
Nun besteht ja die Natur aus diskreten Bausteinen und chemische Reaktionen laufen in diskreten Schritten ab. Nur wenn sehr viele davon miteinander reagieren, bekommen wir davon nix mit und es erscheint uns als ein kontinuierlicher Prozess. In einer lebenden Zelle sieht das mitunter schon wieder anders aus, da es von vielen Proteinen nur wenige hundert oder tausend Kopien gibt und die Produktion dieser Proteine ist fuer die Zelle deshalb auch ein diskreter Process, sie muss erst das Gen transkribieren, die mRNA zurechtschneiden und dann in ein Protein uebersetzen und das dann eventuell noch zurechtfalten oder sonstwie modifizieren.
Wie laesst sich nun ein solcher Prozess stochastisch beschreiben? An welcher Stelle tritt eigentlich die Zufaelligkeit zu Tage? Wenn wir uns wieder unser Gen von eben hernehmen und wir wissen, dass es aktiv ist, was kann als naechstes passieren? Es kann ein einzelnes Produkt produziert werden (das Gen also transkribiert in ein mRNA Molekuel), allerdings koennen wir nicht mit Sicherheit sagen wann. Zerfallen kann erstmal noch nix. Wenn wir nun das erste Produkt Molekuel haben, koennen zwei Dinge passieren, es kann wieder zerfallen oder es wird noch ein weiteres produziert. Es gibt also in diesem Zusammenhang zwei zufaellige Ereignisse, naemlich wann eine Reaktion passiert und welche (Produktion oder Zerfall). Die dazugehoerigen Wahrscheinlichkeiten werden dabei von den Raten und
bestimmt. Solche Prozesse heissen Birth-Death-Prozesse. Dabei folgt die Anzahl von Reaktionen, die in einem Intervall
stattfinden einem Poisson Prozess, d.h. die Laenge der Zeitintervalle zwischen zwei Reaktionen ist Exponentiell verteilt.
Um das ganze nun zu simulieren braucht man einen Zufallszahlen Generator, allerdings gibt es da sehr gute und entsprechende Funktionen bzw. Methoden in allen gaengigen Programmiersprachen. Desweiteren braucht man natuerlich die beiden Raten und einen Startwert. Als erstes generiert man sich nun also ein Zeitinterval als Sample einer exponentiell verteilten Zufallsvariable, der Parameter haengt dabei natuerlich von den beiden Raten ab und der aktuellen Anzahl von Molekuelen:
. Wenn man das Zeitintervall hat, braucht man eine zweite Zufallszahl um festzulegen, welche der beiden moeglichen Reaktionen Produktion und Zerfall stattfindet. Das kann man z.B. durch eine gleichverteilte Zufallszahl aus dem Intervall [0,1] machen und wenn die Zahl kleiner als
ist, wird ein weiteres Molekuel produziert. So ist die Wahrscheinlichkeit fuer einen Zerfall gering, so lange es wenige Produktmolekuele y gibt. Andersrum wird ein Zerfall umso wahrscheinlicher, je mehr produziert wurde. Wenn man weiss welche Reaktion stattfindet, passt man die Anzahl an Produktmolekuelen y entsprechend an und mach weiter mit dem naechsten Zeitintervall. Und so weiter…
In einem solchen Modell haengt die Zahl von Reaktionen pro Zeit und damit die Annaeherung an ein Gleichgewicht nicht mehr nur vom Zerfallsparameter ab, sondern auch von der Produktionsrate
. Das sieht man recht gut in den Plots fuer zwei Simulationen solcher Systeme. Im ersten Fall ist
und
und damit der Gleichgewichtswert 10. Die farbigen Linien zeigen einzelne Realisationen fuer ein solches System und die dicke schwarze Linie den deterministischen Fall mit identischen Paramtern. Darunter ist ein Histogram der Anzahl der mRNA Molekuele nachdem das Gleichgewicht erreicht wurde und als drittes folgt ein Histogramm der Laengen der Zeitintervalle zwischen zwei aufeinanderfolgenden Reaktionen. Letztere sind wie beschrieben exponentiell verteilt.

Im zweiten Fall wurden beide Reaktionsraten um den Faktor 10 verringert, wobei natuerlich ihr Verhaeltnis und damit der Gleichgewichtswert bei 10 blieb. Das Histogramm der Molekuelzahlen sieht auch genauso aus, wie im ersten Beispiel. Nur sind die Zeitintervalle zwischen den Reaktionen deutlich laenger und damit auch die Fluktuationen viel langsamer.
Der hier vorgestellte Birth-Death-Prozess ist ein Spezialfall eines Markov-Prozesses mit stetiger Zeit, d.h. die Reaktionen finden nicht in festen Zeitintervallen statt, sondern die Zeitintervalle sind variabel. Bei solchen Prozessen haengt der zukuenftige Systemzustand nur vom aktuellen Zustand ab und nicht von Zustaenden weiter in der Vergangenheit.
Im naechsten Post will ich einen Algorithmus zur Simulation von stochastischen Prozessen genauer Vorstellen, der von Daniel Gillespie in den 70er Jahren entwickelt wurde. Dazu gibt es dann sicher wieder ein oder zwei Beispiele.
Zum Abschluss wiedermal eine Buchempfehlung: vor einigen Wochen bin ich ueber Darren Wilkinsons Buch Stochastic modelling for systems biology gestolpert und war spontan begeistert. D. Wilkinson ist Professor fuer Statistik an der Uni in Newcastle und schreibt nicht nur Buecher zur stochastischen Modelierung, sondern auch Software. Zum Buch gibt es ein R Paket, mit allen Beispielen. Damit laesst sich schoen rumspielen und man kann direkt auch eigene einfache Modelle basteln.











