Evoluční model pesimizace prostředí v Maple Jiří Kalina Masarykova univerzita, Institut biostatistiky a analýz, Kamenice 126/3, Brno, 625 00 kalina@mail.muni.cz Abstrakt Příspěvek se zabývá popisem a diskuzí populačního modelu se zahrnutím náhodné evoluční složky, založeného na jednoduchém modelu chemostatu popsaném v [2]. Model je rozšířen na neomezený počet populací a je z něj vyňata omezující podmínka rovnosti mortality a podílu nevyužitého substrátu podléhajícího zkáze, což umožňuje jeho aplikaci i na složitější společenstva. Jádrem příspěvku je implementace modelu v prostředí Maple, ve kterém je pomocí numerického řešiče diferenciálních rovnic dsolve spočten vývoj velikosti populací a na výsledných grafech je demonstrována oprávněnost předem vyslovených předpokladů. Abstract This paper deals with the description and discussion of the population model including a random evolutionary component, based on a simple model of chemostat described in [2]. The model is extended to an unlimited number of populations and a condition of equality of mortality and the proportion of unused perishable substrate, which allows its application to more complex systems. The core contribution is an implementation of the model in a Maple environment with an use of a numerical differential equation solver "dsolve". The evolution of population sizes is counted and plots of results demonstrate the legitimacy of former assumptions. Klíčová slova Populační model, pesimizace. Keywords Population model, pesimization. 1 Základní vlastnosti modelu Předpokládejme uzavřený systém složený z prostředí a navzájem přímo neinteragujících populací, které ho obývají, ∈ . Jedince náležející k jedné populaci budeme souhrnně nazývat druhem. Úživnost prostředí, tj. veličinu, která charakterizuje podmínky pro přežití a reprodukci jedinců, označme a považujme ji zjednodušeně za množství potravy (substrátu, resp. energie), které je k dispozici společně pro všechny populace v prostředí. budeme považovat za funkci času . Přírůstek substrátu v prostředí za jednotku času (růst, resp. transport energie do systému) považujme za konstantní a označme jej . Vzhledem k předpokladu, že populace v modelu spolu přímo neinteragují, omezuje se jejich vzájemné ovlivnění na soupeření o dostupné množství potravy, jde tedy o ryze kompetitivní vztah [1]. Vzhledem k jednoduchosti modelu budeme druhy odlišovat pouze schopností jejich zástupců spotřebovávat společnou potravu a takto přijatý substrát využít k vlastnímu přežití a reprodukci. Tento proces popíšeme dvojicí veličin a , kde představuje množství substrátu přijaté jedním jedincem -té populace za jednotku času. Tuto veličinu budeme dále považovat za funkci množství dostupného substrátu . Naproti tomu hodnota označuje efektivitu přeměny substrátu na živé jedince -tého druhu jako poměr spotřebovaného substrátu vůči počtu nově narozených jedinců. V první variantě modelu budeme hodnoty a považovat za vzájemně nezávislé. Dále budeme předpokládat konstantní úmrtnost , stejnou pro všechny populace, jako podíl jedinců daného druhu uhynulých za jednotku času a nakonec zbytkový podíl nespotřebovaného substrátu , který podlehne zkáze, aniž by byl spotřebován některou z uvažovaných populací. Na počátku modelovaného období budeme předpokládat jednu (nebo větší konstantní počet) populaci v prostředí, která spotřebovává dostupný substrát pro vlastní reprodukci. Snadno lze ukázat, že po určité době dojde k ustálení počtu jedinců populace (uvědomme si, že jde o silně zjednodušený idealizovaný model) na hodnotě odpovídající úživnosti prostředí, která se v dalším průběhu času již nebude měnit. Je nutno poznamenat, že pro jednoduchost model opomíjí veškeré komplikace, které s sebou přináší rozmanitost pohlavního rozmnožování jako je samo pohlaví, hledání partnera či nestejný počet potomků každého jedince. Budeme tedy nadále přepokládat, že dostatečným požadavkem pro další rozmnožování populace je přítomnost alespoň jednoho jedince daného druhu v systému, že potomkem jedince určitého druhu je (až na výjimky) jedinec téhož druhu a že křížení mezi jedinci různých druhů není možné. Situace se nicméně začne komplikovat ve chvíli, kdy se v prostředí objeví zástupce nového druhu. V modelovém uzavřeném systému není možný průnik takového jedince zvenčí, proto jej budeme považovat za mutanta některé z již existujících populací, který se bude odlišovat parametry a , poruší tedy předpoklad z předcházejícího odstavce. Výskyt mutantů bude relativně řídký a bude se odehrávat pouze v určitých okamžicích. Takový jedinec se v našem modelu může, bude-li úspěšný ve vzájemné kompetici, stát zakladatelem populace nového druhu, v opačném případě může ovšem nová populace záhy vyhynout a uvolnit tak prostředky pro zástupce jiných druhů. Časovou změnu množství substrátu dostupného v systému je na základě uvedených předpokladů možné vyjádřit následujícím vztahem: ∙ ∑ ∙ (1) kde je dostupné množství substrátu v čase , je přírůstek substrátu za jednotku času, je podíl nevyužitého substrátu za jednotku času, je počet populací přítomných v systému, je spotřeba substrátu pro jednoho jedince -tého druhu za jednotku času a je počet jedinců -té populace. Popišme nyní tvar závislosti množství substrátu spotřebovaného jedincem -té populace za jednotku času. Opět zjednodušeně předpokládejme, že dostupnost substrátu v prostředí je stejná pro všechny jedince všech populací a její vyhledávání proto nehraje v modelu žádnou roli. V případě, že bude množství dostupného substrátu nulové, 0, bude přirozeně rovněž 0. Jakmile začne v prostředí narůstat množství dostupné potravy, poroste přímo úměrně také množství zkonzumované jedinci všech populací. Vzhledem k fyziologickým možnostem organizmů je však maximální množství potravy přijaté za jednotku času shora omezené, proto se bude s narůstajícím množstvím substrátu v prostředí křivka funkce oddělovat od přímky přímé úměry, až plynule konverguje ke konstantní funkci. Uvedený předpoklad dokonale splňuje v biologii a příbuzných oborech hojně využívaná logistická křivka následujícího tvaru: ∙ ∙ ∙ (2) kde parametr určuje strmost počáteční části křivky a rovněž rychlost dosažení následující asymptotické fáze, zatímco parametr její výšku. V první, jednodušší, variantě modelu budeme předpokládat jednak, že je totožné pro všechny populace (což není sám o sobě předpoklad, který by výrazně ovlivnil chování modelu), současně ale také, že je konstantní a nezávislé na ostatních parametrech modelu. Rovněž parametr budeme v první variantě modelu považovat za konstantu. Přejděme nyní k rovnici popisující vývoj počtu jedinců -té populace v čase. Změna počtu jedinců, bude vyjádřena jako množství substrátu zkonzumovaného jedinci populace, vynásobeného efektivitou přeměny substrátu na živé jedince, minus počet uhynulých jedinců: ∙ ∙ ∙ (3) kde je počet jedinců populace v čase , je spotřeba substrátu pro jednoho jedince -tého druhu za jednotku času, je efektivita přeměny substrátu na živou hmotu -tou populací a je relativní úmrtnost, pro jednoduchost opět totožná pro všechny populace. Uvedený postup nám poskytuje soustavu 1 diferenciálních rovnic, jejímž řešením lze získat informaci o vývoji počtu jedinců jednotlivých populací v čase. Vzhledem k tomu, že v první variantě modelu uvažujeme totožné pro všechny populace v systému, je řešení závislé především na hodnotách pro jednotlivé populace. V případě, že se v systému nevyskytne mutant, lze ukázat, že v dlouhodobém horizontu se v závislosti na hodnotách konstantních parametrů , , a a na počátečních podmínkách 0 a 0 v systému ustálí stav, ve kterém je množství dostupného substrátu v prostředí konstantní a navíc počet jedinců nejvýše jedné populace nekonverguje k nule. Platí totiž: ∑ ∑ ∑ (4) a tedy po integraci podle potom: ∑ ∙ ∑ ∙ ∙ (5) a po úpravě: ∑ ∑ ∙ ∙ (6) za předpokladu, že pro → ∞ konvergují 2. a 3. sčítanec k nule, dostáváme výsledek: lim → ∑ (7) Tedy, pokud bude úživnost prostředí dostatečná, zvítězí v kompetici o omezený zdroj potravy v uzavřeném systému pouze jediný druh (viz (8) a (9)), zatímco ostatní postupně vyhynou. Lze ukázat (viz [2]), že půjde o druh s nejvyšší hodnotou , tj. ten, který je schopen nejefektivněji zpracovávat stravu ve prospěch přírůstku vlastních jedinců. Zkomplikujme nyní situaci možností, že se v určitém čase v systému objeví mutant, tj. jedinec, jehož se bude lišit od jeho rodiče a půjde tak vlastně podle našich kritérií o zakladatele nového druhu. V modelu předpokládáme dobu mezi výskytem dvou mutantů náhodnou s rovnoměrným rozdělením pravděpodobnosti, což sice neodpovídá reálnému prostředí, nicméně zřejmě samotná perioda výskytu mutantů nemá na stav systému z dlouhodobého hlediska vliv. Hodnotu mutovaného potomka pak budeme uvažovat v intervalu 〈 ; ∙ 〉, kde odpovídá hodnotě efektivity přeměny substrátu na živou hmotu jeho rodiče. Rozdělení pravděpodobnosti využijeme pro jednoduchost opět rovnoměrné. V modelu budeme uvažovat stejnou pravděpodobnost zplození mutanta pro všechny žijící jedince, tj. pravděpodobnost vzniku mutanta -té populace bude přímo úměrná velikosti této populace. 2 Implementace modelu v prostředí Maple Základem modelu v softwarovém prostředí Maple 16 je cyklus, jehož jedna iterace odpovídá období mezi výskytem dvou mutantů. Časová délka období reprezentovaného iterací je proměnlivá. Jednotlivé populace jsou v modelu reprezentovány řádky populační matice , jejíž první sloupec udává pro každou populaci efektivitu , druhý sloupec pak odpovídá počtu jedinců dané populace. celkemcyklu := 250; a := 2.5; b := 0.8; d := 12; m := 0.1; z := 0.1; cas := 0 q := proc (x) options operator, arrow; a*x/(1+a*b*x) end proc Před spuštěním cyklu je vygenerována dvouřádková matice se dvěma náhodnými populacemi, které vstupují do prvního období. Délka období je volena náhodně z uživatelem zadaného intervalu. Na začátku každého období je provedena kontrola všech řádků matice a „neúspěšné“ populace s počtem jedinců nižším než 1 (tj. populace nesplňující zadanou podmínku pro další rozmnožování)1 jsou z modelu odstraněny. Následně probíhá vygenerování náhodného mutanta, který je zapsán jako nový řádek do matice (tj. populace s jediným jedincem). # Matice se dvema nahodnymi populacemi; h1 := rand(1 .. 10); h2 := (1/100)*rand(1 .. 100); Matice := Matrix(2, 2, [[h2(), h1()], [h2(), h1()]]); # Hlavni cyklus; for cyklus from 1 to celkemcyklu do # Cyklus pro odstraneni populaci s mene nez 1 jedincem; v:=0: for i from 1 by 1 to RowDimension(Matice) do if (Matice[i-v,2]<1) then Matice:=DeleteRow(Matice,i-v); v:=v+1 end if: end do: # Cyklus pro vznik jednoho jedince s nahodnou mutaci; jedincu:=add(Matice[k,2],k=1..RowDimension(Matice)); h:=(rand(1..10000* floor(jedincu)))/(10000): budemutant:=h(): soucet:=0: j:=0: for j from 1 by 1 while soucet: V dalším kroku je pomocí vnořeného cyklu vytvořena soustava 1 diferenciálních rovnic postupně pro jednotlivé populace z matice a navíc pro celkové dostupné množství substrátu v systému, z velikosti populací a proměnné , udávající množství dostupného substrátu po skončení předešlého období jsou vygenerovány počáteční podmínky pro řešení a následně je soustava řešena pomocí numerické varianty příkazu . Řešení je z důvodu omezení časové náročnosti hledáno pouze pro období mezi výskyty mutantů, po jehož uplynutí se změní počáteční podmínky. # Vytvoreni soustavy diferencialnich rovnic; rovnice[0]:=diff(S(t),t)=d-z*S(t)- q(S(t))*add(P[i](t),i=1..RowDimension(Matice)): for i from 1 by 1 to RowDimension(Matice) do rovnice[i]:=diff(P[i](t),t)=-m*P[i](t)+q(S(t))*P[i](t)*Matice[i,1]: end do: soustava:=seq(rovnice[i],i=0..RowDimension(Matice)): # Vytvoreni pocatecnich podminek; if cas=0 then podminka[0]:=S(0)=1: else podminka[0]:=S(cas)=Y: end if: for i from 1 by 1 to RowDimension(Matice) do podminka[i]:=P[i](cas)=Matice[i,2]: end do: podminky:=seq(podminka[i],i=0..RowDimension(Matice)): # Cas do pristi mutace; minule:=cas: h:=rand(1..40): cas:=cas+h(1..40): # Reseni a vykresleni do grafu; reseni:=dsolve({soustava,podminky},numeric, maxfun=10000000): krivky:=seq([t,P[i](t)],i=1..RowDimension(Matice)): graf[cyklus]:=odeplot(reseni,([krivky],minule..cas)): assign(reseni(cas)): # Zapis novych hodnot do matice; for i from 1 by 1 to RowDimension(Matice) do Matice[i,2]:=P[i](t): unassign('P[i](t)'): end do: Y:=S(t): unassign('t'): unassign('P'): unassign('S'): end do: Vzhledem k tomu, že konverguje v období mezi dvěma mutacemi (jak bylo výše ukázáno) ke konstantě, konverguje podle definice ke konstantní hodnotě rovněž funkce příjmu potravy a tedy můžeme provést následující výpočet: →∞ ∙ ∙ →∞ ∙ (8) Tedy zřejmě velikost populace v rámci období mezi mutacemi konverguje k monotónnímu chování. Předpokládejme nyní: →∞ . ⇒ ∙ ∙ →∞ 0 →∞ 0 ⋁ →∞ (9) Protože pro libovolné populace , : platí , konverguje v období mezi dvěma mutacemi velikost všech populací, případně vyjma populace s nejvyšším , k nule. Příkazem display(seq(graf[i], i = 1 .. celkemcyklu)) si nyní můžeme zobrazit přes sebe graf vývoje populací ve všech obdobích (obr. 1) a na první pohled je zřejmé, že skutečně v každém období mezi dvěma mutacemi je jediná populace, jejíž velikost roste „na úkor“ ostatních populací. Ve variantě bez mutací by se po určité době velikost této populace ustálila na hodnotě úživnosti prostředí, zatímco velikost konkurenčních populací by klesla k nule. Náhodně se vyskytující mutace však do systému vnáší dynamické chování, neboť v případě, kdy se objeví mutant s vyšším, než bylo dosavadní nejvyšší , rozšíří se tato -tá populace a postupně zapříčiní vyhynutí všech populací včetně -té, pokud není ještě dříve vystřídána jinou populací, vzešlou z mutanta s vyšším . Vzhledem k tomu, že v této variantě modelu není výše efektivity přeměny potravy na živou hmotu ničím omezena, dochází k tomu, že se neustále objevují mutanti s vyšší hodnotou , která může v rozporu s realitou překročit dokonce hodnotu 1 a růst neomezeně do nekonečna. Velikost momentální největší populace se tak bude v modelu pohybovat v čase po exponenciální křivce donekonečna. Takový proces je samozřejmě zcela nerealistický a proto je třeba model upravit. Zvyšování efektivity přeměny potravy na živou hmotu vlastního druhu se v modelu ukazuje jako správná strategie vedoucí k přežití. Je ovšem patrné, že s rostoucím vítězné populace dochází k poklesu substrátu v prostředí, neboť vítězná populace ostatní jednoduše „vyhladoví“. V našem nerealistickém případě z dlouhodobého hlediska dochází až ke konvergenci k nule. Tento pokles úživnosti a obecně evoluční zhoršení podmínek prostředí na ty nejhorší možné se nazývá pesimizací prostředí. Obr. 1: Exponenciální růst populací v modelu ve variantě bez omezení efektivity přeměny substrátu na živou hmotu. 3 Varianta s omezením efektivity Fyzikální podstata živých tvorů neumožňuje růst efektivity nad hodnotu 1 (ve skutečnosti je číslo vzhledem k chemickým procesům v organizmu samozřejmě ještě nižší), vlastní fyziologie organizmů však znamená ještě výraznější pokles efektivity, neboť v reálném prostředí nemůže věnovat žádný jedinec maximum prostředků na hospodárné nakládání se získanou energií (tj. vlastní růst a reprodukci), ale musí rovněž zajistit řadu dalších procesů, mezi které kromě jiných strategických činností patří samotný příjem potravy. V předchozí zjednodušené variantě modelu jsme předpokládali závislost množství přijaté potravy pouze na množství substrátu v prostředí, zbývající parametry a jsme považovali za konstantní pro všechny populace. Definujme nyní vztah mezi efektivitou přeměny substrátu na živou hmotu a množstvím potravy, kterou je jedinec schopen přijmout za jednotku času , který bude reprezentovat strategii druhu ve smyslu rozhodování mezi investicí energie do většího množství získané potravy nebo lepším využitím potravy pro vlastní růst a reprodukci. Nahradíme tedy v rovnici logistické křivky (2) konstantu kvadratickou funkcí takto: (10) tj. s rostoucí efektivitou 0 bude injektivně klesat množství substrátu , které je jedinec daného druhu schopen zkonzumovat. Tato vazba zamezí nekontrolovanému růstu a přiblíží model realitě. Optimální strategií pak bude dosáhnout takového , aby byla hodnota ∙ q (viz (3)) maximální, tj: ∙ ∙ ∙ ∙ ∙ (11) a po derivaci podle : ∙ ∙ 1 ∙ ∙ 2 ∙ ∙ 1 ∙ ∙ 2 ∙ ∙ 2 . (12) Po úpravě, protože jmenovatel je nenulový, získáváme vztah pro optimální : 1 ∙ . (13) Ve zdrojovém kódu Maple stačí změnit definici funkce na funkci dvou proměnných a a upravit cyklus určující tvar diferenciálních rovnic: q := proc (x, y) options operator, arrow; a*x/(1+a*(y^2+c)*x) end proc # Vytvoreni soustavy diferencialnich rovnic; rovnice[0] := diff(S(t), t) = d-z*S(t)-add(q(S(t), Matice[i, 1])*P[i](t), i = 1 .. RowDimension(Matice)); for i to RowDimension(Matice) do rovnice[i] := diff(P[i](t), t) = -m*P[i](t)+q(S(t), Matice[i, 1])*P[i](t)*Matice[i, 1]: end do; soustava := seq(rovnice[i], i = 0 .. RowDimension(Matice)) Podle očekávání je po vykreslení grafu jasně patrný nárůst celkového momentálně nejúspěšnějších populací, který však zpomaluje až ke konstantní hodnotě odpovídající úživnosti prostředí pro populaci s hodnotou blízkou podle rovnice (13). Z dlouhodobého hlediska by za neměnných podmínek v systému opravu po řadě mutací převládla jediná populace s . Protože je nejvyšší ze všech dlouhodobě možných efektivit přeměny substrátu na živou hmotu, dojde postupem času v prostředí opět k poklesu dostupného množství potravy na nejnižší možnou mez – vzájemný souboj populací o zdroje tedy stejně jako v předchozím případě nastoluje ty nejhorší podmínky pro přežití, tj. dochází k evoluční pesimizaci. Obr. 2: Omezený růst populací v modelu s omezením efektivity přeměny substrátu na živou hmotu. Po dosažení hodnoty se růst zastavuje na hodnotě úživnosti prostředí. 4 Závěr Ve variantě modelu, kdy není efektivita přeměny substrátu na živou hmotu shora omezena, dochází vlivem prostředí za vzniku náhodných mutací k postupnému zvyšování přibližně exponenciálním tempem, přičemž každá populace po určité době vyhyne, protože úživnost prostředí bude snížena populací s vyšší efektivitou pod přípustnou mez. Pokud je do modelu vložena zpětná vazba omezující růst efektivity jeho klesající úměrností vůči množství zkonzumovaného substrátu , dochází postupnými mutacemi k přibližování maximální efektivity k mezní hodnotě úživnosti prostředí ̂, která je daná fyziologickými vlastnostmi organizmů. V tomto případě platí odvozené vztahy dokazující, že nejvýše jedna populace v dlouhodobém měřítku nevyhyne a obsadí celý systém. Pro tuto populaci platí podmínka ̂, tj. není ohrožena žádným mutantem. Vlivem evolučních změn vedoucích k postupnému hynutí populací s nižší hodnotou než mají jejich konkurenční mutanti, dochází k tzv. pesimizaci prostředí, tj. snížení úživnosti na nejmenší přístupnou mez. V krajním případě může tento jev vést až k úplnému vyhynutí všech populací. 5 Literatura [1] Hřebíček, J., Pospíšil, Z., Urbánek J.: Úvod do matematického modelování s využitím Maple. Akademické nakladatelství CERM, s.r.o., Brno, 2010. [2] Diekmann, O. A beginner’s guide to adaptive dynamics. Mathematical modelling of population dynamics, Banach Center Publications, vol. 63. Institute of Mathematics Polish Academy of Sciences, Warzsawa, 2004.