Lyhyt johdatus matemaattiseen mallinnukseen osa 2: bifurkaatiot ja kaaosteoria

Isokliinit, virtaus- ja vaihekartat

Lotka-Volterran tasapainopisteitä ovat alfa/beta ja gamma/sigma. Piirretään näiden pohjalta isokliinit (alla), ja huomataan, että kun tasapainopiste 1:n yläpuolella populaatio A pienenee ja alapuolella kasvaa, populaatio B taas vähenee tasapainopiste 2:n vasemmalla ja kasvaa oikealla puolella. Kahden populaation dynamiikan yhteisvaikutusta on kuvattu mustalla nuolella:

Luciano da F. Costa on piirtänyt Lotka-Volterran dynamiikasta virtaus- ja vaihekartan, josta huomataan virtaussuunnan seuraavan isokliinien avulla saatua vaikutelmaa. Costan mallissa kaikki parametrit ovat 1, jolloin myös kahden isokliinin leikkauspiste on (1,1) ja kyseessä on epästabiili tasapainopiste. Koska pisteessä ominaisarvot olivat imaginäärisiä ilman reaalisia osia, on järjestelmä vaimentamaton oskillaattori.

Bifurkaatiot

Bifurkaatioilla tarkoitetaan väljimmillään yllättävää muutosta funktion laadullisessa käytöksessä. Satula-solmu-bifurkaatiossa kaksi tasapainopistettä kohtaavat ja eliminoivat toisensa, transkriittisessä bifurkaatiossa tasapainopisteiden stabiilisuusominaisuudet vaihtuvat, ja kahdentumisbifurkaatiossa tasapainopisteiden määrä vaihtelee arvojen muutoksen mukaan. Näiden lisäksi on myös Hopfin bifurkaatiossa, jossa tasapainopisteiden stabiilisuus- sekä jaksollisuusominaisuudet vaihtuvat. Aiheesta löytyy lisää informaatiota sekä muodolliset todistukset Parvisen luentomonisteesta (kevät 2022).

Seuraavassa nähtävän logistisen kartan ja bifurkaation tarkemmat laskennalliset perusteet löytyvät täältä. Aiemmin käsiteltiin logistista kasvumallia
\begin{equation}\frac{dN}{dt}=F(N)=rN(1-\frac{N}{K}). \end{equation}
Kirjoitetaan tässä logistinen funktio rekursiomuodossa
\begin{equation}x_{t+1}=rx_t(1-x_t),\end{equation}
jolla on haarukkabifurkaatiot kun r= 3,3.449, 3.544... Hyvin yksinkertaisena ideana on laskea funktion rekursiokaavan avulla uusia ilmentymiä ajassa t, laskea kiintopisteet ja arvioida niiden stabiilius sekä käyttäytyminen eri kasvuvakio r:n arvoilla. Kuten aiemmassa linkissä huomataan, alkaa rekursiofunktio monimutkaistuessaan saada yhä useampia juuria, joiden vuoksi koko yhtälö alkaa hajota toistuvien bifurkaatioiden sekä kiintopistedynamiikan vuoksi yhä kaoottisemmaksi ja hallitsemattomammaksi sotkuksi.Seuraavassa nähtävän logistisen funktion dynamiikan bifurkaatioiden Feigenbaumin kuvaajan R-koodi löytyy täältä

Kiehtovinta logistisessa funktiossa on, että vaikka sillä onkin kaoottisia vaiheita ja näennäisen loputtomasti bifurkaatioita, löytyy r=3 ja r<4 väliltä myös vakaita, ei-kaoottisia ja ei-bifurkaatioita aiheuttavia poikkeamia. Kun r >= 4, populaatio romahtaa. Logistisen bifrukaatiokasauman kiehtovista poikkeamista, kaaoksen teoriasta ja säännönmukaisuuksista on myös erinomainen Wikipedia-artikkeli myöhempää käyttöä varten.

Kaaos

Arkikielessä kaaoksella tarkoitetaan löyhästi jonkinlaista hallitsematonta epäjärjestystä. Matematiikassa kaaosta koskeva teoria on huomattavasti "tylsempää" ja määritelmät suurpiirteisempiä. Määritelmistä puhutaankin tässä kohtaa monikossa, sillä niitä riittää aiheelleen uskollisesti pökerryttävän suuri määrä. Seuraavassa muutama esimerkki:

(Kaaos Lin ja Yorken mukaan): Olkoon J reaaliakselin väli. Kuvaus f : J → J
on kaoottinen Lin ja Yorken mielessä, jos se on jatkuva ja sillä on kolmejaksoinen rata. Toisinaan vaaditaan lisäksi herkkä riippuvuus alkuehdoista.

(Kaaos Devaneyn mukaan). Olkoon X metrinen avaruus. Jatkuva kuvaus
f : X → X on kaoottinen Devaneyn mielessä joukossa X, jos seuraavat kolme ehtoa ovat voimassa:

1. f on topologisesti transitiivinen
2. kuvauksen f jaksolliset pisteet ovat tiheässä joukossa X
3. kuvauksella f on herkkä riippuvuus alkutilasta

Yllä mainitut määritelmät löytyvät tarkempine lisätietoineen Parvisen luentomonisteesta (kevät 2022).

Näiden lisäksi Feldman antaa kirjassaan Chaos and Dynamical Systems vielä seuraavan löyhän määritelmän: dynaaminen systeemi on kaoottinen, jos

1. sen aikariippuvainen kulku on deterministinen funktio
2. sen radat ovat rajoitettuja
3. sen radat ovat herkkiä lähtöoletuksille
4. sen radat ovat epäsäännöllisiä

Erikseen on mainittava, että vaikka kolme eri määritelmää tietyssä mielessä risteävätkin, ei yhden määritelmän mukainen kaaos välttämättä ole sitä kahden muun mukaan.

Ehkä kuuluisin kaikista kaoottisista epälineaarisita dynaamisista järjestelmistä on Lorenzin systeemi
\begin{equation}\frac{dx}{dt}\sigma (y-x)\end{equation}
\begin{equation}\frac{dy}{dt}x(\rho-z)-y\end{equation}
\begin{equation}\frac{dz}{dt}xy-\beta z,\end{equation}
joka kuvaa tasaisen, alta lämpenevän ja yltä viilenevän nestekerroksen dynamiikkaa. Kyseessä on epälineaarinen, epäjaksollinen, kolmiulotteinen deterministinen systeemi, jolla on kokonaiset kolme tasapainopistettä. Kun p=1, on funktiossa haarukkabifurkaatio. Kun taas p=28, sigma=10 ja beta = 8/3, on systeemillä kaoottisia ratkaisuja. Tällöin systeemi jää kiertämään äärettömästi Lorentzin attraktorin, kirjallisuudesta riippuen myös "oudon attraktorin" ympärille, ja systeemi hakeutuu kaikkiin kolmeen tasapainopisteeseen. Tällöin myös yllä annetut määritelmät pätevät, ja ratojen kuvaajaksi saadaan Wikipedian koodilla

Kaaosteorian attraktorit ovat oma kiehtova ryhmänsä, josta voisi kirjoittaa vaikka kuinka paljon. Tässä kuitenkin vain koskettiin aiheeseen näin yleisellä tasolla, joten ei siitä sen enempää. Ehkä joskus myöhemmin kirjoittanen aiheesta lisää.

Kirjasuosituksia

Viimeiset pari päivää on pitkälti mennyt seuraavaa aihetta opiskellessa ja kirjoja lukiessa. En vielä vähään aikaan varmaankaan saa kirjoitettua kokonaista matemaattista postausta, joten yhtä hyvin voinen kertoa lähemmin lukemastani.

Sarah P. Otto, Troy Day: A Biologist's Guide to Mathematical Modeling in Ecology and Evolution

Tämä likemmäs tuhatsivuinen järkäle on oikea sovelletun matematiikan aarreaitta. Ensiksikin: tästä kirjasta löytyy perusteet kaikkeen, mitä tähänkin asti olen blogissani käsitellyt. Alkaen ekologiasta sekä muutamista perusmalleista aina Markov-ketjuihin ja stokastisiin malleihin asti, kaikkea löytyy, mutta ilman kaikenlaista ylimääräistä teoriaa. Suosittelenkin tätä kirjaa lämpimästi kevyeksi sisäänheittotuotteeksi kaikille sovelletusta matematiikasta ja ekologiasta. Henkilökohtaisesti olisin kaivannut enemmän tietoa peliteoriasta sekä sen biologisista sovelluksista, mutta tätä varten löysin toisen kirjan....

John Maynard Smith: Evolution and the Theory of Games

Tässä on kokonainen kirja täynnä ekologista peliteoriaa. Smithin kirja on 70-luvulta, joten siitä puuttuu monia nykyajan hienouksia, mutta aiheensa johtavana ytimekkäänä kokonaisesityksenä Evolution and the Theory of Games puolustaa paikkaansa yhä edelleen.


David P. Feldman: Chaos and Dynamical Systems

Chaos and Dynamical Systems alkaa siitä, mihin aiempi biologinen kirja jäi: determinististen funktioiden kaoottiseen hajoamiseen. Näitä juttuja käsiteltiin lyhyesti yliopistokursseilla, muttei ikinä menty pintaa syvemmälle. Tämä kirjanen on erittäin selkeä ja helppolukuinen opas aiheeseen ilman kaikenlaista syvempää matemaattista salatieteilyä. Verrattuna erääseen toiseen lukemaani saman aiheen oppikirjaan tämä on ehkä hivenen verran yliyksinkertaistettu.

SIS ja SIRVI-malli

SIS-malli

Lähde löytyy täältä.

SIS-malli kuvaa tilannetta, jossa taudin sairastaminen ei tuotakaan pysyvää immuniteettia, vaan henkilö palaakin takaisin lähtötilanteeseen eli taudille alttiiksi.
\begin{equation}I(t+\triangle t) = I(t)+\beta \triangle t S(t) I(t) - \gamma \triangle t I(t), \end{equation}
kun kolmio-t lähestyy nollaa, saadaan differentiaaliyhtälö
\begin{equation}\frac{dI}{dt} =\frac{\beta SI}{N} - \frac{\gamma I}{N}. \end{equation}
sekä vastineeksi
\begin{equation} \frac{dS}{dt}=-\frac{\beta SI}{N}+\frac{\gamma I}{N}.\end{equation}

Oletus I+S=N saadaan käännettyä muotoon S=N-I, sekä sijoittamalla saatu muoto edelliseen differentiaalityhtälöön saadaan
\begin{equation}\frac{dI}{dt} \beta (N-I)I - \gamma I. \end{equation}
sekä edelleen muotoilemalla
\begin{equation}\frac{dI}{dt} = (\beta N - \gamma)I(1-\frac{\beta I}{N-\gamma}), \end{equation}
joka on siis logistinen kasvumalli.

Edellinen differentiaaliyhtälö voidaan edelleen ratkaista. Ratkaisun pituuden ja monimutkaisuuden vuoksi en käy ratkaisumenetelmää yksityiskohtaisesti tässä läpi, täydellinen ratkaisu löytyy tästä Khanin, Sabbirin ja Sadiqin artikkelista. Differentiaaliyhtälön ratkaisu lähtee ylempänä muotoillun yhtälön integroinnista. Kun muuttujan I differentiaaliyhtälö on saatu ratkaistua, saadaan muuttujan S ratkaisu yksinkertaisesti käyttämällä identiteettiä S = N-I, jolloin ratkaisuiksi saadaan
\begin{equation}S=N - \frac{\beta N - \gamma}{\beta - (\beta N - \gamma)\frac{\beta N - \gamma - i_0 \beta}{(\beta N - \gamma)i_0}e^{-(\beta N - \gamma )t}}. \end{equation}
sekä muuttujan I yhtälöksi
\begin{equation}I=\frac{\beta N - \gamma}{\beta - (\beta N - \gamma)\frac{\beta N - \gamma - i_0 \beta}{(\beta N - \gamma)i_0}e^{-(\beta N - \gamma )t}}. \end{equation}

Skipataan tässä yksinkertaisemmassa mallissa tasapainopisteiden haku sekä niiden stabiiliusominaisuuksien analyysi.  Aiheesta löytyy erinomainen esitelmä täältä.
 
Tartuntoja kuvaavien differentiaaliyhtälöiden kuvaaja on

SIS-mallissa perustartuntaluku on muotoa
\begin{equation}R_0 = \frac{\beta N}{\gamma}, \end{equation}
eli jos R_0 > 1, epidemia käynnistyy. Tällöin sairastuneella on vähintään yksi kontakti, ja tauti voi siten lähteä leviämään.

SIRVI-malli 

SIRVI-malli lähtee pitkälti samoista oletuksista kuin kaksi edellistäkin, sillä erotuksella, ettei rokotus tällä kertaa annakaan pysyvää immuniteettiä. SIRVIn asetelma muistuttaa siis vahvasti koronarokotteen vaikutuksia, ja on siksi hyvin ajankohtainen. Koronarokote ei siis antanut kestävää immuniteettia tautia vastaan, estänyt ketään sairastamasta tai levittämästä tautia, mutta vahvisti hetkellisesti immuniteettia tautia vastaan. SIRVIin löytyy erinomaisen perusteellinen lähde täältä, ja käytän sitä myös tämän blogipostauksen pohjana. Kaikki olennainen tieto on peräisin Turkyilmazoglun artikkelista, enkä ota tästä kunniaa itselleni.

SIRVI-malli muodostuu yhtälöistä:

\begin{equation} \frac{dS}{dt}=-\beta SI -vS\end{equation}

\begin{equation} \frac{dI}{dt}=\beta SI + \alpha VI -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I \end{equation}

\begin{equation} \frac{dV}{dt}=v S - \alpha VI \end{equation}

missä alfa, beta, gamma ja v ovat rokotettujen tartuntamäärä, epäiltyjen tartuntamäärä, parantumismäärä ja epäiltyjen rokotusmäärä.

Hetkellä t=0 annettujen funktioiden muuttujille pätee:

\begin{equation} S_0 = 1-\tau \end{equation}

\begin{equation} I_0 = \tau \end{equation}

\begin{equation} R_0 = 0 \end{equation}

\begin{equation} V_0 = 0 \end{equation}.

Tällöin on voimassa myös yhtälö

\begin{equation}S+I+R+V = 1 \end{equation}.
Mallin perustartuntaluku on
\begin{equation}R_0=\frac{S_0 \beta -V_0 \alpha}{\gamma}\end{equation}
sekä tasapainopisteet
\begin{equation}(S,I)=(0,0)\end{equation}
\begin{equation}(V,R)=(0 \leq V, R \leq 1)\end{equation}

SIRVI-mallin Jacobin matiiris on muotoa
\begin{equation}J=\begin{bmatrix}-v & -SB & 0 & 0 \\ 0 & V \alpha + S \beta - \gamma & 0 & 0 \\ v & -V \alpha & 0 & 0 \\ 0 & \gamma & 0 & 0  \end{bmatrix}\end{equation}

SIRVI-mallin stabiilisuusominaisuuksista voidaan ainakin sanoa, että nollasta eroavat ominaisarvot
\begin{equation}\lambda_1 = -v \end{equation}
sekä
\begin{equation}\lambda_2 = V \alpha + S \beta - \gamma, \end{equation}
mistä huomataan, että stabiilisuus pätee rokotetulla ja taudista vapaalla populaatiolla kun
\begin{equation}V \leq \frac{\gamma}{\alpha}\end{equation}

Kuvaaja on siis muotoa

Vaikka tässä käsiteltiinkin vain kolmea erityistä epidemiologista mallia, riittää vastaavia sovellutuksia maailmassa vaikka kuinka paljon. Epidemiologia on yksi sovelletun matematiikan merkittävimmistä haaroista.

Käytetyt koodit, perusten aiempaan SIR-mallissa käytettyyn Python-koodiin:

# SIS-malli
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0 = 1
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 0.09, 2.0
# A grid of time points (in days)
t = np.linspace(0, 160, 160)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma):
    S, I = y
    dSdt = -beta * S * I / N + gamma * I/N
    dIdt = beta * S * I / N - gamma * I/N
    return dSdt, dIdt

# Initial conditions vector
y0 = S0, I0,
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I = ret.T

# Plot the data on three separate curves for S(t), I(t) and R(t)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()

#################################################################

# SIRVI-malli
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0, V0 = 1, 0, 0
# Everyone else, S0, is susceptible to infection initially.
S0= N - I0 - R0 -V0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma, v ,alpha= 0.2, 2.10, 0.03, 0.03
# A grid of time points (in days)
t = np.linspace(0, 160, 160)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma, v, alpha):
    S, I, R, V = y
    dSdt = -beta * S * I/N - v * S
    dIdt = beta * S * I/N + alpha*V*I/N - gamma * I/N
    dRdt = gamma * I/N
    dVdt = v * S-alpha*V*I/N
    return dSdt, dIdt, dRdt, dVdt

# Initial conditions vector
y0 = S0, I0, R0, V0
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma,v,alpha))
S, I, R, V = ret.T

# Plot the data on three separate curves for S(t), I(t) and R(t)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
ax.plot(t,V/1000, 'y', alpha=0.5, lw=2, label = 'vaccinated' )
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()

Vaihtoehtoja Lotka-Volterralle

Diclaimer: Bloggerin matematiikkatyökalu on järjettömän buginen ja epäluotettava. Siksi olenkin tekstissäni pyrkinyt mahdollisuuksien mukaan välttämään ylimääräisten symbolien käyttöä.

Edellisessä, Lotka-Volterran peto-saalis -mallia käsittelevässä postauksessa käytiin läpi mallin historiaa sekä yleisiä ominaispiirteitä ja rajoituksia. Huomattiin, että kyseinen malli perustuu hyvin epärealistisille oletuksille peto- ja saaliseläinten kantojen ominaisuuksista. Käytännössä onkin huomattu, ettei Lotka-Volterra aina pädekään erityisen tarkasti kaikissa tilanteissa, minä vuoksi mallille onkin vuosien varrella kehitetty vaihtoehtoisia, oletuksiltaan joustavampia peto-saalis -malleja. Tässä postauksessa käydään läpi kaksi keskeisintä ja ehkäpä yleisimmin käytettyä: Arditi Ginzburg sekä Rosenzweig-MacArthur.

Mistä Lotka-Volterraa siis tarkalleen ottaen kritisoitiin? Ensinnäkin: Lotka-Volterra lähti oletuksesta, että petoeläimillä on rajaton ruokahalu, mikä ei reaalimaailmassa voi pitää eikä pidäkään paikkanasa. Samoin Lotka-Volterrassa petoeläinten puuttumisen seurauksena saaliseläimen kanta kasvaa Malthusin lain mukaisesti rajattomasti - tämäkään oletus ei välttämättä kestä päivänvaloa. Tässä  erinomaisesa luennossa on käsitelty ekologisten mallien ominaisuuksia syvällisemmin ja kattavammin kuin tässä blogipostauksessa.

Arditi-Ginzburg

Arditi-Ginzburgin peto-saalis -malliksi kutsutaan yhtälöparia

\begin{equation} \frac{dx}{dt}=f(x)x-g(\frac{x}{y})y \end{equation}

\begin{equation} \frac{dy}{dt}=\alpha g(\frac{x}{y})y -\beta y \end{equation}
jossa funktio f() kuvaa saaliseläimen kannan luonnollisia muutoksia, funktio g() liittyy saaliseläinten sekä petojen kantojen suhteeseen, alfa on vakio petojen saaliiksi saamasta lisääntymisedusta sekä beta petojen kuolemaa kuvaava suhdeluku. Arditi-Ginzburg ottaa siis huomioon monipuolisemmin populaatioiden muutoksen muutenkin kuin näiden keskinäisen riippuvuuden kautta. Malli onkin vahvasti suhdelukupainotteinen. Arditin ja Ginzburgin päätyminen mainittuihin suhdelukuihin on pitkä tarina, josta voi lukea tarkemmin täällä, referoin keskeisiä kohtia tänne muutamin täydennyksin ja huomioin. 

Tässä kohtaa on syytä huomata, ettei funktioita f() ja g() oleteta yleisesti annetuiksi, vaan niiden paikalle voidaan sovittaa tarpeen mukaan erilaisia suhdeluvuilla toimivia troofisia funktioita eriävine ominaisuuksineen. Tässä blogipostauksessa on tarkoituksella käytetty mahdollisimman yksinkertaisia troofisia funktioita

\begin{equation}f(x)=\gamma x \end{equation}

sekä

\begin{equation}g(\frac{x}{y})= \frac {\gamma x}{y} \end{equation}

jolloin yhtälöpari voidaan kirjoittaa muotoon

\begin{equation} \frac{dx}{dt}=\gamma x² - \gamma x \end{equation} 

\begin{equation} \frac{dy}{dt}=\alpha \gamma x -\beta y \end{equation}

ja asettamalla derivaatat nolliksi, saadaan tasapainopisteiksi (0,0) , (1,gamma*alfa/beta), joiden muodostama Jacobin matriisi on muotoa

\begin{bmatrix} 2 \gamma x - \gamma & 0 \\ \alpha \gamma & -\beta \end{bmatrix}

mistä huomataan, että mikäli beta ja gamma oletetaan aina nollaa suuremmaksi, ovat matriisin ominaisarvot origossa aina samanmerkkiset, eli järjestelmä on annetuilla arvoilla ja oletuksilla lähtökohtaisesti epävakaa. Toisen tasapainopisteen kohdalla itseisarvoista on toinen positiivinen ja toinen negatiivinen, jolloin kyseessä on epävakaa satulapiste. Ominaisarvoista ja stabiiliudesta lisää tämän linkin takaa.

Demonstroisin mieluusti mallia soveltuvalla kuvaajalla, mutta annetuilla oletuksilla ja valinnoilla kannat nousevat aluksi vain pudotakseen nopeasti nollaan, joten vaikka esimerkkinä käytetyt funktiot malliin sopivatkin, eivät ne yhdessä näytä tuottavan juurikaan kelvollista jälkeä. Tässä yksi pätevä opetus siitä, miksi lähtöoletusten kanssa on oltava tarkkana!

Rosenzweig-MacArthur

Rosenzweig-MacArthur -malliksi kutsutaan yhtälöparia

\begin{equation} \frac{dx}{dt}=rx(1-\frac{x}{K})-yh(x) \end{equation}

\begin{equation} \frac{dy}{dt}=y(-\delta+\gamma h(x)) \end{equation}

 Rosenzweig-MacArthur -mallista voi lukea lisää täällä. Yhtenä keskeisenä tavoitteena Rosenzweigilla ja MacArthurilla oli saaliseläimen populaation malthuslaisen kasvun rajoittaminen logistisella kasvufunktiolla. Yksinkertaisuudessaan ylempi differentiaaliyhtälö kuvaa saaliseläimen lisääntymistä ja kuolemaa muista kuin petoeläimestä johtuvista syistä vähennettynä petoeläimen saalistamisesta johtuva väheneminen. Alempi differentiaaliyhtälö kuvaa petoeläinkannan kasvua ja vähenemistä. Tässä mallissa funktio h(x) kuvaa petoeläimen saalistamien saaliseläinten lukumäärää, K populaatioiden koko yhteensä, delta petoläinten poistumaa ja gamma lisääntymistä.

Yllä annetussa artikkelissa on erinomaisen syvällistä analyysia aiheesta, ja tässäkin blogissa tullaan ehdottomasti myöhemmin käsittelemään mainittua teoriaa, mutta tässä vaiheessa tyydymme käymään lyhyesti mallin pääkohdat läpi. Smith yksinkertaistaa yllä annetun kaavan Holling II -tyyppiseksi peto-saalis -malliksi, jolloin sen muodoksi saadaan
\begin{equation} \frac{dx}{dt} = x(1-\frac{x}{K})-\frac{mxy}{1+x} \end{equation}
\begin{equation} \frac{dy}{dt} = y(-c + \frac{mx}{1+x})\end{equation}

Tutkitaan seuraavaksi lyhyesti mallin stabiilisuusominaisuuksia. Tasapainopisteet saadaan yhtälöparista
\begin{equation} 0 = x(1-\frac{x}{K}-\frac{my}{1+x})\end{equation}
\begin{equation} 0 = y(-c + \frac{mx}{1+x})\end{equation}
Tässä tapauksessa tasapainopisteet ovat (0,0), (k,0). Smithin \ määritellään myös eläinten välisen yhteiselon tasapainopiste, mutta emme ota sitä tällä kertaa huomioon. Jacobin matriisiksi saamme tasapainopisteellä (k,0)
\begin{equation}J=\begin{bmatrix} -1 & -\frac{mk}{1+k} \\ 0 & \frac{mk}{1+k}-c \end{bmatrix}\end{equation}
mistä huomataan, että piste on epästabiili, sillä toinen ominaisarvoista on aina positiivinen, ja alempikin on negatiivinen vain jos pätee
\begin{equation}\frac{mk}{1+k}<c.\end{equation}
Tasapainopisteen (0,0) Jacobin matriisi on
\begin{equation}J=\begin{bmatrix} 1 & -m \\ 0 & -c \end{bmatrix},\end{equation}
mistä taas huomataan, että pisteet ovat epästabiileja, sillä toinen ominaisarvo on välttämättä positiivinen, sillä oletuksena c > 0.

Ensimmäisessä kuvaajassa on voimassa
\begin{equation}\frac{mk}{1+k}<c,\end{equation} jolloin piste (k,0) ja (0,0) ovat epästabiileja

Alemmassa epäyhtälö on käännetty, jolloin pisteet (k,0) ja (0,0) edelleen epästabiileja, mutta käyrien muoto on täysin erilainen.

Käytetyt koodit, muokattu Lotka-Volterran mallista:

import numpy as np
import matplotlib
matplotlib.use('tkagg')
import matplotlib.pyplot as plt
from scipy.integrate import odeint

y0 = [10,1] # fish, bears, units in hundreds

t = np.linspace(0,50,num=1000)

m = 1.1
c = 2
k = 100

params = [m,c,k]

def sim(variables, t, params):
    # fish population
    x = variables[0]
    # bear population
    y = variables[1]
    
    m = params[0]
    c = params[1]
    k = params[2]
    
    dxdt = x*(1-x/k) - m*x*y/(1+x)
    dydt = -c*y + m*x*y/(1+x)
    
    return([dxdt,dydt])

y = odeint(sim, y0, t, args=(params,))

f,(ax1,ax2) = plt.subplots(2)

line1, = ax1.plot(t,y[:,0], color="b")

line2, = ax2.plot(t,y[:,1], color="r")

ax1.set_ylabel("Fish (hundreds)")
ax2.set_ylabel("Bears (hundreds)")
ax2.set_xlabel("Time")

plt.show()

Hyvin lyhyt johdatus matemaattiseen mallinnukseen

Mallien johtaminen

Oletetaan, että olemme tehneet joukon havaintoja ilmiöstä A suhteessa B. Huomataan, että havaintojen perusteella voidaan kirjoittaa kaksiulotteinen suora. Havaintojemme perusteella voidaan siis kirjoittaa suoralle yhtälö
\begin{equation}y-y_0 = k(x-x_0) \end{equation}
missä k on suoran kulmakerroin, eli
\begin{equation}k = \frac{\triangle y}{\triangle x}=\frac{y_2-y_1}{x_2-x_1}.\end{equation}
Näin siis yksinkertaisten havaintojen avulla voidaan kirjoittaa ilmiötä yleisesti tai monimutkaisempaa likimain hetkellisesti kuvaava suoran yhtälö, johon havaintoarvoja x asettamalla voidaan arvioida y suuruutta sekä toisin päin.

Otetaan seuraavaksi tarkasteluun logistista mallia käsittelevästä postauksesta tuttu esimerkki eksponentiaalisen kasvun kaavan johtamisesta. Olkoon siis x(t) populaation koko hetkellä t, b syntyneiden määrä ja d kuolleiden määrä hetkellä t. Saadaan
\begin{equation}x(t+\triangle t)=x(t)+b \triangle t x(t) -d \triangle t x(t)  \end{equation}\begin{equation}\frac{x(t+\triangle t)-x(t)}{\triangle t}=(b-d)x(t)  \end{equation}
kun t lähestyy nollaa, saadaan differentiaaliyhtälöksi
\begin{equation}\frac{dx}{dt}=(b-d)x(t)  \end{equation}
sekä differentiaaliyhtälön ratkaisuksi
\begin{equation}x(t)=x_0 e^{rt},  \end{equation}
missä r on kasvuvauhti. Differentiaaliyhtälön ratkaisua kutsutaan eksponentiaaliseksi yhtälöksi. 

Otetaan vielä yksinkertainen diskreetti esimerkki yleisen funktion muotoilemisesta, kun ilmiön tila tunnetaan hetkellä t. Kirjoitetaan
\begin{equation}x(t+1)=x(t)+\alpha, \end{equation}
menemällä ajassa taaksepäin ja sijoittamalla menneitä muuttujia x(t) sisäkkäin, saadaan ensin
\begin{equation}x(t+1)=(x(t-1)+\alpha)+\alpha\end{equation}
\begin{equation}x(t+1)=(x(t-2)+\alpha)+2\alpha,\end{equation}
ja niin edespäin. Menemällä tarpeeksi kauas taaksepäin saadaan yhtälöiden perusteella muotoiltua yleinen muoto:
\begin{equation}x(t+1)=x_0 + (t+1) \alpha \end{equation}

Järjestelmän stabiilisuus

Järjestelmän stabiiliudesta lähde täällä sekä täällä.

Aloitetaan Ljaputinovin stabiiliudesta. Olkoon autonominen dynaaminen systeemi
\begin{equation}x' = f(x(t)), x(0)=x_0, \end{equation}
jolla on olemassa tasapainopiste 
\begin{equation}x_{\epsilon}, f(x_{\epsilon})=0. \end{equation}
Tasapainopistettä sanotaan Ljapunov-stabiiliksi, jos jokaista nollaa suurempaa epsilonia kohden on olemassa myös nolla suurempi delta, jolla pätee kaikilla t >0
\begin{equation} ||x(0)-x_{\epsilon}|| < \delta \Longrightarrow ||x(t)-x_{\epsilon}|| < \epsilon \end{equation}
kaikilla t >0.

Jos tasapainopiste toteuttaa ehdon
\begin{equation} ||x(0)-x_{\epsilon}|| < \delta \Longrightarrow \lim_{t \rightarrow \infty} ||x(t)-x_{\epsilon}|| < \epsilon, \end{equation}
on kyseessä asymptoottisesti stabiili piste.

Asymptoottista tasapainopistettä sanotaan eksponentiaalisesti stabiiliksi, jos mielivaltaisilla muuttujilla alfa,beta ja delta >0 pätee
\begin{equation} ||x(0)-x_{\epsilon}|| < \delta \Longrightarrow ||x(t)-x_{\epsilon}|| < \epsilon \leq \alpha ||x(t)-x_{\epsilon}||e^{- \beta t} \end{equation}
Ehdoista ensimmäinen tarkoittaa, että mikäli yhtälöt lähtevät tasapainopisteen läheisyydestä, ne jäävät kiertämään sen välitöntä ympäristöä kykenemättä erkaantumaan siitä. Asymptoottinen stabiilius tarkoittaa, että yhtälöt eivät vain kierrä tasapainopistettä, vaan myös suppenevat sitä kohti. Eksponentiaalisessa tapauksessa suppenemisnopeus on annettua ylärajaa nopeampaa.

Esimerkkejä Ljapunovin stabiiliuden selvittämisestä löytyy enemmän annetuista lähteistä.

Järjestelmän stabiilisuus voidaan myös selvittää seuraavasti. Olkoon mielivaltainen yhtälöpari
\begin{equation} \frac{dy}{dt}=ax-bz. \end{equation}
\begin{equation} \frac{dz}{dt}=-ax+bz+d. \end{equation}
Ensin selvitetään differentiaaliyhtälöiden tasapainopisteet, eli nollakohdat yllä mainituille yhtälölle:
\begin{equation} 0=ax+b. \end{equation}
\begin{equation} 0 =-ax+bz+d. \end{equation}
Seuraavaksi kirjoitetaan samaisen differentiaaliyhtälökaksikon Jacobin matriisi
\begin{equation}J=\begin{bmatrix} a & 0 \\ a & b \end{bmatrix}\end{equation}
josta laskettujen ominaisarvojen
\begin{equation} det(J-\lambda I) = 0 \end{equation}
avulla saadaan selvitettyä funktion keskeiset stabiilisuusominaisuudet. Eli siis jos jatkuvan dynaamisen järjestelmän ominaisarvot ovat reaalisia ja
\begin{equation}\lambda_{1,2} < 0,\end{equation}
on järjestelmä stabiili, kun taas jos on voimassa
\begin{equation}\lambda_{1,2} > 0\end{equation}
\begin{equation}\lambda_{1} < 0, \lambda_{2}>0,\end{equation}
on järjestelmä epästabiili. Lisäksi jos reaaliosa on nolla, on järjestelmä silloinkin epästabiili. Jos ominaisarvot ovat imaginäärisiä, ja reaalinen osa on negatiivinen, on järjestelmä stabiili. Muutoin kyseessä on epästabiili järjestelmä. Diskreetissä tapauksessa järjestelmä on asymptoottisesti stabiili jos ja vain jos ominaisarvon reaaliselle osalle pätee
\begin{equation}|\lambda_{1,2}|<1.\end{equation}
Jatkuvassa tapauksessa järjestelmä on asymptoottisesti stabiili jos ja vain jos ominaisarvojen reaaliosat ovat negatiivisia.

Nämä tiedot riittävät tässä vaiheessa mallinnusta koskevien blogipostausten ymmärtämiseen. Tilastollista mallintamista, optimointia ja todennäköisyyslaskentaa koskevia osuuksia käsitellään erikseen myöhemmin.

SIR- ja SIRV-mallit sekä rokotuksen vaikutuksen mallinnus

SI-mallin johtaminen

Aiheen johdannoksi esitellään Charnovin monisteiden esittelemä tapa tartuntatautimallin määrittelemiseksi. Aloitetaan seuraavasti tunnetuksi oletetusta kaavasta
\begin{equation}I(t+\triangle t) = I(t)+ \beta \triangle t S(t) I(t)\end{equation}
missä t kuvaa aikaa, I(t) tartuntojen määrää ja S(t) alttiiden määrää hetkellä t sekä beta kontaktien määrää aikayksikössä, eli leviämisintensiteettiä. Viemällä I(t) funktion vasemmalle puolelle, jakamalla molemmat puolet kolmio-t:llä, kun kolmio-t lähestyy nollaa, saadaan tartuntojen määrän muutosta kuvaava differentiaaliyhtälö
\begin{equation} \frac{dI}{dt}=\beta SI, \end{equation}
joka on tuttu myöhemmistä aiheen variaatioista. Myöhempiä malleja ei johdeta yhtä perusteellisesti, vaan käytännössä aloitetaan valmiista ja annetuiksi oletetuista differentiaaliyhtälöistä.

SIR-malli

SIR-malli kuvaa epidemian dynamiikkaa olettamalla olevaksi kolme tilaa: alttiiden S, todetun tautitapauksen I sekä parantuneen R. Tässä mallissa ei oteta huomioon syntymiä eikä kuolemia. Kerran saatu immuniteetti ei häviä, eikä samaan tautiin voi sairastua toistamiseen. Näistä oletuksista on erilaisia vapaampia muotoiluja, kuten SIRI-malli, jossa henkilölle ei muodostu kestävää immuniteettiä, SIRD-malli, jossa henkilön kuolema otetaan huomioon, sekä joukko erilaisia rokotusmalleja, joista tässä käsiteltävässä immuniteetti on pysyvä. Vaihtoehtoisia muotoiluja oletuksineen käsitellään erikseen myöhemmin.

Yksinkertainen SIR-malli muodostuu differentiaaliyhtälöistä
\begin{equation} \frac{dS}{dt}=\frac{-\beta SI}{N} \end{equation}

\begin{equation} \frac{dI}{dt}=\frac{\beta SI}{N} -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I. \end{equation}
SIR-mallin parametreista beta kuvaa tartuntaintensiteettiä ja gamma poistumis- tai parantumisintensiteettiä riippuen miten mallia tulkitaan. Kahden parametrin avulla voidaan edelleen määritellä perustarttuvuusluku
\begin{equation} R_0 = \frac{\beta}{\gamma}, \end{equation}
joka karkeasti tarkoittaa odotettujen tartuntojen määrän. Seuraavaksi huomataan, että sijoittamalla perustarttuvuusluku tartuntojen muutosta kuvaavaan differentiaaliyhtälöön, saadaan
\begin{equation}\frac{dI}{dt}=(R_0 \frac{S}{N}-1)\gamma I, \end{equation}
mistä edelleen huomataan, että kun pätee
\begin{equation}R_0 S > N, \end{equation}
on infektioiden muutos samalla kasvava, eli tartuntojen määrä lisääntyy ja epidemia on kiihtymisvaiheessa. Käänteisessä tapauksessa tartuntojen määrä vähenee ja epidemia on asettumassa. Sama voitaisiin todeta toisin, eli olettamalla likimain jokaisen populaatioon kuuluvan yksilön olevan altis eli S = N, seuraisi siitä, että epidemian käynnistymiseen riittää ehto
\begin{equation}R_0 > 1, \end{equation}
sekä talttumiseen käänteinen tulos.

SIR-mallin kuvaajassa näkyvät eri ryhmien muutosten suunnat sekä kasvunopeudet. Seuraavan SIR-mallin kuvaajaan löytyy koodi täältä.

Kuvaajasta huomataan, että alttiiden määrä vähenee taudin leviämisen sekä siitä parantuneiden kasvamisen mukana. Kaikki muutosnopeudet tasaantuvat ajan mittaan. Samaista koodia sovelletaan mahdollisuuksien mukaan myös saman mallin eri sovelluksissa, kuten seuraavassa SIRV-mallissa.

SIRV-malli

Määritellään johdannoksi rokotusmalliin seuraava tulos. Olkoon S_0 alttiiden prosentuaalinen osuus hetkellä 0, olkoon lisäksi p rokotettujen osuus väestöstä, ja olkoon lisäksi voimassa S_0=(1-p). Jos oletetaan rokotuksen tuoma immuniteetti pysyväksi, voidaan kirjoittaa edellisen kohdan perustarttuvuuslukua hyödyntäen seuraava ehto epidemian käynnistymiselle:
\begin{equation} R_0 (1-p) > 1. \end{equation}
Annetun ehdon avulla voidaan määritellä rokotettavien minimäärä p* epidemian estämiseksi, jolloin voidaan kirjoittaa
\begin{equation} p^* = 1-\frac{1}{R_0}.\end{equation}
SIRV-malli on perusteiltaan kuin aiemmin kuvailtu yksinkertainen SIR-malli, mutta siinä on otettu huomioon rokotuksen vaikutus. Rokotus antaa tässä mallissa pysyvän immuniteetin, eikä tautiin kuolleita oteta tässäkään huomioon.
\begin{equation} \frac{dS}{dt}=\frac{-\beta SI}{N} -vS\end{equation}

\begin{equation} \frac{dI}{dt}=\frac{\beta SI}{N} -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I \end{equation}

\begin{equation} \frac{dV}{dt}=v S \end{equation}
Tässä mallissa parametrit tarkoittavat samaa kuin edellisessä mallissa, mutta uutena SIRV-malliin on lisätty rokotusintensiteettiä kuvaava parametri v. Kun muut parametrit ovat samat kuin ensimmäisessä ja rokotusintensiteetti v = 0.005, saadaan kuvaaja

Kuvasta siis näemme, että rokotuksilla saadaan taudin leviäminen hallintaan ja immuniteetti puolelle samasta määrästä ihmisiä, jotka saivat saman sairastamalla taudin. Infektoiden määrän kasvuvauhtia saatiin näin loivennettua ja huippu madallettua.

SIRV-mallissa käytetty Python-koodi, muokattu aiemmin linkatusta esimerkistä:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

N = 1000
I0, R0, V0 = 1, 0, 0
S0 = N - I0 - R0 -V0
beta, gamma, v = 0.2, 1./10, 0.005
t = np.linspace(0, 160, 160)

def deriv(y, t, N, beta, gamma, v):
    S, I, R, V = y
    dSdt = -beta * S * I / N - v * S
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    dVdt = v * S
    return dSdt, dIdt, dRdt, dVdt

y0 = S0, I0, R0, V0

ret = odeint(deriv, y0, t, args=(N, beta, gamma,v))
S, I, R, V = ret.T

fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
ax.plot(t,V/1000, 'y', alpha=0.5, lw=2, label = 'vaccinated' )
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()

Logistiset kasvumallit

Diclaimer: Bloggerin matematiikkatyökalu on järjettömän buginen ja epäluotettava. Siksi olenkin tekstissäni pyrkinyt mahdollisuuksien mukaan välttämään ylimääräisten symbolien käyttöä.

Lähdetään oletuksesta, että N(t) on populaation koko hetkellä t, b on syntyneiden määrä ja d kuolleiden määrä. Saadaan

\begin{equation}N(t+\triangle t)=N(t)+b \triangle t N(t) -d \triangle t N(t)  \end{equation}

\begin{equation}\frac{N(t+\triangle t)-N(t)}{\triangle t}=(b-d)N(t)  \end{equation}

kun t lähestyy nollaa, saadaan differentiaaliyhtälöksi

\begin{equation}\frac{dN}{dt}=(b-d)N(t)  \end{equation}

sekä differentiaaliyhtälön ratkaisuksi

\begin{equation}N(t)=N_0 e^{rt},  \end{equation}
missä r on kasvuvauhti. Kyseessä on toisin sanoen eksponentiaalisesti kasvava yhtälö. Malthusin lähde, Chasnovin Mathemathical Biology -moniste löytyy ilmaiseksi täältä. Seuraavassa yksinkertainen diskreetti esimerkki eksponentiaalisen kasvun mallin kuvaajasta:

Perusominaisuuksiltaan eksponentiaalinen malli sopii erinomaisesti taloudellisen kasvun ja kassaantumisen kuvaamiseen, mutta resurssien luonnossa asettamien rajoitusten vuoksi erityisen huonosti populaatioiden kasvun mallintamiseen. Tästä syystä Malthusin kasvumallin pohjalta kehitettiin Verholstin malli eli logistinen kasvumalli, jossa otetaan huomioon kasvun yläraja K:

\begin{equation}\frac{dN}{dt}=F(N)=rN(1-\frac{N}{K}). \end{equation}
Verhulstin mallissa N kuvaa populaation kokoa, r kasvuvakiota ja K populaation koon ylärajaa. Logistisen kasvumallin kasvunopeus hidastuu ajan ja populaation koon kasvun mukana. Logistisen kasvumallin differentiaaliyhtälölle voidaan muotoilla seuraavanlainen kuvaaja populaation muutosnopeuden ja koon mukaan:

Logistinen kasvumalli ratkaisemalla differentiaaliyhtälöiden ratkaisumenetelmiä hyödyntäen päädytään kaavaan

\begin{equation}N(t)=\frac{P_0*K*e^{r*t}}{(K-P_0)+P_0*e^{rt}} \end{equation}

jonka pohjalta saadaan muotoiltua populaation koon kuvaaja:

 

Logistiselle kasvumallille saadaan tasapainopisteiksi N=0 ja N=K. Tasapainopisteiden stabiilius saadaan tarkistettua derivoimalla logistinen malli, jolloin saadaan
\begin{equation}F'(N)=r(1-\frac{2N}{K}). \end{equation}

Sovittamalla tasapainopisteet derivoituun yhtälöön, huomataan, että kun N=0, yhtälö saa muodon F'(0)=r ja kun F'(K)=-r. Kyseessä on toisin sanoen satulapiste, jossa nollapiste on epästabiili sekä K stabiili. Populaatio voi siis lähestyä nollapistettä mielivaltaisen lähelle kyeten kuitenkin kasvamaan siitä pois, mutta lähestyminen kapasiteetin ylärajaa (voi) johtaa ajautumiseen kohti ääripistettä.

Logistisella mallilla on olemassa myös Allee-efektillä täydennetty muoto

\begin{equation}\frac{dN}{dt}=-rN(1-\frac{N}{A})(1-\frac{N}{K}), \end{equation}
missä populaation kasvuvauhti on negatiivinen kun 0<A<N sekä positiivinen kun N<A<K, oletuksella että on aina voimassa A<K. Ensimmäisessä tapaksessa Allee-efektiä sanotaan "vahvaksi", jälkimmäisessä "heikoksi".  Tällöin populaation koolla on siis vaikutusta kasvunopeuteen, mutta vaikutuksen suuruus ja suunta riippuvat muuttujan A arvoista. Allee-efektillä varustetun logistisen kasvumallin kuvaajat ovat vahvalla Allee-efektillä kun K=102 ja A=40:

sekä heikolla kun A = 10 ja N = 11-100:

Allee-efektillä täydennetyn mallin stabiilisuusominaisuudet seuraavat pitkälti tavallisen logistisen mallin vastaavia, joten niitä ei tässä erikseen käydä läpi. Aiheesta on syvällisemmin kirjoittanut Lv, Chen ja Chen artikkelissaan.

Wallace ja Tsoularis muotoilivat esseessään seuraavan yleisen muodon logistisen kasvun funktiolle:

\begin{equation} \frac{dN}{dt}=rN^{\alpha}(1-(\frac{N}{K})^{\beta})^{\gamma} .\end{equation}

Huomataan, että kun alfa = 1 ja gamma = 0, saadaan eksponentiaalinen kasvumalli. Kun taas alfa = 1, beta=1 ja gamma=1, saadaan tavallinen logistinen kasvumalli. Muuttujien eri arvoilla voidaan edelleen muotoilla suoraan ainakin Gompertzin malli, Blumbergin malli, Hyper-Gompertz... käytännössä kaikki kuvaavat samanlaista käyrää, mutta pienin erilaisin painotuksin. Mainittuja malleja ei tässä käsitellä enempää.

Käytetyt R-koodit:

eksponentiaalifunktio <- function(r,t){
    funktio <- r*exp(t)
    return(funktio)
}

plot(x=1:7,y=eksponentiaalifunktio(0.1,1:7),type = "l", col="red",xlab="aika",ylab = "populaation koko")

logistinen_funktio <- function(r,N,K){
    funktio <- r*N*(1-N/K)
    return(funktio)
}

plot(x=0:100,y=logistinen_funktio(0.5,0:100,100),type = "l", col="red",xlab="N",ylab = "F(N)")

logistinen_funktio1 <- function(p_0,K,t,r){
    funktio <- p_0*K*exp(r*t)/((k-p_0)+p_0*exp(r*t))
    return(funktio)
} 

plot(x=0:30,y=logistinen_funktio1(1,100,0:30,0.5),type = "l", col="red",xlab="aika",ylab = "populaation koko")

allee_funktio <- function(r,N,K,A){
    funktio <- -r*N*(1-N/A)*(1-N/K)
    return(funktio)
}

plot(x=0:100,y=allee_funktio(0.5,0:100,102,40),type = "l", col="red",xlab="N",ylab = "F(N)")

lines(x, y2, type = "l", col = "black")   

plot(x=0:89,y=allee_funktio(0.5,11:100,102,10),type = "l", col="red",xlab="N",ylab = "F(N)")

lines(x, y2, type = "l", col = "black")  

Lotka-Volterra ja kahden lajin logistinen kilpailumalli

Diclaimer: Bloggerin matematiikkatyökalu on järjettömän buginen ja epäluotettava. Siksi olenkin tekstissäni pyrkinyt mahdollisuuksien mukaan välttämään ylimääräisten symbolien käyttöä.

Eräänä päivänä joskus ensimmäisen maailmansodan melskeiden jälkeisinä vuosina italialainen matemaatikko Vito Volterra sai tyttärensä ystävän vaikutuksesta tehtäväkseen selittää kalakantojen mystistä vaihtelua. Sodan vuoksi kalastaminen oli jäänyt maassa pitkäksi aikaa vähemmälle, ja Italian vesissä kaiken järjen mukaan eli enemmän saaliskalaa kuin aikoihin. Käytännössä kuitenkin huomattiin kalastajien haaveihin tarttuvan enemmän petoja kuin ravintoketjun pohjilla olevaa saalista. Volterra julkaisi tuloksensa vuonna 1926. Tämän tietämättä samaiseen tulokseen oli kuitenkin tullut Alfred J. Lotka soveltaessaan kemiallisiin reaktioihin kehittämiään funktioita biologisiin kohteisiin. Näin mallin nimeksi siis vakiintui Lotka-Volterran peto-saalis -malli. Kyseessä on yksi biologisen matematiikan perustuloksista.

Lotka-Volterran malli perustuu seuraaville yksinkertaisille oletuksille:

1. Saaliseläimet löytävät aina tarpeeksi ruokaa.

2. Petojen ravinto riippuu täysin saaliseläimistä.

3. Populaation muutosnopeus riippuu populaation koosta.

4. Tarkasteluaikana ympäristö ei muutu.

5. Petoeläimillä on loputon ruokahalu.

6. Kumpaakin populaatiota on vain yhtä "tyyppiä", eikä lajien yksilöiden välillä ole merkittävää eroa.

Käytännössä mallin oletukset eivät ole kovinkaan realistisia, minkä vuoksi Lotka-Volterrassa on perinteisesti ollut paljon kehittämisen ja soveltamisen varaa. Lotka-Volterran peto-saalismallia sekä logistista kasvumallia on edelleen jatkokehitelty ja yhdistelty spesifimpiä käyttökohteita silmällä pitäen, kuten Arditi-Ginzburgin sekä Rosenzweig-MacArthurin mallit, mutta nämä ansaitsevat omat blogipostauksensa joskus toiste.Tässä on käsitelty lähinnä yksinkertaiseen Lotka-Volterraan sekä kahden samasta ekologisesta lokerosta kilpailevan lajin malliin.

Lotka-Volterra -peto-saalis-malli

Lotka-Volterran peto-saalis-malliksi kutsutaan seuraavaa yhtälöparia:

\begin{equation} \frac{dx}{dt} = \alpha x - \beta xy \end{equation}

\begin{equation} \frac{dy}{dt} = \delta xy - \gamma y, \end{equation}

missä x on saaliseläinten esiintymistiheys, y saalistajien esiintymistiheys, alfa saaliseläimen maksimitiheys, beta saalistajan vaikutus saalistettavan lisääntymisnopeuteen, delta saaliseläimen vaikutus saalistajan lisääntymisnopeuteen sekä gamma saalistajan kuolemisnopeuteen.

Otetaan seuraavaksi esimerkki, jossa alfa on 1.1, beta 0.4, delta 0.1 ja gamma 0.4. Tällöin saadaan seuraavat ajasta riippuvat kuvaajat:

Kuvassa ylempänä saaliseläimen, tässä tapauksessa kalan, dynamiikka suhteutettuna saalistajan, eli karhun, vastaavaan suhteessa aikaan. Kuvasta huomataan petoeläimen nopean kasvun johtavan saaliseläimen ja pienellä viiveellä myös saalistajan itsensä kannan romahtamiseen. Saaliseläinten kannan kasvaessa seuraa saalistajan kanta perässä. 

Esimerkki ja sen visualisointiin käytetty python-koodi löytyvät Mike Saint Antoinen Youtube-kanavalta, josta löytyy enemmänkin ekologiseen mallintamiseen paneutuvaa materiaalia käyttökelpoisine ja selkeine koodeineen. Näitä tullaan soveltamaan enemmänkin jatkossa sikäli kuin mahdollista.

Jatketaan kuitenkin mallimme ominaisuuksien tutkimista. Haetaan ensin Lotka-Volterran tasapainopisteet, eli kun derivaatta on nolla, kirjoitetaan

\begin{equation} 0 = \alpha x - \beta xy \end{equation}

\begin{equation} 0 = \delta xy - \gamma y, \end{equation}

mistä nähdään, mallilla on kaksi tasapainopistetta, ensimmäinen kun x=0 ja y=0 sekä toinen, kun y = alfa/beta ja x = gamma/delta. Ottamalla mallista Jacobin matriisi

\begin{bmatrix}\alpha-\beta y & -\beta x \\ \delta y & \delta x - \gamma \end{bmatrix}

huomataan, että pisteessä (0,0) matriisin ominaisarvot ovat alfa ja -gamma,
eli funktion nollapiste on satulapiste ja epästabiili, sillä kyseiset vakiot ovat mallissa 
oletusarvoisesti nollaa suurempia, eivätkä ominaisarvot siksi voi missään kohtaa olla 
samanmerkkiset. Näin siis lajit voivat kutistua kuinka pieneksi tahansa ja kasvaa silti entiseen
kokoonsa. Toisella tasapainopisteellä saadulla Jacobin matriisilla ei ole reaalisia ominaisarvoja. 

Kahden lajin logistinen Lotka-Volterra -kilpailumalli

Johdetaan seuraavaksi kahden lajin välinen kilpailumalli yhdistelemällä Lotka-Volterraa ja logistista kasvufunktiota

\begin{equation} \frac{dN_i}{dt} = r_i N_i (\frac{K_i-N_i}{K_i}) \end{equation}

missä N kuvaa populaation kokoa, K kapasiteettia sekä r kasvunopeutta. Yllä annettu kaava kuvaa lajin sisäisen kilpailun vaikutusta populaation kasvuun. Kun i = 1,2, saadaan kahden eri populaation toisistaan erillinen sisäinen dynamiikka. Yhdistämällä kahden populaation logistiset kasvumallit Lotka-Volterraan saadaan

\begin{equation} \frac{dN_1}{dt} = r_1 N_1 (\frac{K_1-N_1-\alpha_{12}N_2}{K_1}) \end{equation}
\begin{equation} \frac{dN_2}{dt} = r_2 N_2 (\frac{K_2-N_2-\alpha_{21}N_1}{K_2}) \end{equation}

Kaavat kuvaavat kahden populaation kilpailullista tilannetta, alfat toisen populaation vaikutusta toisen kasvuun. Edelleen voidaan laskea tasapainoehdot:

\begin{equation} N_1 = K_1 - \alpha_{12}N_2 \end{equation}

\begin{equation} N_2 = K_2 - \alpha_{21}N_1 \end{equation}

mistä huomataan, että kun kun ylemmässä kaavassa N_1 = 0, saadaan N_2 = K_1/alfa_12 sekä alemmassa N_1 = K_2/alfa_21 kun N_2 = 0. Tulosta vastaavat kuvat virtaussuuntineen ja  perusteluineen löytyvät täältä. Ottamalla seuraavaksi Jacobin matriisi, jonka diagonaalialkiot ovat muotoa

\begin{equation} J_{ii} = r_i(\frac{K_i-2N_i-\alpha_{ij}N_j}{K_i}) \end{equation}

sekä päätyalkiot muotoa

\begin{equation} J_{ij} = \frac{-r_i N_j \alpha_{ji}}{K_i} \end{equation}

Asettamalla saadut tasapainopisteet (K_1/a_21,0) ja (0,K_2/a_12), huomataan, että päätyalkioista toinen menee aina nollaan, toinen diagonaalialkioista on aina negatiivinen sekä matriisin ominaisarvot ovat samanmerkkisiä vain jos pätee K_1/a_21 > K_2 ja K_2/a_12 > K_1.

Esimerkkinä kahden kilpaleivan lajin vaikutuksesta toisiinsa saa helposti luotua kuvaajat muokkaamalla aiempaa koodia. Satunnaisilla annetuilla arvoilla saadaan susien ja karhujen välisestä kilpailusta seuraavat kuvaajat:

Kuvasta nähdään, että koska kasvuvakio r_1 > r_2, susien lisääntyminen on nopeampaa, mutta karhujen ilmaantuessa kilpailemaan elintilasta näiden populaatioiden koot tasaantuvat sinne 600-650 tietämille kummassakin. Kyseessä on kahden lajin välinen stabiili tasapainotila, sillä annetuilla arvoilla K_(1,2) = 1000 sekä alfa_(1,2) = 0.5 stabiilin tasapainopisteen ehdot täyttyvät.


Kilpailumallissa käytetty Python-koodi, muokattu aiemmasta koodista:

import numpy as np
import matplotlib
matplotlib.use('tkagg')
import matplotlib.pyplot as plt
from scipy.integrate import odeint

y0 = [10,1]
t = np.linspace(0,50,num=1000)

alpha_12 = 0.5
alpha_21 = 0.5
r_1 = 0.8
r_2 = 0.7
K_1 = 1000
K_2 = 1000

params = [alpha_12, alpha_21,r_1,r_2,K_1,K_2]
    x = variables[0]
    y = variables[1]
    
    alpha_12 = params[0]
    alpha_21 = params[1]
    r_1 = params[2]
    r_2 = params[3]
    K_1 = params[4]
    K_2 = params[5]
    
    dxdt = r_1 * x * (K_1 - x - alpha_12 * y)/K_1
    dydt = r_2 * y * (K_2 - y - alpha_21 * x)/K_2
    
    return([dxdt,dydt])

y = odeint(sim, y0, t, args=(params,))

f,(ax1,ax2) = plt.subplots(2)

line1, = ax1.plot(t,y[:,0], color="b")

line2, = ax2.plot(t,y[:,1], color="r")

ax1.set_ylabel("Wolves (hundreds)")
ax2.set_ylabel("Bears (hundreds)")
ax2.set_xlabel("Time")

plt.show()