Zagadnienia

11.1 Równania ewolucyjne: równanie Fishera – Kołmogorowa
- Równanie Skellama
- Równanie Fishera – Kołmogorowa
11.2 Wzory Turinga

11. Modele z zależnością przestrzenną: dyfuzja w procesach biologicznych

Do tej pory rozważaliśmy modele populacyjne, w których nie uwzględnialiśmy zależności od położenia osobników w przestrzeni, a tylko średnie zagęszczenia i ich zmiany w czasie. Takie podejście do modelowania dynamiki populacji kryje za sobą dodatkowe założenia dotyczące równomiernego rozłożenia przestrzennego osobników danej populacji oraz braku istotnego wpływu zmian położenia na opisywany proces. Oczywiste wydaje się, że w wielu przypadkach położenie i ruch w przestrzeni odgrywają dużą, jeśli nie zasadniczą rolę w przebiegu danego procesu. W związku z tym kolejną część wykładu poświęcimy omówieniu jednego z możliwych sposobów modelowania ruchu, tzw. procesu dyfuzji. W przypadku procesów chemicznych dyfuzja opisuje ruch substancji w kierunku od większych do mniejszych stężeń. Możliwe są oczywiście także przepływy innego typu — np. chemotaksja działa w pewnym sensie przeciwnie do dyfuzji.

Wyprowadzimy najpierw równanie dyfuzji substancji chemicznej w przestrzeni jednowymiarowej (dla przestrzeni dwu i trójwymiarowej rozumowanie przebiega analogicznie). Załóżmy, że mamy bardzo wąską rurkę o długości $L$ i przekroju $S$ , wzdłuż której przemieszcza się pewna substancja $\mathcal{A}$ . Rurka jest na tyle wąska, że ruch cząsteczek odbywa się tylko wzdłuż niej. Niech stężenie cząsteczek tej substancji w chwili $t$ w miejscu $x$ wynosi $A(t,x)$ i załóżmy, że $A$ jest dostatecznie gładką funkcją obu zmiennych.

Z punktem $x$ związujemy wycinek rurki $(x,x+\Delta x)$ o objętości $\Delta V=S\Delta x$ . Zmiana stężenia substancji $\mathcal{A}$ w tym wycinku następuje wyniku przepływu — mamy z jednej strony napływ do wycinka, z drugiej odpływ z niego. Zgodnie z prawem Ficka masa substancji przepływająca przez wąski przekrój w ciągu krótkiego przedziału czasu $(t,t+\Delta t)$ jest proporcjonalna do gradientu stężenia (w naszym jednowymiarowym przypadku do $\tfrac{\partial A}{\partial x}$ ), długości przedziału czasu i pola przekroju rurki, zatem

$Q_{x}=DS\frac{\partial A(t,x)}{\partial x}\Delta t$

w punkcie $x$ oraz

$Q_{{x+\Delta x}}=DS\frac{\partial A(t,x+\Delta x)}{\partial x}\Delta t$

w punkcie $\Delta x$ , gdzie $D>0$ jest pewnym współczynnikiem proporcjonalności. Zatem zmiana masy w objętości $\Delta V$ wynosi $\Delta Q=Q_{{x+\Delta x}}-Q_{x}$ , czyli zmiana stężenia substancji w wycinku wyraża się wzorem

$\Delta A=\frac{\Delta Q}{\Delta V}=\frac{\left(\frac{\partial A(t,x+\Delta x)}{\partial x}-\frac{\partial A(t,x)}{\partial x}\right)DS\Delta t}{S\Delta x}=D\frac{\frac{\partial A(t,x+\Delta x)}{\partial x}-\frac{\partial A(t,x)}{\partial x}}{\Delta x}\Delta t.$

Stąd, dzieląc stronami powyższą równość przez $\Delta t$ i przechodząc do granicy $\Delta x\to 0$ oraz $\Delta t\to 0$ dostajemy

$\frac{\partial A}{\partial t}=D\frac{\partial^{2}A}{\partial x^{2}},$

(11.1)

czyli jednowymiarowe równanie dyfuzji. W ogólnym przypadku równanie dyfuzji ma postać

$\frac{\partial A}{\partial t}=D\triangle A,$

gdzie

$\triangle A=\sum\limits _{{i=1}}^{n}\frac{\partial^{2}A}{\partial x_{i}^{2}},$

$n$ oznacza wymiar przestrzeni (w przypadku procesów obserwowanych w naturze mamy $n=1,2$ lub $3$ ) i $\triangle$ nazywamy operatorem Laplace'a lub krócej laplasjanem, a $D>0$ — współczynnikiem dyfuzji.

Z naszego populacyjnego punktu widzenia będziemy rozpatrywać dyfuzję jako proces losowego przemieszczania się osobników w przestrzeni. Najprostsze podejście do tego typu zagadnień reprezentowane jest przez tzw. ruchy Browna. Ponownie dla uproszczenia omówimy sytuację jednowymiarową. Rozważamy prostą podzieloną na odcinki o długości $\Delta x$ . Niech osobnik znajduje się w chwili $t$ w odcinku $[x,x+\Delta x)$ i w przedziale czasu o długości $\Delta t$ może się przemieścić do odcinka na lewo, na prawo, albo pozostać w tym samym odcinku. Załóżmy, że prawdopodobieństwo przemieszczenia się do odcinka $[x-\Delta x,x)$ wynosi $p_{l}$ , a do odcinka $[x+\Delta x,x+2\Delta x)$ — $p_{p}$ . Niech $a(t,x)$ będzie liczbą osobników w chwili $t$ w odcinku $[x,x+\Delta x)$ . Po upływie czasu $\Delta t$ w tym samym odcinku mamy

$a(t+\Delta t,x)=a(t,x)+p_{p}a(t,x-\Delta x)-p_{p}a(t,x)+p_{l}a(t,x+\Delta x)-p_{l}a(t,x),$

(11.2)

gdzie $p_{p}a(t,x-\Delta x)$ oznacza średnią liczbę osobników, które w czasie $t$ były na lewo od rozpatrywanego odcinka i przemieściły do niego, $p_{l}a(t,x+\Delta x)$ to osobniki, które były na prawo i przeszły na lewo, natomiast $p_{p}a(t,x)$ , $p_{l}a(t,x)$ odpowiadają tym osobnikom, które były w omawianym odcinku w czasie $t$ i w przedziale czasu $(t,t+\Delta t)$ przemieściły się odpowiednio na prawo lub lewo. Załóżmy teraz, że funkcję $a$ możemy rozwinąć w szereg Taylora w punkcie $(t,x)$ zarówno jako funkcję zmiennej $t$ z parametrem $x$ , jak i zmiennej $x$ z parametrem $t$ . Otrzymujemy odpowiednio

	$\displaystyle a(t+\Delta t,x)$	$\displaystyle=a(t,x)+\frac{\partial a(t,x)}{\partial t}\Delta t+\frac{1}{2}\frac{\partial^{2}a(t,x)}{\partial t^{2}}(\Delta t)^{2}+\ldots$
	$\displaystyle a(t,x+\Delta x)$	$\displaystyle=a(t,x)+\frac{\partial a(t,x)}{\partial x}\Delta x+\frac{1}{2}\frac{\partial^{2}a(t,x)}{\partial x^{2}}(\Delta x)^{2}+\ldots$

Wstawiając powyższe rozwinięcia do równania (11.2) dostajemy

$\begin{split}\frac{\partial a(t,x)}{\partial t}\Delta t+\frac{1}{2}\frac{\partial^{2}a(t,x)}{\partial t^{2}}(\Delta t)^{2}&=p_{p}\left(-\frac{\partial a(t,x)}{\partial x}\Delta x+\frac{1}{2}\frac{\partial^{2}a(t,x)}{\partial x^{2}}(\Delta x)^{2}\right)+\\ &+p_{l}\left(\frac{\partial a(t,x)}{\partial x}\Delta x+\frac{1}{2}\frac{\partial^{2}a(t,x)}{\partial x^{2}}(\Delta x)^{2}\right)+\ldots\\ \end{split}$

Załóżmy teraz, że $p_{l}=p_{p}=\tfrac{1}{2}$ , czyli osobnik z jednakowym prawdopodobieństwem przemieszcza się zarówno w lewo, jak i w prawo. Mamy wtedy

$\frac{\partial a(t,x)}{\partial t}\Delta t+\frac{1}{2}\frac{\partial^{2}a(t,x)}{\partial t^{2}}(\Delta t)^{2}=\frac{1}{2}\frac{\partial^{2}a(t,x)}{\partial x^{2}}(\Delta x)^{2}+\ldots$

Dzieląc to równanie przez $\Delta t$ , przechodząc do granicy $\Delta t\to 0$ i $\Delta x\to 0$ w taki sposób, że $\tfrac{(\Delta x)^{2}}{2\Delta t}=D$ dostajemy równanie (11.1).

Otrzymaliśmy więc związek losowego przemieszczenia się osobników z równaniem dyfuzji z jednej strony i pewnym łańcuchem Markowa z drugiej.

11.1. Równania ewolucyjne: równanie Fishera – Kołmogorowa

W przypadku rzeczywistych populacji oprócz przemieszczania się rozpatrujemy zwykle także inne procesy, np. rozrodczości/śmiertelności. W ogólnym przypadku otrzymujemy równanie reakcji-dyfuzji

$\frac{\partial N(t,x)}{\partial t}=f(t,N(t,x))+\triangle N(t,x)\,,$

(11.3)

zwane także równaniem ewolucyjnym, które opisuje zmiany zagęszczenia populacji w chwili $t$ w miejscu o położeniu $x$ w przestrzeni. Czasami także w odniesieniu do funkcji $f$ używa się nazwy kinetyka oddziaływań (albo reakcji w kontekście biochemicznym), a odpowiadające równaniu (11.3) równanie zwyczajne $\dot{N}=f(t,N)$ zwane jest wtedy równaniem kinetycznym.

Równanie Skellama

Najprostsze tego typu równanie z liniową kinetyką zostało użyte przez Skellama w 1951 roku do opisu rozprzestrzeniania się populacji. Dokładniej, za pomocą równania (11.3) z funkcją $f(t,N)=\alpha N$ opisał on rozprzestrzenianie się populacji piżmaków w Europie. Gatunek ten został sprowadzony do Europy w 1905 roku przez węgierskiego posiadacza ziemskiego i szybko rozprzestrzenił się wokół. Mimo że obszar Europy jest ograniczony, Skellam rozpatrywał to równanie w całym $\mathbb{R}^{2}$ z warunkiem początkowym, który odpowiada introdukcji kilku osobników w punkcie $(0,0)$ , czyli $N(0,x,y)=M\delta(x,y)$ , gdzie $M$ to początkowa liczba osobników, a $\delta(x,y)=\delta(x)\delta(y)$ , gdzie $\delta(x)$ oznacza deltę-Diraca, czyli $\delta(x)=0$ dla $x\ne 0$ oraz $\int\limits _{{-\infty}}^{{\infty}}\delta(x)dx=1$ . Aby prześledzić rozprzestrzenianie się osobników zgodnie z modelem Skellama zbadajmy najpierw zagadnienie prostsze.

Rozpatrzmy równanie dyfuzji (11.1) w całym $\mathbb{R}$ z warunkiem początkowym $A(0,x)=\delta(x)$ , gdzie $\delta$ oznacza jednowymiarową deltę-Diraca. Można pokazać, że rozwiązaniem tego zagadnienia początkowego jest funkcja

$A(t,x)=\frac{1}{\sqrt{4\pi Dt}}\exp\left(-\frac{x^{2}}{4Dt}\right),$

którą nazywamy rozwiązaniem fundamentalnym. Oczywiście spełnienie warunku początkowego należy rozumieć w sensie granicznym. Ogólniej, jeśli rozważamy dyfuzję w przestrzeni $n$ -wymiarowej z warunkiem początkowym $\delta(\mathbf{x})=\Pi\delta(x_{i})$ , $\mathbf{x}=(x_{1},x_{2},\ldots,x_{n})$ , to rozwiązanie fundamentalne ma postać

$A(t,\mathbf{x})=\frac{1}{(4\pi Dt)^{{n/2}}}\exp\left(-\frac{|\mathbf{x}|^{2}}{4Dt}\right).$

Teraz połączymy znajomość rozwiązania fundamentalnego równania dyfuzji dla $n=2$ i rozwiązania równania kinetycznego $\dot{N}=\alpha N$ w celu otrzymania rozwiązania równania Skellama

$N(t,r)=\frac{M}{4\pi Dt}\exp\left(\alpha t-\frac{r^{2}}{4Dt}\right),\quad r=\sqrt{x^{2}+y^{2}}\,,$

które zadaje falę inwazyjną. W punktach o ustalonym promieniu $r$ w pewnej chwili $t$ pojawiają się osobniki, których tam wcześniej w zasadzie nie było (należy to interpretować w taki sposób, że fizycznie tych osobników nie było, mamy więc do czynienia z przełączeniem, przy czym funkcja $A$ w sposób ciągły przybliża taką dynamikę), a choć ich liczebność wraz z upływem czasu maleje, ponieważ stale rozprzestrzeniają się dalej, to populacja jako taka w tym miejscu egzystuje, por rys. 11.1.

Rys. 11.1. Przykładowe rozwiązanie równania Skellama.

Równanie Fishera – Kołmogorowa

Równanie Fishera [3] lub inaczej Fishera – Kołmogorowa jest z kolei najprostszym nieliniowym równaniem typu ewolucyjnego i stanowi naturalne uogólnienie równania logistycznego na przypadek populacji złożonej z osobników, które mogą się przemieszczać i ruch ten ma znaczenie dla opisywanego procesu. Jest to także prototypowe równanie, w którym występują rozwiązania w postaci fal biegnących. Typowo jest to równanie rozpatrywane w całej przestrzeni $\mathbb{R}^{n}$ , gdzie $n$ zwykle odpowiada naturalnym wymiarom $1$ , $2$ lub $3$ . Zostało zaproponowane przez Fishera w 1937 roku do opisu rozprzestrzeniania się genu dominującego w populacji.

Równanie to ma postać

$\frac{\partial N(t,x)}{\partial t}=D\triangle N(t,x)+rN(t,x)\left(1-{N(t,x)}\right),$

(11.4)

przy czym w najprostszym przypadku, który omówimy, $x\in\mathbb{R}$ , $\triangle N=\tfrac{\partial^{2}N}{\partial x^{2}}$ , a $D>0$ jest współczynnikiem dyfuzji.

Pokażemy, że równanie (11.4) ma rozwiązania w postaci fal biegnących. Zaczniemy od przeprowadzenia ubezwymiarowienia w celu zredukowania liczby parametrów. Przeprowadzając zamianę zmiennych

$\hat{t}=rt,\quad\hat{x}=x\left(\frac{r}{D}\right)^{{1/2}}$

i wracając do standardowych oznaczeń zmiennych dostajemy

$\frac{\partial N(t,x)}{\partial t}=\triangle N(t,x)+N(t,x)\left(1-{N(t,x)}\right).$

(11.5)

Odpowiadające równaniu (11.5) równanie kinetyczne, czyli równanie logistyczne, ma dwa stany stacjonarne, niestabilny $N\equiv 0$ i stabilny $N\equiv 1$ . Będziemy więc poszukiwać fali biegnącej od $N=0$ do $N=1$ .

Definicja 11.1

Falą biegnącą nazywamy nietrywialne, ograniczone rozwiązanie równania (11.3) postaci

$N(t,x)=U(x-ct)=U(z),\quad z=x-ct,$

dla pewnego $c>0$ , które nazywamy prędkością fali.

Zauważmy, że

Uwaga 11.1

Równanie (11.1) nie ma rozwiązań w postaci fali biegnącej.

Załóżmy, że takie rozwiązanie $U(z)$ istnieje. Skoro tak, to spełnia równanie różniczkowe

$D\frac{d^{2}U}{dz^{2}}+c\frac{dU}{dz}=0,$

zatem wartości własne tego równania spełniają równanie kwadratowe $D\lambda^{2}+c\lambda=0$ , czyli rozwiązania mają postać $U(z)=A+B\exp\left(-\tfrac{c}{D}z\right)$ , gdzie $A$ , $B$ są stałymi. Ponieważ fala biegnąca jest ograniczona dla wszystkich $z\in\mathbb{R}$ , więc $B=0$ , czyli $U(z)=A$ , zatem nie jest to rozwiązanie w postaci fali.

∎

Dla równania Fishera – Kołmogorowa (11.5) rozwiązanie w postaci fali biegnącej musi spełniać równanie

$\frac{d^{2}U}{dz^{2}}+c\frac{dU}{dz}+U(1-U)=0$

(11.6)

i postulujemy, że $\lim _{{z\to-\infty}}U(z)=0$ , $\lim _{{z\to+\infty}}U(z)=1$ , zatem fala łączy dwa stany stacjonarne — jest odpychana od niestabilnego i przyciągana przez stabilny.

Równanie (11.6) jest równoważne układowi równań

$\begin{array}[]{lcl}U^{{\prime}}&=&V,\\ V^{{\prime}}&=&-cV-U(1-U),\end{array}$

(11.7)

gdzie prim ( ${}^{{\prime}}$ ) oznacza różniczkowanie po $z$ . W przestrzeni fazowej $(U,V)$ układ (11.7) ma oczywiście dwa punkty krytyczne $(0,0)$ oraz $(1,0)$ . Macierz Jacobiego tego układu ma postać

$\text{MJ}(U,V)=\left(\begin{array}[]{cc}0&1\\ -1+2U&-c\end{array}\right).$

Badając lokalną stabilność wyznaczamy wielomian charakterystyczny

$\lambda^{2}+c\lambda+1-2U=0$

i wartości własne

	$\displaystyle\lambda _{{1,2}}=\frac{1}{2}\left(-c\pm\sqrt{c^{2}-4}\right)\quad\text{dla}\quad(0,0),$
	$\displaystyle\lambda _{{1,2}}=\frac{1}{2}\left(-c\pm\sqrt{c^{2}+4}\right)\quad\text{dla}\quad(1,0),$

skąd wynika, że $(1,0)$ jest zawsze punktem siodłowym, natomiast $(0,0)$ jest stabilnym ogniskiem dla $c^{2}<4$ i węzłem stabilnym dla $c^{2}>4$ .

Stwierdzenie 11.1

Dla $c\geq 2$ równanie (11.5) ma rozwiązanie w postaci fali biegnącej o prędkości $c$ .

Dla układu (11.7) przy $c\geq 2$ z kierunku pola wektorowego łatwo wnioskujemy, że istnieje orbita łącząca punkt krytyczny $(0,0)$ z punktem krytycznym $(1,0)$ , przy czym można pokazać, że pozostaje ona w dodatniej półpłaszczyźnie, por. rys. 11.2. Orbita ta tworzy rozwiązanie w postaci fali biegnącej łączącej punkty $(0,0)$ oraz $(1,0)$ .

∎

$Portret fazowy układu~\eqref{fala-uklad} dla $c=3$ (na górze) oraz dla $c=1$ (na dole).$

Rys. 11.2. Portret fazowy układu (11.7) dla $c=3$ (na górze) oraz dla $c=1$ (na dole).

Uwaga 11.2

Dla $c<2$ także istnieje orbita łącząca dwa punkty krytyczne $(0,0)$ i $(1,0)$ , ale ze względu na typ punktu $(0,0)$ nie jest zachowana nieujemność rozwiązań, zatem takie rozwiązania są nieistotne z biologicznego punktu widzenia, por. rys. 11.2.

Zagadnienie występowania fal biegnących (bądź też fal innego typu) jest często badane w kontekście tego typu równań. Należy jednak zdawać sobie sprawę z tego, że nie zawsze rozważana dziedzina jest na tyle duża, że można założyć iż jest to cała przestrzeń. Wtedy, przy założeniu ograniczonej dziedziny, rozważamy zagadnienia początkowo-brzegowe.

$Przykładowe rozwiązanie równania~\eqref{fisher-skal} w postaci fali biegnącej o prędkości $c=3$.$

Rys. 11.3. Przykładowe rozwiązanie równania (11.5) w postaci fali biegnącej o prędkości $c=3$ .

Równanie (11.4) było także rozpatrywane w kontekście rozwoju nowotworu przy założeniu warunku brzegowego. Załóżmy dla uproszczenia jednowymiarową propagację komórek na odcinku $[0,1]$ z warunkiem początkowym

$N(0,x)=N_{0}(x)$

(11.8)

i warunkiem brzegowym Neumanna, czyli

$\tfrac{\partial N}{\partial x}\left|{}_{{x=0,1}}=0\right.$

(11.9)

Dla zagadnienia początkowo-brzegowego (11.4) (11.8) (11.9) lokalne istnienie i jednoznaczność wynika z własności kinetyki.

Korzystając z twierdzenia porównawczego 11.1 i własności rozwiązań równania logistycznego możemy wykazać, że dla warunku początkowego $N_{0}(x)\in(0,1)$ , $x\in[0,1]$ rozwiązanie zagadnienia początkowo-brzegowego (11.4) (11.8) (11.9) pozostaje w przedziale $[0,1]$ .

Twierdzenie 11.1

Niech $u(t,x)$ , $v(t,x)$ będą ograniczonymi funkcjami ciągłymi spełniającymi nierówności

$\begin{split}&A\leq u(t,x)\,,\ v(t,x)\leq B\,,\\ &\frac{\partial u}{\partial t}-\frac{\partial^{2}u}{\partial x^{2}}-f(u)\geq\frac{\partial v}{\partial t}-\frac{\partial^{2}v}{\partial x^{2}}-f(v)\quad w\quad[0\,,T]\times(a,b)\,,\\ &A\leq v(0,x)\leq u(0,x)\leq B\quad w\quad(a,b)\,,\\ \end{split}$

gdzie $-\infty\leq a<b\leq\infty$ i $0<T\leq\infty$ oraz $f$ jest funkcją klasy $\mathbf{C}^{1}$ . Załóżmy także, że

jeśli $a>-\infty$ , to

$A\leq v(t,a)\leq u(t,a)\leq B\quad w\quad[0,T];$
jeśli $b<\infty$ , to

$A\leq v(t,b)\leq u(t,b)\leq B\quad w\quad[0,T].$

Wtedy albo $u\equiv v$ w $[0\,,T]\times[a,b]$ , albo $u>v$ w $[0\,,T]\times(a,b)$ .

Badając dynamikę modelu zwykle zaczynamy od zbadania stabilności rozwiązań stacjonarnych. Dla modeli niejednorodnych przestrzennie, oprócz rozważanych do tej pory rozwiązań stacjonarnych niezależnych od przestrzeni możemy mieć także rozwiązania stacjonarne, które zależą od położenia w dziedzinie równania. Jednak typowo analizę takiego zagadnienia zaczynamy od zbadania istnienia i stabilności rozwiązań stacjonarnych jednorodnych przestrzennie, a potem przechodzimy do znacznie trudniejszych zagadnień związanych z niejednorodnością przestrzenną. Jednym z takich zagadnień jest powstawanie tzw. wzorów Turinga — od strony analitycznej wiąże się to z bifurkacją stabilnego przestrzennie stanu stacjonarnego do niestabilnego, kiedy w otoczeniu stanu stacjonarnego powstają skupiska osobników — grupy o różnym zagęszczeniu — tworzące wzory przestrzenne.

Niech $\eta=N-\bar{N}$ będzie odchyleniem rozwiązania równania (11.4) od jednorodnego przestrzennie stanu stacjonarnego $\bar{N}$ (gdzie $\bar{N}=0$ lub $\bar{N}=1$ ). Równanie zlinearyzowane wokół $\bar{N}$ ma postać

$\eta _{t}=a\eta+D\eta _{{xx}},$

(11.10)

gdzie $a=r$ dla stanu $\bar{N}=0$ i $a=-r$ dla stanu $\bar{N}=1$ , przy czym równanie (11.10) badamy z warunkiem brzegowym $\tfrac{\partial\eta}{\partial x}\left|{}_{{x=0,1}}=0\right.$ i początkowym $\eta(0,x)=\eta^{0}$ . Z twierdzenia o linearyzacji wynika, że jeśli kinetyka jest funkcją różniczkowalną, to ze stabilności rozwiązania równania zlinearyzowanego wynika stabilność rozwiązania równania wyjściowego. Badając stabilność równania liniowego szukamy rozwiązań postaci $\eta(t,x)=\eta _{0}\exp(\lambda t)W(x)$ , gdzie $W^{{\prime\prime}}(x)+k^{2}W(x)=0$ , $k\in\mathbb{Z}$ , jest rozwiązaniem fundamentalnym równania dyfuzji. Istnienie takiego rozwiązania wymaga spełnienia zależności $\lambda=a-Dk^{2}$ . Zatem dla $a=-r$ mamy $\lambda=-r-Dk^{2}<0$ , czyli $\bar{N}=1$ jest rozwiązaniem stabilnym, gdyż wszystkie wartości własne są ujemne, zaś dla $\bar{N}=0$ mamy $\lambda=r-Dk^{2}$ i przy $k=0$ wartość własna $\lambda _{0}=r>0$ , zatem to rozwiązanie jest niestabilne. Dokładnie tak samo kształtuje się stabilność w przypadku bez dyfuzji, zatem dla równania Fishera – Kołmogorowa formowanie się wzorów Turinga nie jest możliwe. Co więcej widać, że aby takie wzory powstały, model musi być opisany za pomocą układu co najmniej dwóch równań.

W przypadku równania (11.4) możemy udowodnić znacznie więcej, czyli globalną stabilność $\bar{N}=1$ przy założeniu, że $N$ jest oddzielone od zera.

Twierdzenie 11.2

Jeśli $N_{0}(0)\geq d>0$ dla pewnego $d\in(0,1)$ oraz $t\geq 0$ , to $N(t,x)\to 1$ przy $t\to\infty$ .

Twierdzenie porównawcze pokazuje, że $N(t,x)\geq d$ dla $t\geq 0$ i $x\in[0,1]$ . Jeśli $\eta=N-1$ , to

$\dot{\eta}=-r\eta(\eta-1)+D\triangle\eta,\quad\frac{\partial\eta}{\partial x}\left|{}_{{x=0,1}}\right.=0$

(11.11)

oraz $\eta(t,x)\geq-(1-d)$ dla $t\geq 0$ i $x\in[0,1]$ . Zdefiniujmy

$L(t)=\frac{1}{2}\int\limits _{{[0,1]}}\eta^{2}(t,x)dx.$

Wtedy

$\dot{L}(t)=\int\limits _{{[0,1]}}\dot{\eta}(t,x)\eta(t,x)dx.$

Mnożąc równanie (11.11) przez $\eta$ i całkując przez części dostajemy

$\dot{L}(t)=-r\int\limits _{{[0,1]}}\eta^{2}(t,x)dx-r\int\limits _{{[0,1]}}\eta^{3}(t,x)dx-D\int\limits _{{[0,1]}}\left(\frac{\partial\eta(t,x)}{\partial x}\right)^{2}dx.$

Widzimy, że jeśli $N>1$ , to $\eta>0$ oraz $\dot{L}(t)<-2rL$ . Jeśli zaś $N<1$ , to $0>\eta>-(1-d)$ i $\eta^{3}>-(1-d)\eta^{2}$ , skąd wynika $\dot{L}(t)<-2rdL$ . W obu przypadkach $L(t)\to 0$ przy $t\to\infty$ , skąd wynika $N(t,x)\to 1$ .

∎

11.2. Wzory Turinga

Jak już wspomnieliśmy, w przypadku jednego równania kinetycznego formowanie się wzorów Turinga nie jest możliwe, gdyż stabilność nie zmienia się pod wpływem dyfuzji. Natomiast dla $n=2$ istnieje taka możliwość. W szczególności można w taki sposób opisać tworzenie się wzorów umaszczenia u ssaków i próbować na tym gruncie wyjaśnić, dlaczego zebry czy żyrafy mają ustalony wzór umaszczenia, a u kotów domowych ten wzór jest zmienny. Podobnie — dlaczego wzór na tułowiu różni się od wzoru na kończynach lub ogonie. Taki model nazywamy modelem melanogenezy, ponieważ zakłada się, że za umaszczenie odpowiada stężenie barwnika zwanego melaniną. Melanina powstaje w reakcji substratu (tyrozyny) z enzymem (tyrozynazą). Model może być zatem opisany za pomocą układu dwóch równań reakcji-dyfuzji. Ponieważ w przypadku tułowia w uproszczeniu możemy powierzchnię skóry rozpatrywać jako pewną figurę płaską, to taki układ rozważa się najczęściej w ograniczonym obszarze w $\mathbb{R}^{2}$ utożsamianym zwykle z prostokątem. Zadaje się zatem warunek początkowy oraz warunek brzegowy, który tu jest warunkiem Neumanna. Jeśli natomiast rozpatrujemy kończyny lub ogon, to rozpatrujemy taki układ np. na walcu lub stożku.

Badając możliwość występowania wzorów Turinga przeanalizujemy lokalną stabilność jednorodnych rozwiązań stacjonarnych w układzie dwóch równań reakcji-dyfuzji

$\begin{array}[]{lcl}s^{1}_{t}&=&f_{1}(s^{1},s^{2})+D_{1}\triangle s^{1},\\ s^{2}_{t}&=&f_{2}(s^{1},s^{2})+D_{2}\triangle s^{2},\end{array}$

(11.12)

gdzie $s^{i}$ , $i=1,2$ , są stężeniami substancji chemicznych, $f_{i}$ oznaczają tempa produkcji tych reakcji, a $D_{i}$ są ich współczynnikami dyfuzji. W przypadku melanogenezy funkcje produkcji mają postać

$f_{1}(s^{1},s^{2})=s^{1}_{0}-s^{1}-\frac{\varrho s^{1}s^{2}}{1+s^{1}+k(s^{1})^{2}},\quad f_{2}(s^{1},s^{2})=\alpha(s^{2}_{0}-s^{2})-\frac{\varrho s^{1}s^{2}}{1+s^{1}+k(s^{1})^{2}},$

gdzie $s^{1}$ to stężenie tyrozyny, a $s^{2}$ — tyrozynazy. Zauważmy, że dla $k=0$ nieliniowy składnik $f_{i}$ ma postać funkcji Michaelisa – Menten, natomiast dla $k>0$ wraz z rosnącą ilością substratu reakcja ulega zahamowaniu.

Dla uproszczenia układ (11.12) będziemy rozpatrywać z jednorodnym warunkiem brzegowym Neumanna na $[0,\pi]$ — analiza w przypadku kwadratu $[0,\pi]\times[0,\pi]$ przebiega podobnie. Także dla uproszczenia założymy, że $(0,0)$ jest jednorodnym stanem stacjonarnym układu (11.12) (jeśli stanem stacjonarnym jest $(\bar{s}^{1},\bar{s}^{2})$ , to dokonujemy zamiany zmiennych $s^{i}\rightarrow s^{i}-\bar{s}^{i}$ ).

Badając stabilność stanu stacjonarnego zlinearyzujemy układ (11.12)

$\begin{array}[]{lcl}s^{1}_{t}&=&\alpha _{{11}}s^{1}+\alpha _{{12}}s^{2}+D_{1}\triangle s^{1},\\ s^{2}_{t}&=&\alpha _{{21}}s^{1}+\alpha _{{22}}s^{2}+D_{2}\triangle s^{2},\end{array}$

(11.13)

przy czym $\alpha _{{ij}}=\tfrac{\partial f_{i}}{\partial s^{j}}$ , $i,j\in\{ 1,2\}$ . Wzory Turinga mogą pojawić się w sytuacji, gdy w modelu bez dyfuzji, czyli dla $D_{1}=D_{2}=0$ stan stacjonarny jest stabilny, a dla pewnych $D_{i}>0$ destabilizuje się.

Okazuje się, że liniowe równanie (11.13) można rozwiązać metodą rozdzielenia zmiennych, czyli $s^{i}(t,x)=F(t)G(x)$ , gdzie $F(t)$ to rozwiązanie układu kinetycznego (dla $D_{i}=0$ ), a $G(x)$ jest rozwiązaniem zagadnienia własnego operatora Laplace'a. W ogólnym przypadku na odcinku $[0,\pi]$ takie rozwiązanie jest kombinacją sinusów i kosinusów, jednak ponieważ badamy układ (11.13) z warunkiem brzegowym Neumanna, to spełnia go tylko kosinus. Stąd

$s^{i}(t,x)=A_{i}\text{e}^{{Bt}}\cos(Cx),\quad i=1,\ 2.$

(11.14)

Podstawiając postać rozwiązania (11.14) do układu (11.13) otrzymujemy

$\begin{array}[]{lcl}A_{1}B&=&\alpha _{{11}}A_{1}+\alpha _{{12}}A_{2}-D_{1}C^{2}A_{1},\\ A_{2}B&=&\alpha _{{21}}A_{1}+\alpha _{{22}}A_{2}-D_{2}C^{2}A_{2},\\ \end{array}$

(11.15)

skąd dostajemy układ równań na współczynniki $A_{1}$ , $A_{2}$

$\begin{array}[]{lcl}(B-\alpha _{{11}}+D_{1}C^{2})A_{1}-\alpha _{{12}}A_{2}&=&0,\\ -\alpha _{{21}}A_{1}+(B-\alpha _{{22}}+D_{2}C^{2})A_{2}&=&0.\\ \end{array}$

(11.16)

Rozwiązanie dostajemy przyrównując wyznacznik układu (11.16) do 0

$B^{2}+B(-\alpha _{{11}}-\alpha _{{22}}+C^{2}(D_{1}+D_{2}))+(\alpha _{{11}}-D_{1}C^{2})(\alpha _{{22}}-D_{2}C^{2})-\alpha _{{12}}\alpha{21}=0.$

(11.17)

W przypadku bez dyfuzji równanie (11.17) ma postać

$B^{2}-B(\alpha _{{11}}+\alpha _{{22}})+\alpha _{{11}}\alpha _{{22}}-\alpha _{{12}}\alpha{21}=0$

i stabilność stanu stacjonarnego zależy od śladu i wyznacznika macierzy Jacobiego układu kinetycznego. Precyzyjniej, jeśli

$\alpha _{{11}}+\alpha _{{22}}<0$ ,
$\alpha _{{11}}\alpha _{{22}}-\alpha _{{12}}\alpha{21}>0$ ,

to stan stacjonarny dla $D_{1}=D_{2}=0$ jest asymptotycznie stabilny. Jeśli teraz chcemy, by był on niestabilny dla $D_{i}>0$ , to któraś z powyższych nierówności powinna mieć przeciwny znak. Nierówność I w przypadku z dyfuzją przyjmuje postać $\alpha _{{11}}+\alpha _{{22}}-C^{2}(D_{1}+D_{2})<0$ , więc dla dowolnych $D_{i}>0$ jest spełniona przy założeniu $\alpha _{{11}}+\alpha _{{22}}<0$ . Zatem tylko nierówność II może zmienić znak dla $D_{i}>0$ . Postulując

$(\alpha _{{11}}-D_{1}C^{2})(\alpha _{{22}}-D_{2}C^{2})-\alpha _{{12}}\alpha{21}<0$

dostajemy nierówność na $C^{2}$ , skąd otrzymujemy warunki destabilizacji stanu stacjonarnego pod wpływem dyfuzji

$\alpha _{{11}}+\alpha _{{22}}<0$ ,
$\alpha _{{11}}\alpha _{{22}}-\alpha _{{12}}\alpha{21}>0$ ,
$D_{2}\alpha _{{11}}+D_{1}\alpha _{{22}}>2\sqrt{D_{1}D_{2}(\alpha _{{11}}\alpha _{{22}}-\alpha _{{12}}\alpha{21})}$ .

W analogiczny sposób możemy badać tworzenie się wzorów przestrzennych w modelu melanogenezy, przy czym wielkość i różnorodność formujących się wzorów zależy oczywiście od parametrów modelu. Jeśli dla danego gatunku dominują składniki kinetyczne, to powstające wzory są powtarzalne (jak np. u żyraf czy zebr), a jeśli składniki losowe, to wzory mogą się bardzo różnić w zależności od osobników (jak w przypadku kotów domowych).

$Graficzne przedstawienie rozwiązań układu~\eqref{tur} dla różnych układów parametrów. Ciemne plamy odzwierciedlają stężenie substratu $s^{1}$ przewyższające wartość rozwiązania stacjonarnego, natomiast obszary z jasnymi plamami - stężenie poniżej tej wartości.$

Rys. 11.4. Graficzne przedstawienie rozwiązań układu (11.12) dla różnych układów parametrów. Ciemne plamy odzwierciedlają stężenie substratu $s^{1}$ przewyższające wartość rozwiązania stacjonarnego, natomiast obszary z jasnymi plamami - stężenie poniżej tej wartości.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Modelowanie matematyczne w biologii i medycynie