5  Zweidimensionale Differentialgleichungsmodelle

5.1 Räuber-Beute-Modelle

Beim Beispiel des Bakterienwachstums sind wir davon ausgegangen, dass die Abnahme der Wachstumsrate von der Bakterienbiomasse selbst abhängt. Je mehr Bakterien da sind, umso schlechter kann die Population wachsen. Dies nennt man auch intraspezifische Konkurrenz. Im folgenden wollen wir überlegen was passiert, wenn eine andere Art die Population dezimiert und deren Biomasse also von der ersten Art abhängt. Solche Modelle werden auch Räuber-Beute-Modelle genannt.

5.2 Lotka-Volterra-Modell

In diesem Abschnitt soll also das Zusammenspiel zwischen einer Beutepopulation, den Hasen, und einer Räuberpopulation, den Luchsen, simuliert werden (Abbildung 5.2). Das Modell geht auf Lotka und Volterra (Lotka 1925; Volterra 1926) zurück. Mit Hilfe dieses sehr einfachen Modells kann man die Schwankungen in der Hasen- und Luchspopulation in Nord-Kanada beschreiben (Abbildung 5.1).

Abbildung 5.1: Fellzahlen von Hasen und Luchsen, die bei der Hudson’s Bay Company von Trappern zwischen 1850 bis 1900 abgegeben wurden.

Wir machen folgende Annahmen

1. Hasen wachsen proportional zur aktuellen Hasenzahl

2. Luchse fressen einen Teil der Hasen, umso mehr, je mehr Hasen da sind.

3. Ein Teil der Luchse stirbt.

Abbildung 5.2: Diagramm des Räuber-Beute-Systems. Die Pfeile geben den Biomassenfluss an.

Sei \(B(t)\) (Beute) die Anzahl der Hasen und \(R(t)\) (Räuber) die Anzahl der Luchse zur Zeit \(t\). Wir nehmen in Kauf, dass auch nicht-ganzzahlige Werte vorkommen können. Wen das stört, der kann als Einheit auch kg Hasen und kg Luchse wählen. Das Wachstum der Hasen wird wie in Abschnitt 3.2 durch

\[\dot B = r \cdot B\]

beschrieben. Hierbei ist \(r\) die Wachstumsrate der Hasen. Die Abnahme der Räuber ist durch

\[\dot R = - s \cdot R\]

beschrieben.Hierbei ist \(s\) die Sterberate der Räuber.

Nun müssen wir noch berücksichtigen, wie viel ein Luchs frisst. Klar ist, dass er nichts zu fressen findet, wenn es keine Hasen gibt, und dass er umso mehr frisst, je mehr Hasen da sind (er muss dann nicht so lange jagen). Das Wachstum der Luchse wird also von der Hasenzahl abhängen, z.B. kann man annehmen, dass sie linear mit der Beute wächst. Die Wachstumsrate der Räuber sei hier also eine Funktion der Beute \(f(B)=b\cdot B\). Diese Rate könnte man auch als Fressrate oder Aufnahmerate bezeichnen.

\[\begin{aligned} \dot R &=& b \cdot B\cdot R \nonumber \end{aligned}\]

Alles zusammengenommen erhält man dann das Räuber-Beute-Modell von Lotka und Volterra (1926):

\[\begin{aligned} \dot B &= r \cdot B \;-\; b \cdot B \cdot R \\ \dot R &= b \cdot B\cdot R \; - \;s \cdot R \nonumber \end{aligned} \tag{5.1}\]

Simuliert man nun die Anzahl der Hasen- und Luchse (Abbildung 5.3), so fällt folgendes auf: Beide Populationen oszillieren. Zuerst steigt die Zahl der Hasen. Dadurch steht den Luchsen mehr Nahrung zur Verfügung und ihre Zahl steigt kurz darauf ebenfalls. Dies reduziert die Zahl der Hasen wieder und den Luchsen steht wieder weniger Nahrung zur Verfügung.

Abbildung 5.3: Lösungskurven des Räuber-Beute-Modells

5.2.1 Bemerkung zur Struktur des Modells

Das Modell nach Lotka und Volterra (5.1) besteht ja nun bereits aus vielen “Buchstaben”. Hier ist es nun wichtig genau zu verstehen worum es sich jeweils handelt. Auf der linken Seite der Gleichheitszeichen stehen die Änderungen oder Ableitungen (\(\dot B\) und \(\dot R\)) der Zustandsvariablen (B und R). Zustandsvariablen sind also Funktionen der Zeit (B(t) und R(t)), deren Funktionswerte gesucht sind. Da die Funktionen selbst nicht bekannt sind, werden ihre Änderungen (Ableitungen) durch das Modell beschrieben (auf der rechten Seiten der Gleichheitszeichen). Die Änderungen kann man als Funktionen der Zustandsgleichungen auffassen

\[\begin{aligned} \frac{dB(t)}{dt} &= f_1(R(t),B(t))& &\dot B= f_1(R,B)&&\\[-2ex] &&\quad kurz\quad \quad \quad &&\\[-2ex] \frac{dR(t)}{dt} &= f_2(R(t),B(t)) && \dot R= f_2(R,B)&& \end{aligned}\]

Man hat also nun die Funktionen \(f_1\) und \(f_2\) die jeweils von B und R abhängen. Sie sind Funktionen von Funktionen.

5.2.2 Bemerkung zu Parametern, Einheiten und Dimensionen

Die rechten Seiten der Gleichungen von 5.1 beschreiben also die Änderungen der Zustandsvariablen B und R. Diese Änderungen werden meistens mit Hilfe sogenannter Parameter aufgeschrieben. Als Parameter bezeichnet man in der Modellierung Größen, deren Wert man festlegt, deren Wert sich also während einer Simulation nicht ändert. Man könnte also anstelle eines Platzhalters (Buchstaben) auch den jeweiligen Wert schreiben. Dies ist jedoch umständlich, das sich die Werte durch neue Erkenntnisse aus der Literatur verändern und man diesen Wert dann an allen Stellen im Modell verändern müsste. Im Lotka-Volterra-Modell taucht z.B. der Parameter b an zwei Stellen auf.

Der Begriff Parameter bezeichnet in der Modellierung also etwas anderes als in der Biologie. Dort sind Parameter meistens Größen, die gemessen werden.

Zustandsvariablen werden häufig mit großen Buchstaben aber auch mit x oder y bezeichnet. In größeren Modellen erhalten sie aber auch länger Namen. Parameter werden in der Modellierung häufig mit kleinen Buchstaben, auch gerne aus dem griechischen Alphabet, bezeichnet. Dies ist aber nicht streng festgelegt.

Bei Modellen, deren Simulationsergebnisse mit der Realität (Messdaten) verglichen werden sollen, ist es wichtig, dass die Dimensionen und Einheiten aller Größen richtig sind, da man sonst sprichwörtlich Äpfel mit Birnen vergleicht.

In einem System drückt die Dimension einer Größe deren qualitative Eigenschaften aus. So können Räuber und Beute als Anzahl oder (Bio)-Masse angegeben werden. Werden sie als Biomasse angegeben, so muss jeweils eine Einheit z.B. Kilogramm (k) oder auch Kilogramm pro Quadratkilometer (\(kg/m^2\)) festgelegt werden. Günstigstenfalls haben Räuber und Beute dieselbe Einheit. Beim Aufstellen eines Modells ist es wichtig, sich immer über Dimension und Einheiten aller auftretender Größen im Klaren zu sein und bei jeder Gleichung die Einheiten auf Konsistenz zu überprüfen, um Fehler beim Aufstellen der Gleichungen zu vermeiden. Hierbei überprüft man, ob Terme, die aufsummiert werden, dieselbe Einheit haben und ob rechts und links vom Gleichheitszeichen dieselbe Einheit steht.

Nehmen wir hier mal an, dass wir Räuber und Beute jeweils mit der Dimension Biomasse in der Einheit Kilogramm angeben. Da beide Größen Funktionen der Zeit sind, müssen wir auch eine Zeiteinheit vorgeben. In diesen Modell sei die Einheit der Zeit Jahr (a). Auf der linken Seite stehen nun die zeitlichen Änderungen der Zustandsvariablen, also die Änderung der Biomasse mit der Zeit in der Einheit Kilogramm pro Jahr (kg/a). Auf der rechten Seite müssen alle Terme dieselbe Einheit haben, also sind \(r\cdot B\), \(r\cdot B\cdot R\) und \(s\cdot R\) auch in kg/a. daraus ergib sich sofort, dass r in 1/a angegeben werden muss, da B die Einheit kg hat. Dasselbe gilt für die Sterberate s. Schwieriger wird es bei b. Der Term \(b \cdot B \cdot R\) hat die Einheit kg/a. Da B und R die Einheit kg haben, muss b die Einheit \(1/(kg \cdot a)\) haben. Dies ist einfacher zu verstehen, wenn man sich erinnert, dass \(b\cdot B\) die Aufnahmerate war, also \(b \cdot B\) die Einheit 1/a hat.

Beispiel für konsistente Einheiten im Lotka-Volterra-Modell

\[\begin{matrix} \dot B &=& r \cdot B &-& b \cdot B \cdot R \\ [\frac{kg}{a}] &=& [\frac{1}{a}]\cdot [kg] &-& [\frac{1}{kg \cdot a}] \cdot [kg] \cdot [kg] \\ \dot R &= & b \cdot B\cdot R & -& s \cdot R \\ [\frac{kg}{a}] &=& [\frac{1}{kg \cdot a}] \cdot [kg] \cdot [kg]& -& [ \frac{1}{a}]\cdot [kg] \end{matrix}\]

Wir haben im Modell nun also drei Raten, die Wachstumsrate der Beute r, die Sterberate der Räuber s und die Aufnahme- oder Fressrate der Räuber \(b\cdot B\). Alle drei haben die Einheit 1/a. Man erkennt, dass b alleine keine Rate darstellt. Die Wachstumsrate r wird in der Ökologie auch per capita growth rate (Pro-Kopf-Wachstumsrate) genannt.

5.2.3 Fixpunkte

Ähnlich wie beim eindimensionalen Modell kann man die stationären Zustände des Modells bestimmen. Hierzu werden beide Gleichungen gleich Null gesetzt. Man erhält ein Gleichungssystem aus zwei Gleichungen mit den zwei Unbekannten B und R.

\[\begin{aligned} \dot B &= r \cdot B \;-\; b \cdot B \cdot R = (r \;-\; b \cdot R) \cdot B=0 \\ \dot R &= b \cdot B\cdot R \; - \;s \cdot R = (b \cdot B \; - \;s) \cdot R=0\nonumber \end{aligned} \tag{5.2}\]

Man braucht es nicht als Gleichungssystem zu lösen, sondern kann die Lösungen auch durch scharfes Hinsehen finden: Die erste Gleichung ist null, wenn entweder

\[B=0 \quad oder \quad r -b \cdot R = 0\]

gilt. Die zweite Gleichung ist null, wenn entweder

\[R=0 \quad oder \quad b \cdot B -s = 0\]

gilt. Damit gibt es 4 mögliche Kombinationen:

1. B=0 und R=0

2. B=0 und \(b\cdot B -s = 0\)

3. \(r-b \cdot R = 0\) und R=0

4. \(r-b\cdot R = 0\) und \(b\cdot B -s = 0\)

Fall 1 führt zur trivialen Lösung, dass B und R Null sind. Ökologisch bedeutet das, das weder Räuber noch Beute da sind. Wir haben also den Fixpunkt

\[(B_1^*,R_2^*) = (0,0)\]

Der Fixpunkt wird wiederum durch \(\ast\) gekennzeichnet. Betrachtet man den Fall 2 so bedeutet dies das B=0 und \(B=\frac{s}{b}\) gelten soll. Das ist aber nicht möglich, da wir \(r,s > 0\) angenommen haben. Hieraus ergibt sich also kein Fixpunkt. Dasselbe gilt für Fall 3. Einen weiteren Fixpunkt erhalten wir, wenn wir in Fall 4 nach R bzw. B auflösen. Der zweite Fixpunkt lautet dann

\[(B_2^*,R_2^*) = (\frac{s}{b},\frac{r}{b})\]

Wenn das System genau diesen Zustand hat, wird sich nichts mehr ändern. Sowohl Beute als auch Räuber haben eine konstante Populationsgröße. Eine kurze Zusammenfassung zur Analysis mehrdimensionaler Systeme und zur Stabilitätsanalyse zweidimensionaler Systeme ist in Anhang G *und I *beschrieben.

5.2.4 Richtungsfeld

Für ein autonomes 2-dimensionales System kann man ebenfalls ein Richtungsfeld zeichnen. Man betrachtet hierzu

\[\frac{dR}{dB}= \frac{ \frac{dR}{dt}}{\frac{dB}{dt}}\]

d.h. die Änderung des Räubers nach der Beute. Im Falle des Lotka-Volterra-Modells gilt

\[\frac{dR}{dB}= \frac{ b \cdot B\cdot R \; - \;s \cdot R}{ r \cdot B \;-\; b \cdot B \cdot R }\]

Abbildung 5.4: Richtungsfeld des Lotka-Volterra-Modells mit Lösungskurven zu verschiedenen Anfangswerten.

Die R-B-Ebene nennt man auch Zustandsraum oder Phasenraum des Systems, eine Lösungskurve auch Bahn oder Trajektorie¹. Beides zusammen heißt auch Phasendiagramm.

5.2.5 Teilgleichgewichte (Halbgleichgewichte, Isoklinen)

Im Richtungsfeld kann man zwei ausgezeichnete Richtungen erkennen. Die, wo die Richtungslinien, waagerecht verlaufen und die, wo die Richtungslinien senkrecht verlaufen. Diese ergeben sich, wenn man die Gleichungen aus 5.2 einzeln betrachtet. Die Gleichung für die Beute

\[\dot B = r \cdot B \;-\; b \cdot B \cdot R = (r \;-\; b \cdot R) \cdot B=0\nonumber\]

gilt für \(B=0\) und für \(R=\frac{r}{b}\). Für diese Werte gilt, dass sich \(B\) nicht ändert. Die Richtungslinien laufen senkrecht. Das Teilgleichgewicht von B kann man also als die Menge aller Punkte im Phasenraum, für die B keine Änderung erfährt, definieren:

\[\mbox{Teilgleichgewicht von B} = \lbrace \dot B=0\rbrace= \lbrace(B,R) \vert (r \;-\; b \cdot R) \cdot B=0\rbrace\]

Analog kann man das Teilgleichgewicht für die Räuber definieren:

\[\dot R = b \cdot B\cdot R \; - \;s \cdot R = (b \cdot B \; - \;s) \cdot R=0\nonumber\]

gilt für \(R=0\) und \(B=\frac{s}{b}\). Für diese Werte gilt, dass sich \(R\) nicht ändert. Die Richtungslinien laufen waagerecht.

\[\mbox{Teilgleichgewicht von R} = \lbrace \dot R=0 \rbrace= \lbrace(B,R) \vert \; (b \cdot B \; - \;s) \cdot R=0\rbrace\]

Die Teilgleichgwichte sind in Abbildung 5.5 dargestellt. Die Schnittpunkte der Teilgleichgewichte sind gerade die Fixpunkte. Teilgleichgewichte werden auch Halbgleichgewichte oder Isoklinen genannt.

Abbildung 5.5: Teilgleichgewichte für das Lotka-Volterra-Modell. Die blauen Linien zeigen das Teilgleichgewicht \(\lbrace \dot R=0\rbrace\), die grünen Linien das Teilgleichgewicht \(\lbrace \dot B=0\rbrace\).

5.3 Rosenzweig-MacArthur-Modell

Das bisherige Modell kann natürlich verbessert werden. Folgende Verbesserungen fallen einem sofort ein:

• Die Hasen dürfen nicht unbegrenzt wachsen, wenn keine Luchse da sind. Es muss eine Kapazität wie in Abschnitt 3.2 geben.

• Die Luchse müssen satt werden. Die Zunahme der Wachstumsrate der Luchse darf nicht linear von der Zahl der Hasen abhängen.

• Nur ein Teil der Hasen kann als Nahrung verwertet werden (Nahrungsverwertungseffizienz, “Fell und Knochen werden ausgeschieden”).

5.3.1 Kapazität

Verbessern wir das Lotka-Volterra-Modell aus 5.2 nun dahingehend, dass wir eine Kapazität für Hasen berücksichtigen. In diesem Fall wird das Wachstum der Hasen durch Gleichung 4.3 beschrieben, und das vollständige Modell lautet:

\[\begin{aligned} \dot B&=& \; r \cdot B \cdot \left (1- \frac{B}{K} \right )\;-\; b \cdot B \cdot R \\[1ex] \dot R &=& \; b \cdot B \cdot R \; - \;s \cdot R \nonumber \end{aligned}\]

Die Modellergebnisse sind in Abbildung 5.4 für \(K=250\) und \(K=500\) dargestellt (andere Parameter wie bisher). Die Oszillationen sind nun mehr oder weniger stark gedämpft. Die Zahl der Hasen und Luchse strebt letztlich auf konstante Werte zu.

 

 

Abbildung 5.6: Lösungskurven des Räuber-Beute Systems mit einer Kapazität für die Hasen (Beute) von K = 250 (links) und K = 500 (rechts).

5.3.2 Sättigung

Will man nun zusätzlich berücksichtigen, dass die Luchse satt werden, so muss man anstelle der monoton steigenden Fressrate \(b\cdot B\) eine Fressrate annehmen, die durch einen maximalen Wert begrenzt ist. Diese maximale Fressrate gibt an wie viel ein Luchs pro Zeiteinheit maximal fressen kann. Dies ist u.a. von der so genannten “handling time” abhängig, der Zeitspanne, die ein Luchs benötigt, die Nahrung aufzunehmen (den Hasen zu zerlegen). Ein weiterer Faktor ist die Stoffwechselgeschwindigkeit, d.h. wie schnell Biomasse aus der Nahrung aufgebaut werden kann. Die tatsächliche Fressrate soll vom Nahrungsangebot abhängen. Hierbei kann man davon ausgehen, dass die Rate anfangs mehr oder weniger linear ansteigt und sich dann dem Maximalwert nähert.

Hierzu geeignet ist die Monod-Funktion (siehe Anhang A.7):

\[f(B)= b \cdot \frac{B}{B+M}\]

Die Funktion f geht für \(B\rightarrow \infty\) gegen b. Somit ist \(b\) die maximale Wachstumsrate. Der Wert M ist die sogenannte Halbsättigungskonstante. Setzt man M in die Funktion f ein, so gilt \(f(M)=\frac{b}{2}\). Es wird als gerade die Hälfte der maximalen Rate b erreicht. Die Wachstumsrate b ist in der Dimension 1/Zeit und hat z.B. die Einheit 1/Jahr, M hat dieselbe Dimension und Einheit wie B, also zum Beispiel Biomasse in kg.

Abbildung 5.7: Graph der Monod-Funktion mit der Halbsättigungskonstanten \(M\) und dem Grenzwert \(b\).

Der Graph der Monod-Funktion ist in Abbildung 5.8 für große und kleine Werte von \(M\) gegeben. Die Steigung bei \(B=0\) beträgt \(\frac{1}{M}\).

Abbildung 5.8: Graph der Monod-Funktion mit großer (rot) und keiner (schwarz) Halbsättigungskonstanten \(M\) und \(b=1\).

Für sehr kleine Halbsättigungskonstanten M ist die Fressrate f(B) für nahezu alle Beutekonzentrationen konstant. Der Räuber reagiert quasi nicht auf die Beutekonzentration. Beute und Räuber schwingen sehr stark. Für große Halbsättigungskonstanten reagiert der Räuber auf die Beutekonzentration, das System wird stabilisiert (Abbildung 5.9).
Für sehr große Werte von \(M\) (\(M>>K\)) wird die Fressrate nahezu null, der Räuber stirbt aus.

 

 

Abbildung 5.9: Simulationsergebnisse des Modells mit der Formulierung in Gleichung A.11 für M = 150 und M = 200. Weitere Parameterwerte: r = 0.25, K = 250, n = 0.7 s = 0.05, b = 0.3

5.3.3 Effizienz

Um zu berücksichtigen, dass nicht die gesamte durch den Räuber aufgenommene Beutebiomasse zu eigener Biomasse wird, wird ein dimensionsloser Parameter \(n\) eingeführt. Er kann Werte zwischen Null (keine Verwertung) und Eins (100%ige Verwertung) annehmen. Letzteres ist praktisch allerdings unmöglich. Dies geht auf ein physikalisches Grundprinzip zurück, den 2. Hauptsatz der Thermodynamik.

5.3.4 Modellgleichungen und Simulation

Man erhält nun folgendes allgemeinere Räuber-Beute-Modell ähnlich dem Modell von Rosenzweig und MacArthur (Rosenzweig und MacArthur 1963) mit den Zustandsvariablen B und R und den Parametern gemäß Tabelle 5.1:

\[\begin{aligned} \dot B&=& \; r \cdot B \cdot \left (1- \frac{B}{K} \right )\;-\; b\cdot \frac{B}{B+M} \cdot R \\[1ex] \dot R &=& n \cdot b\cdot \frac{B}{B+M} \cdot R \; - \;s \cdot R \nonumber \end{aligned} \tag{5.3}\]

r	Wachstumsrate der Beute
K	Kapazität der Beute
b	maximale Aufnahmerate des Räubers
M	Halbsättigungskonstante des Räubers
n	Nahrungsaufnahmeeffizienz des Räubers
s	Sterberate des Räubers

Tabelle 5.1: Parameter des Rosenzweig-MacArthur-Modells

Eigentlich gibt es das Rosenzweig-MacArthur-Modell nicht. Es handelt sich vielmehr um eine Klasse von Modellen mit ähnlicher Struktur. Hier wird eine typische Realisierung behandelt.

Die Simulationsergebnisse sind in Abbildung 5.10 gegeben. Man sieht, dass das System bei bestimmten Parameterkonstellationen einen sogenannten Grenzkreis erreicht. Ein Grenzkreis stellt selbst eine geschlossen Lösung des Systems dar. Lösungen, die sich dem Grenzkreis nähern, erreichen diesen aber nicht erreichen, ihm aber beliebig nahe kommen.

 

 

 

 

 

 

Abbildung 5.10: Simulationsergebnisse für das allgemeine Räuber-Beute-Modell nach -@ eq-rballg mit den Parametern r = 0.01, K = 250, M = 30, s = 0.05, n = 0.7, b = 0.1 (oben), r = 0.01, K = 250, M = 50, s = 0.05, n = 0.7, b = 0.1 (mitte) und r = 0.01, K = 250, M = 30, s = 0.07, n = 0.7, b = 0.1 (unten). Links ist jeweils das Zeitdiagramm (Beute:schwarz, Räuber:rot), rechts das Phasendiagramm dargestellt.

Lotka, A. J. 1925. Elements of Physical Biology. Williams; Wilkins.

Rosenzweig, M. L., und R. H. MacArthur. 1963. „Graphical representation and stability conditions of predator-prey interactions“. American Naturalist 97 (895): 209–23. https://doi.org/10.1086/282272.

Volterra, V. 1926. „Variazioni e fluttuazioni del numero d’individui in specie animali conviventi“. Mem. Acad. Lincei Roma 2: 31–113.

---

1. Diese Bezeichnung ist etwas ungenau. Streng genommen ist die Trajektorie eine Abbildung der Zeit in den Zustandsraum, eine Lösung des Anfangswertproblem. Die Bahn ist das Bild, enthält also die Zeitinformation nicht mehr.