Učební text k předmětu Matematické modely v biologii Úvod Tento text slouží k výuce předmětu M81B0 Matematické modely v biologii na Přírodovědecké fakultě MU a je zamýšlen především, avšak ne výlučně, jako studijní opora pro studenty oboru Matematická biologie. Obsah Enzymová kinetika, 1. část 3 1 Kinetika chemických reakcí....................... 3 2 Základní princip enzymové reakce................... 7 3 Quasi-steady-state aproximace..................... 11 Enzymová kinetika, 2. část 19 1 Časová škála................................ 19 2 Rychlost reakce .............................. 22 3 Smyslové receptorové systémy..................... 26 Difúze 30 1 Mikroskopický pohled a náhodná procházka ............ 30 2 Makroskopický pohled a Fickův zákon................ 37 Řešení rovnice difúze, 1. část 45 1 Úvod..................................... 45 2 Nekonečná trubice ............................ 47 3 Nekonečná trubice s podélným zdrojem............... 52 4 Trubice s nulovou koncentrací na okraji............... 56 5 Trubice s podélným zdrojem a nulovou koncentrací na okraji . . 59 Lánský, Pokora: 8. září 2014 1 Řešení rovnice difúze, 2. část 64 1 Trubice s nulovým tokem na okraji .................. 64 2 Trubice s podélným zdrojem a nulovým tokem na okraji..... 67 3 Trubice s trvalým zdrojem na okraji.................. 68 4 Trubice s bodovým zdrojem a nulovou koncentrací na okrajích . 70 5 Radiální difúze z mikropipety...................... 74 Vybraná témata z teorie informace 78 1 Úvod do teorie informace........................ 78 2 Axiomatické odvození Shannonovy entropie ............ 81 3 Kvantizace a entropie spojité náhodné veličiny........... 85 Matematické modely neuronu 92 1 Logický neuron .............................. 92 2 Logický neuron a logické funkce.................... 98 3 Model reálného neuronu......................... 99 Použité zdroje a literatura k dalšímu studiu 106 2 Enzymová kinetika, 1. část Předpoklady Pro studium této výukové jednotky je potřebná znalost • vyšetřování průběhu funkce, • lineárních diferenciálních rovnic a metod jejich řešení. Výstupy z učení Studenti se naučí • modelovat kinetiku chemických reakcí pomocí diferenciálních rovnic, • využívat quasi-steady-state aproximace pro hodnocení dynamiky reakcí. 1 Kinetika chemických reakcí Pro kinetický popis biochemických reakcí se jako vhodné (nezbytné) ukázalo matematické modelování. Obzvláště zajímavá je kinetika enzymových reakcí, jejímž principům se v této kapitole budeme věnovat. Vycházíme přitom především z monografií Segel (1984), Murray (2002, kapitola 11) a Keener & Sneyd (2009, kapitola 1), v nichž lze nalézt další modely. Nejdříve si vysvětlíme užívaná schématická značení a způsob sestavování rovnic matematického modelu. Uvažujeme chemickou reakci, v níž molekula látky A reaguje s molekulou látky B za vzniku molekuly látky C. Naopak, molekula C se za určitých podmínek může rozpadnout na molekuly složek A a B. Pro popis množství látek v čase se používá veličina koncentrace. Tu lze uvažovat v několika variantách, např. jako hmotnost látky v jednotkovém objemu, jako molární množství látky v jednotkovém objemu, anebo jako počet molekul látky v jednotkovém objemu. Standardně užívanou koncentrací pro popis biochemických reakcí je molární množství látky v jednotkovém objemu. Ačkoliv odpovídající jednotkou v soustavě SI je mol m-3, v praxi obvykle používanou jednotkou je mol dm-3, někdy zkráceně označovanou symbolem M. Pro matematický popis dynamiky reakce zavedeme nezáporné funkce A(t), B(t), C(t) označující koncentraci látek A, B, C v čase ř > 0. 3 Tradičně se předpokládá, že rychlost tvorby látky C je úměrná koncentracím látek A i B, tento princip je nazýván law ofmass action. Tento předpoklad znamená, že pokud koncentrace látky B vzroste na dvojnásobek (resp. trojnásobek), rychlost tvorby látky C bude také dvojnásobná (resp. trojnásobná). Pokud koncentrace látky B vzroste dvakrát a koncentrace látky A třikrát, rychlost tvorby látky C vzroste dokonce šestkrát. Takový předpoklad bude splněn při nízkých koncentracích, kdy se každá molekula pohybuje nezávisle na ostatních molekulách. V této situaci znamená zvýšení koncentrace látky B úměrný nárůst počtu „srážek" molekul látek A a B, které mohou vést ke vzniku látky C a tedy odpovídajícímu nárůstu koncentrace látky C. Dále předpokládáme, že molekuly C se chovají vzájemně nezávisle, tedy že molekula C se určitou danou a pro všechny molekuly C stejnou pravděpodobností za jednotkový čas rozpadne na molekuly A a B. Při obou předpokladech o tvorbě a rozpadu látky C zavedeme konstanty úměrnosti fcl5 udávající podíl A-B-kolizí vedoucích ke vzniku C za jednotku času, a k_i udávající průměrný počet rozpadů molekul C za jednotku času. Tyto konstanty úměrnosti jsou nazývány rychlostními konstantami (rate constants) a jejich podíl je někdy označován jako tzv. disociační konstanta Kd = -r1. (D h Rychlostní konstanty přitom nemají stejné jednotky. Typicky, jednotkou rychlostních konstant popisujících analýzu je převrácená hodnota jednotky času, zatímco jednotkou rychlostních konstant spojených se syntézou je převrácená hodnota součinu jednotky času a koncentrace. Disociační konstanta proto není bezrozměrnou veličinou, ale má rozměr koncentrace. Uvažovaná reakce se graficky znázorňuje pomocí tzv. kinetického schématu A + B C k-i ' (2) kde znak + naznačuje slučování (syntézu), resp. rozklad (analýzu), látek a dvojitá šipka označuje reverzibilní reakci. Ireverzibilní reakce se označují jen šipkou jednoduchou. Veličina uvedená u každé šipky je pak rychlostní konstantou pro příslušnou reakci, příp. zvlášť pro každý směr reverzibilní reakce. 4 Pro popis dynamiky systému potřebujeme matematické vyjádření změn koncentrací v závislosti na časové změně. Koncentrace A(t), B(t), C(ř) proto často považujeme za spojité a diferencovatelné funkce času ř a dynamiku reakce matematicky vyjadřujeme pomocí (systému) diferenciálních rovnic. Z interpretačních důvodů může být názornější začít sestavením diferenčních rovnic, z nichž pak limitním přechodem obdržíme rovnice diferenciální. Tento postup pro názornost použijeme i my při sestavování první rovnice. Sledujeme schéma (2) a popíšeme změnu koncentrace, AC(ř), látky C za krátký časový interval [t; t + A ŕ] délky Ař > 0. Během uvažované krátké doby Ař koncentrace látky C naroste o fc1A(t)B(t) Ař sloučením látek A a B. Využíváme zde law of mass action, konstanta fcj udává, kolik nových molekul se vytvoří při jednotkových koncentracích látek A a B za dobu jednotkové délky, za interval délky Ař se tedy sloučí fcj Ař molekul. Reakce je však reverzibilní, současně dochází k rozpadu látky C, tedy ke snižování koncentrace C (ŕ). Součin k_i Ař udává počet rozpadů za interval délky Ař při jednotkové koncentraci látky C. Současným použitím obou směrů reakce sestavíme diferenční rovnici AC(ř) = C(ř + Ař)-C(ř) = fc1A(ř)B(ř)Ař-fc_1C(ř) Ař. (3) Tuto rovnici vydělíme výrazem Ař a limitním přechodem dospějeme k derivaci lim -= —-—, (4) At^o+ Ař dř čímž (3) přejde v diferenciální rovnici pro koncentraci látky C, dC(ť dř fc1A(ř)B(ř)-/c_1C(ř). (5) Tuto rovnici využijeme pro určení jednotek konstant fcj a k_i. Pokud koncent-raceA(ř),B(ř), C(ř) látek A, B, C uvažujeme v jednotkách mol dm-3,je rozměr levé strany (5) mol dm-3 s_1. Aby to zůstalo zachováno i pro pravou stranu této rovnice, musí být jednotka konstanty k_i rovna převrácené hodnotě jednotky času, tedy s-1, a rozměr konstanty fcj musí být rovný převrácené hodnotě koncentrace za čas, tzn. mol-1 dm3 s_1. Analogicky se sestaví diferenciální rovnice pro koncentrace látek A a B, detailní postup přenecháváme čtenáři, dA(t) dř -fc1A(ř)B(ř) + fc_1C(ř), (6) 5 (7) A(0)=A0, B(0)=B0, C(0) = C( (8) je pak počáteční úloha (5)-(8) jednoznačně zadána a lze analyticky nebo numericky hledat její řešení. Cílem toho textu však není obecné řešení uvedené úlohy, ale její aplikace na konkrétní enzymovou reakci, kterou probereme v následující části. Na závěr této obecné části si ještě všimněme, že sečtením (5) a (6), obdržíme což znamená, že součet A(t) + C(ř) je konstantní v čase. Podobně, sečtením (5) a (7), resp. (6) a (7), obdržíme další dvě rovnice. Spolu s počátečními podmínkami (8) to znamená, že platí vztahy C(t)+A(t) = C0+A0, C(t)+B(t) = C0+B0,A(t)-B(t) =A0-B0, (10) vyjadřující zákon zachování hmoty látek reagujících v uzavřené soustavě. Poznámka {law of mass action a incidence) V biochemických reakcích standardně užívaný princip law of mass action má blízko pojmu incidence, používanému např. při modelování šíření infekčních chorob nebo v analýze přežití. Vysvětlení souvislosti těchto pojmů proto věnujeme tuto poznámku. Jako tzv. contact rate A(iV) se označuje průměrný počet kontaktů objektu (např. látky A, osoby infikované chorobou) za jednotku času. Konstanta N zde označuje celkový počet zkoumaných objektů (molekul látky, osob), které jsou rozděleny do dvou kategorií, A a B. Incidencí i označujeme průměrný počet nově zformovaných objektů (např. látky C, počet přenosů infekce) všemi objekty za jednotku času. Přitom incidence, i, se počítá jako součin contact rate A(iV), počtu objektů kategorie A a pravděpodobnosti B/N,že kontakt je uskutečněn s objektem jiné kategorie (např. s látkou B, s osobou náchylnou k infekci), dř dř dř (9) N' (ID 6 V matematických modelech se používají dva typy incidencí: (i) standardní incidence (standard incidence, frequency-dependent incidence) odpovídá případu, kdy je A(iV) = fcj konstantní, nezávislá na celkovém počtu objektů N. Pro incidenci dostáváme i = —AB. (12) N (ii) mass-action incidence (density-dependent incidence), modelující případ, kdy A(iV) = kiN je závislá (lineárně) na celkovém počtu objektů N. Odpovídající incidence je rovna i = k1AB. (13) 2 Základní princip enzymové reakce Historie matematického popisu modelů enzymové kinetiky je spojena především se jmény Leonora Michaelise (1875-1949), německého biochemika, fyzikální chemika a fyzika, a Maud Mentenové (1879-1960), kanadské lékařky, působící od roku 1912 v Berlíně. Enzymy umožňují přeměnu molekul jedné látky na molekuly látky jiné. Vstupní látka, která se má transformovat, se nazývá substrát a látka, která reakcí vzniká, se nazývá produkt, budeme pro ně proto užívat zkratky S a P V typických reakcích enzymy (zkratka E) urychlují (katalyzují) proces přeměny látek tím, že se na sebe váží substrát za vzniku tzv. komplexu enzym-substrát (označení C, někdy též ES). Tento komplex C se pak může za určitých podmínek rozpadnout na volný enzym E a produkt R Tato přeměna však může také selhat v tom smyslu že komplex C se rozpadne zpět na volný enzym E a substrát S. Předpokládáme tedy, že vázání enzymu na substrát je reverzibilní, zatímco příp. rozpad komplexu na E a P je ireverzibilní. Kinetické schéma této enzymové reakce je ki k2 S + E C -► E + P ' (14) Rychlostní konstanty fcl5 k_i a k2 jsou pak po řadě spojeny se vznikem komplexu (syntéza), rozpadem komplexu na substrát (analýza) a s formováním 7 produktu z komplexu (analýza). Konstanta k2 je, stejně jako konstanta k_i spojena s analýzou komplexu C, má tedy i stejnou jednotku, s-1. Použitím principu law of mass action lze pak pro koncentrace jednotlivých látek, tedy komplexu enzym-substrát, C (ŕ), volného enzymu, E(t), substrátu, S(t), a produktu P(t), jakožto funkci času ř psát čtyři diferenciální rovnice ^^ = k1E(t)S(t)-(k_1 + k2)C(t), (15) dE(t) dt dS(t) dt -fc1£(ř)S(ř) + (fc_1 + /c2)C(ř), (16) = -k1E(t)S(t) + k_1C(t), (17) dP(t) ■k2C{t), (18) dř spolu s počátečními podmínkami C(0) = C0, E(0)=E0, S(0)=S0, P(0) = 0. (19) Na začátku reakce je množství produktu nulové, v praxi je často nulová i počáteční koncentrace komplexu (C0 = 0). Při známých konstantách fcl5 fc_l5 k2 a hodnotách počátečních koncentrací C0, E0, S0 lze vývoj koncentrací určit numericky. Pro numerické výpočty a grafická zobrazení budeme v celé této kapitole uvažovat následující hodnoty parametrů: fci = 0,001 mol"1 dm3 s"1; k, = 0,0005 s"1; k2 = 0,2 s"1; (20) S0 = 1000 mol dm"3; E0 = 100 mol dm"3; C0 = 0 mol dm"3. Numericky získaná řešení počáteční úlohy (15)-(19) s těmito parametry jsou vykreslena v Obr. 1 a 2. Vidíme, že koncentrace S(t) substrátu s rostoucím časem klesá, a to přibližně exponenciálně, a koncentrace P(t) produktu roste. Koncentrace E(t) volného enzymu nejdříve klesá, poté se zvolna vrací k počáteční koncentraci E0. Koncentrace C(t) komplexu naopak nejdříve roste a poté zvolna klesá k nule. Pro ř —> oo dojde ke stabilizaci systému, všechen substrát bude postupně přeměněn na produkt a enzym po ukončení této reakce zůstane opět výhradně ve volné formě. 8 O 20 40 60 80 100 120 t Obrázek 1: Numericky spočítané průběhy (čas v sekundách) koncentrací (mol dm-3) substrátu S (ŕ) (modře), produktu P(t) (zeleně), volného enzymu E(t) (černě) a komplexu enzym-substrát C(t) (červeně) podle kinetického schématu (14) s parametry (20). Koncentrace C(t) a E(t) mají odlišné šká-lování podle měřítka na pravé straně grafu. Pro analýzu úlohy (15)-(19) bychom kromě numerických výsledků potřebovali i přesná řešení v analytickém tvaru. Jejich výpočet je však i pro uvažovanou reakci velmi složitý, běžnými postupy systém vyřešit neumíme, systém diferenciálních rovnic totiž není lineární. Ukážeme proto tzv. quasi-steady-state aproximaci, která poskytuje jednoduchý a užitečný, ale přitom relativně přesný, pohled na interakci enzym-substrát-produkt. Kromě odvození přibližného řešení systému diferenciálních rovnic ukážeme i obecný postup, jak složitý systém zjednodušit na aproximovaný systém a jak zároveň zhodnotit přesnost aproximace a ověřit podmínky pro její platnost. Než se budeme aproximaci věnovat, všimněme si ještě některých závislostí v úloze (15)-(19). Koncentrace produktu P (ŕ) vystupuje pouze v rovnici (18) a k jejímu výpočtu potřebujeme znát průběh koncentrace C(t) komplexu. 9 Obrázek 2: Detail z Obr. 1 zobrazující koncentrace v počáteční fázi reakce. Nejdříve je tedy nutné vyřešit systém první tří diferenciálních rovnic pro koncentrace C(t), E(t), S(t). koncentraci produktu P(t) pak nalezneme integrováním koncentrace komplexu, r t P(t) C(u) du. (21) Sečtením rovnic (15) a (16) dostaneme dE(t) dC(t) 0, dř dř tzn. součet koncentrací C(t) a E(t) je konstantní, podle (19) rovný E(t) + C(t) =E0 + C0. (22) (23) Enzym se tedy v reakci vlastně vyskytuje ve dvou formách: jako volný E a jako vázaný v komplexu C. To je patrné i na časových průbězích koncentrací E(t) a C (ŕ) v Obr. 1 a 2. Z (23) vyjádříme koncentraci E(t) a dosadíme do rovnic (15) a (17). Tím nám v systému zůstanou pouze dvě diferenciální rovnice dC(t) dř fc1(£0 + C0)S(t)-[fc-i + fc2-fciS(t)]C(t)J (24) 10 -^ = k_1C(t)-k1[E0 + C0-C(t)} S (t) (25) pro koncentrace C (t) a. S (t) s počátečními podmínkami C(0) = C0, S(0) = S0. Součet koncentrací substrátu a produktu S (t) +P(t) však obecně konstantní není. Všechen substrát se sice přemění na produkt, ale součet jejich koncentrací je roven S0 pouze na začátku a na konci reakce, jak je patrné z Obr. 3. Obrázek 3: Součet koncentrací (mol dm 3) S(ř) + P(t) substrátu a produktu v čase (sekundy) podle kinetického schématu (14) s parametry (20). 3 Quasi-steady-state aproximace Z původních čtyř rovnic zůstala soustava dvou diferenciálních rovnic (24), (25) s podmínkami (19). Ani tento systém však není řešitelný ve smyslu nalezení analytického tvaru formulí pro obě koncentrace. V takové situaci, jenž je v matematickém modelování celkem běžná, se snažíme pokročit s řešením tím, že některé závislosti v modelu potlačíme a úlohu tak zjednodušíme. Po nalezení řešení se vrátíme do místa, kde jsme model zjednodušili a prodiskutujeme, jak výrazný dopad má zvolené omezení na reálnou situaci a jak by bylo řešení ovlivněné, kdybychom pracovali s původním nezjednodušeným modelem. 11 V Obr. 2 si všimněme toho, že za dobu 5 s od začátku reakce klesne koncentrace substrátu (modrá křivka) asi o 15 %, zatímco koncentrace enzymu (černá křivka) se za stejnou dobu sníží přibližně o 80 %. Vzhledem k rychlosti změn koncentrace enzymu E(t), a tím i koncentrace C (ŕ) komplexu, je tedy koncentrace substrátu S(t) skoro neměnná. Představme si tedy, že reakce je na začátku svého průběhu a předpokládejme poněkud radikální zjednodušení, že koncentrace S(t) = S0 substrátu se v čase nemění, např. proto, že jeho množství je velké, nebo je substrát neustále doplňován. Rovnice (25) je nyní irelevantní a (24) pro jedinou neznámou koncentraci C (ŕ) komplexu přejde do tvaru = fcj(£0 + C0)S0 - (fcj S0 + fc_! + fc2) C(t) (26) s počáteční podmínkou C(0) = C0. Řešením této diferenciální rovnice je funkce C(t) = Cm + (C0-Cm)e-Kt, (27) kde (En + Cn ) Sn Cm = ^-(28) + So je asymptotická hodnota koncentrace komplexu, které je pro ř —> oo dosaženo nezávisle na počáteční hodnotě C0, ^m = ^-- = Ká + -í (29) je tzv. Michaelisova konstanta a K = k1SQ + k_1 + k2. (30) Zdůrazněme fakt, že asymptotická koncentrace komplexu Cm závisí na rychlostních konstantách pouze prostřednictvím Michaelisovy konstanty Km, tedy prostřednictvím vhodného poměru rychlostních konstant. Pro uvažované hodnoty parametrů (20) dostáváme K = 1,2 s_1, hodnotu Michaelisovy konstanty Km = 200,5 mol dm-3 a asymptotickou koncentraci Cm = 83,3 mol dm-3. Vývoj hodnoty koncentrace C(t) za podmínky konstantní koncentrace substrátu S(t) = S0 během prvních 5 s reakce je vykreslena v Obr. 4 (čárkovaná křivka). Pro porovnání je přikresleno i přesné řešení (plnou čarou) z Obr. 1. Vidíme, že 12 koncentraci komplexu C (t) lze velmi dobře aproximovat hodnotou (27). Přitom si můžeme dovolit považovat koncentraci substrátu za konstantní, rovnu S(t) = S0, na intervalu až do doby 5 s po začátku reakce. 0 1 2 3 4 5 t Obrázek 4: Koncentrace komplexu C(t) (mol dm-3) v závislosti na čase ř (sekundy) pro počáteční fázi reakce. Čárkovaná křivka znázorňuje přibližnou hodnotu dle (27) za podmínky konstantní koncentrace substrátu, plnou čarou je vyznačený přesný, numericky spočítaný, průběh koncentrace. Analogicky lze tuto aproximaci aplikovat i v průběhu reakce, jak je ukázáno v Obr. 5, kde je předpoklad konstantní koncentrace substrátu použit po 60 s od počátku reakce. Je však patrné, že aproximace je v tomto případě méně přesná, a lze ji rozumně použít asi po dobu 3 s, pro delší časové intervaly se již hodnota koncentrace Cm značně odlišuje od přesného řešení. S rostoucím časem se koncentrace komplexu mění pomaleji (derivace klesá k nule) a (27) se blíží limitě Cm dle (28). Není tedy překvapením, že Cm je řešením rovnice fc2 (£0 + C0)S0 - (fcj S0 + fc_j + fc2) Cm = 0 (31) pro rovnovážný stav, kterou obdržíme tak, že v (26) uvažujeme dC(t) /dt = 0. 13 u 46 44 42 - 40 38 - 60 61 62 63 64 65 Obrázek 5: Koncentrace komplexu C (t) (mol dm-3) v závislosti na čase ř (sekundy) po 60 s od počátku reakce. Čárkovaná křivka znázorňuje přibližnou hodnotu dle (27) za podmínky konstantní koncentrace substrátu, plnou čarou je vyznačený přesný, numericky spočítaný, průběh koncentrace. Vraťme se nyní zpět ke vztahu (25) pro koncentraci substrátu a uvažujme, co se na výše uvedeném přístupu změní, když S (ŕ) nebude konstantní, ale bude se v čase měnit. Tvrdíme, že v tom případě bude dobrou aproximací koncentrace komplexu C(t) v libovolném čase ř asymptotická hodnota Cm podle (28), v níž namísto S0 dosadíme okamžitou hodnotu koncentrace S (ř). Představíme-li si (28) jako funkci koncentrace substrátu S0, můžeme tedy psát . . (En -\- Cn ) Sn Cm(S0) = V ° "; °. (32) Pokud se koncentrace S(t) substrátu bude měnit pomalu vzhledem k rychlosti změn koncentrace komplexu C (ŕ), bude koncentrace C(t) dobře aproximována tzv. quasi-steady-state rovnicí ÍE0 + C0)S(t) Cit, cm[s(t)} Km+S(t) (33) Tento důležitý výsledek se nazývá quasi-steady-state (pseudo-steady-state) řešením, neboť jej obdržíme také tím, že levou stranu (24) položíme rovnou 0 14 a rovnici vyřešíme pro C(t). Předpokládáme tedy, že koncentrace komplexu je ustálená, dC(ř)/dř = 0, přestože C(t) ve skutečnosti konstantní není. Stačí jen zajistit, aby se koncentrace S(t) měnila tak pomalu, že koncentrace C (ŕ) se chová, jako by S(t) bylo po určitou (krátkou) dobu konstantní. Jaká koncentrace substrátu odpovídá quasi-steady-state aproximaci koncentrace komplexu? Dosazením (33) do (25) obdržíme diferenciální rovnici dS(ř) [k_1 + k1S(t)][E0 + C0]S(t) ,rWrx foA, —7— =-v --fc^Eo + CojSCt). (34) dř Km + S(t) Vytkneme výraz z pravé strany (33) a po úpravě dostaneme dS(t) (E0 + C0)S(t) (E0 + C0) —7— = -h v , g/.x =-fc2-j^- (35) sít: Porovnáním zjistíme, že dS(ř)/dř = —k2C(t), tzn. quasi-steady-state aproximace C (ŕ) udává, po vynásobení rychlostní konstantou k2> rychlost změny koncentrace substrátu. Diferenciální rovnici pro odpovídající časový vývoj koncentrace produktu obdržíme dosazením (33) do pravé strany (18), dP(t) (E0 + C0)S(t) —;— = ^2-—• (36) dř Km + S(t) Řešení rovnice (35) lze nalézt v implicitním tvaru S(t)+KmlnS(t)=S0 + KmlnS0-k2(E0 + C0)t. (37) To nám umožňuje graficky znázornit (např. pomocí parametricky zadaných křivek) průběhy koncentrací substrátu a produktu na čase při quasi-steady-state aproximaci. Výsledné průběhy jsou vykresleny (čárkované křivky) v Obr. 6. Na Obr. 7 je potom vykresleno (čárkovanou křivkou) quasi-steady-state řešení koncentrace komplexu (33). Vidíme, že při splnění výše uvedených podmínek je aproximace pro globální popis průběhu reakce dostačující. Porovnáním derivací (35) a (36) zjistíme, že důsledkem quasi-steady-state aproximace je konstantní součet S (ŕ) + P (t) = S0 koncentrací substrátu a produktu, což se na obrázku projeví vertikální symetrií čárkovaných křivek. Připomeňme, že pro přesná řešení systému diferenciálních rovnic však toto neplatí, viz Obr. 3. 15 O 20 40 60 80 100 120 t Obrázek 6: Koncentrace substrátu S(t) (modře) a produktu P (ŕ) (zeleně) (mol dm-3) v závislosti na čase ř (sekundy). Čárkované křivky znázorňují quasi-steady-state řešení dle (35) a (36), plnou čarou jsou nakresleny přesné, numericky spočítané, průběhy koncentrací podle kinetického schématu (14) s parametry (20). 16 O 20 40 60 80 100 120 t Obrázek 7: Koncentrace komplexu C(t) (mol dm-3) v závislosti na čase ř (sekundy). Čárkovaná křivka znázorňuje quasi-steady-state řešení dle (33), plnou čarou je nakresleno přesný, numericky spočítaný, průběh koncentrace podle kinetického schématu (14) s parametry (20). 17 Úkoly 1. Ukažte, že řešením diferenciální rovnice (26) pro hledanou koncentraci C(t) s počáteční podmínkou C(0) = C0 je funkce (33). 2. Vysvětlete, proč v reakci (14) není součet koncentrací S(t) +P(t) konstantní v čase, jak je zřejmé z Obr. 3. 18 Enzymová kinetika, 2. část Předpoklady Pro studium této výukové jednotky je potřebná znalost • základních principů enzymové kinetiky • diferenciálního a integrálního počtu funkcí jedné proměnné. Výstupy z učení Studenti se naučí • pracovat s pojmem časové škály, • kvantifikovat rychlost chemické reakce, • aplikovat metodologii enzymové kinetiky na podobné biologické modely. 1 Časová škála V předešlé kapitole jsme jako argument pro aproximaci používali předpoklad, že jedna funkce se mění pomaleji než funkce jiná. Tento vágně formulovaný pojem si nyní definujeme přesně. Matematicky je rychlost změny funkce kvantifikována derivací dané funkce. Porovnávání derivací funkcí je však v praxi problematické z důvodu, že hodnoty funkcí a tím i jejích derivací mohou být vyjádřeny v různých jednotkách. Proto neporovnáváme přímo derivace funkcí, ale délky časových intervalů, během nichž dochází k největším změnám v hodnotách funkcí. Dostáváme se tak k pojmu časové škály, time scale. Poznamenejme, že tento pojem je užíván zejména v matematické analýze, kde umožňuje spojit a zobecnit teorii diferenčních a diferenciálních rovnic. Pro náš účel zavedeme pojem time scale, Tf, funkce / (ř) jako nejkratší dobu potřebnou k tomu, aby se funkční hodnota co nejvíce změnila, a to maximální rychlostí. Jde tedy o podíl maximálního možného rozdílu funkčních hodnot ku maximální absolutní hodnotě derivace, \f'(t) \ = |d//dř|, tzn. sup/(ř)-inf/(ř) Tf sup|/'(t)l ' UJ 19 přičemž uvedená infíma a suprema se počítají přes interval časů t, který je předmětem modelování. Obrázek 1: Grafické určení time scale z f (délka modré úsečky) funkce / (ŕ) (červená křivka) tvaru klesající exponenciály. Z matematické analýzy je známo, že derivace f'(t0) funkce v bodě ř0 má geometrický význam směrnice tečny ke grafu funkce / (t) v tomto bodě. Výraz sup \f'\ tedy udává směrnici tečny ke grafu funkce / (t) v bodě, kde je sklon (tzn. absolutní hodnota derivace) funkce největší. Směrnice tečny je dále rovna tangens úhlu a, který tečna svírá s kladnou vodorovnou poloosou. Rovnici (1) tedy můžeme přepsat do tvaru \ť\ * sup/-inf/ sup 1/ | = tg a =-. (2) Tf Pro geometrickou představu dále sledujme schéma na Obr. 1. Červeně vyznačený graf funkce / (ř) tvaru klesající exponenciály má zřejmě největší sklon na svém levém konci. Tam zároveň nabývá svého maxima, tzn. i suprema. Infi-mum klesající exponenciální funkce je rovno nule, tedy výraz v čitateli zlomku v (2) nalézáme jako vyznačenou svislou vzdálenost. V bodě největšího sklonu sestrojíme tečnu ke grafu funkce f(t). Ta svírá s kladnou poloosou ř úhel a, 20 jak je v obrázku vyznačeno. Z definice funkce tangens jakožto poměru délek protilehlé a přilehlé odvěsny v pravoúhlém trojúhelníku a ze zlomku v (2) identifikujeme time scale z f jakožto délku vodorovné odvěsny (vyznačeno modře). Spočítáme time scale tc pro koncentraci komplexu při quasi-steady-state aproximaci. Z vyjádření (EkI.27) vidíme, že C (ŕ) roste z počáteční koncentrace C0 k asymptotické koncentraci Cm, a její derivace C'(t) = K(Cm — C0) e~Kt zřejmě nabývá svého maxima pro ř = 0. Dosazením do (1) spočítáme, že time scale koncentrace komplexu je Cm — Q) 1 t c =--— = —. (3) (Cm — C0)K K Pokračujeme výpočtem time scale ts koncentrace substrátu při quasi-steady-state aproximaci. Ze záporného znaménka derivace (EkI.35), z počátečních podmínek a asymptotického chování podle Obr. 6 je zřejmé, že koncentrace substrátu klesá z počáteční hodnoty S0 na asymptoticky nulovou hodnotu, a že maximum derivace (EkI.35) nastává v čase ř = 0, tedy pro S(t) = S0. Dosazením do (1) tak obdržíme time scale koncentrace substrátu T — SQ — 0 _ Km + S0 S k2(Eo + C0)S0 k2(E0 + C0)' Km + S0 Druhou možností je použít (EkI.25), tedy neaproximovanou diferenciální rovnici pro koncentraci substrátu. Rychlost poklesu koncentrace S(t) je i tomto případě největší pro ř = 0 a dostáváme tak jinou verzi time scale Sq-0 T* = 7 " = , _ u , ^ . (5) V uvažovaném příkladu obdržíme po dosazení parametrů (Ek1.20) následující hodnoty time-scale: tc = 0,833 s, ts = 60,025 s, t* = 10 s. (6) Nyní se vrátíme k podmínce quasi-steady-state aproximace. Požadovali jsme, aby se koncentrace substrátu S(t) měnila mnohem pomaleji, než se mění koncentrace komplexu C(t). Jinak řečeno, čas potřebný pro významnou změnu 21 S (t) má být mnohem delší, než čas potřebný pro významnou změnu C (t). Pomocí time scale obou koncentrací tedy tento požadavek matematicky vyjádříme podmínkou ^<1, resp. ^<1. (7) V uvažované enzymové reakci dostáváme podíly time scale tc/ts = 0,014 a Tc/Ts = 0,083, tzn. koncentrace komplexu C se mění zhruba 10x-100x rychleji než koncentrace substrátu S. Lze ukázat, že vždy platí tJj < ts, druhá podmínka v (7) je tedy silnější než podmínka první. Dosazením odvozených time scale (3), (4), (5) do (7) dostáváme podmínky pro platnost quasi-steady-state aproximace ve tvaru k2(E0 + C0) E0 KáC0 ^ < 1, resp.------ < 1. (8) ^i (^m + So)2 Km + S0 S0(Km + S0 Tyto podmínky jsou v praxi splněny např. pokud je počáteční koncentrace sub stratu řádově větší než koncentrace enzymu, E0 + C0 = -nETŇET (14) Graf této závislosti je vyobrazen v Obr. 6. Pozorujeme, že pro nízké koncentrace L ligandu je nárůst koncentrace (počtu) Cm ligand-receptorových komplexů strmý, s rostoucí koncentrací L rychlost růstu Cm klesá a hodnota Cm se asymptoticky blíží k N, koncentraci (počtu) všech receptoru. 100 F U 4000 6000 8000 10000 Obrázek 6: Závislost (14) rovnovážné hodnoty koncentrace C ligand-receptorových komplexů (mm-2) na koncentraci L ligandu (mol dm-3) při quasi-steady-state aproximaci modelu (13). Koncentrace všech receptoru je N = 100 mm-2, rychlostní konstanty mají hodnoty fcj = 0,001 dm3 s-1 mol-1, fc_j = 0,0005 s-1 a k2 = 0,2 s-1. V praxi se však tato závislost častěji vykresluje s logaritmickou škálou na vodorovné ose, jak je ukázáno v Obr. 7. Vnímání vstupního signálu (zde L) na logaritmické škále je přitom typickým znakem smyslových signálů. Aditivní změnu (posun) na vstupu (fyzikálně měřitelnou) obvykle pomocí smyslového vnímání pociťujeme jako změnu multiplikativní. Připomeňme v této souvislosti 27 -5 0 5 10 15 ln L Obrázek 7: Graf závislosti (14) z Obr. 6 při logaritmickém měřítku na vodorovné ose má tvar sigmoidální křivky (dose-response curve). např. vnímání zvuku, kdy hlasitost zvuku místo jeho intenzity a její SI jednotky W m-2 vyjadřujeme právě jako jako logaritmus (desítkový) poměru intenzity k jisté referenční hodnotě s pomocí bezrozměrné „jednotky" dB. Grafem je při tomto zobrazení sigmoidální křivka (S-křivka). Je patrné, že při velmi nízkých koncentracích L ligandu je vnímaná odpověď, tedy koncentrace Cm ligand-receptorových komplexů, skoro neměnná a nulová. V určitém intervalu koncentrací L se odpověď mění, senzorický systém dokáže změny v hodnotách L rozlišovat. Rozmezí takových koncentrací ligandu se označuje jako coding range. Z hlediska hodnocení schopnosti smyslového vnímání je rozhodující poloha coding range a rychlost změny Cm, tedy hodnota derivace, pro L uvnitř tohoto intervalu. Pro velmi vysoké koncentrace ligandu pak dochází k tzv. nasycení (saturaci), Cm konverguje k N a senzorický systém již na další změny v hodnotách L nereaguje. S podobným modelem i výsledkem se lze setkat také u biofyzikálních či biochemických modelů, kdy je na organismu sledován efekt podání nějaké látky, např. diety či léku, nebo vystavení určitým změnám v prostředí. Sigmoidální křivka znázorňující tuto závislost je proto nazývána také jako křivka dávka- 28 odpověď (dose-response curve). Úkoly 1. Dokažte, že závislost (12) z Obr. 3 má v Lineweaverově-Burkově grafu na Obr. 4 tvar lineární funkce. Jakým hodnotám jsou rovny souřadnice průsečíků přímky s osami? 2. Zapište závislost (12) ve tvaru S-křivky v Obr. 5, tzn. s logaritmickým měřítkem (uvažujte přirozený logaritmus) pro počáteční koncentraci substrátu S0. Spočítejte asymptoty této funkce, souřadnice inflexního bodu a maximální hodnotu derivace této funkce. 29 Difúze Předpoklady Pro studium této výukové jednotky je potřebná znalost • základních fyzikálních principů • náhodné procházky a normálního rozdělení pravděpodobnosti Výstupy z učení Studenti se naučí • porozumět pojmu difúze, a to z pohledu fyzikálně-biologického i matematického • odvodit rovnici difúze Molekuly či buňky ponořené ve vodním prostředí jsou v neustálém pohybu. Veškeré látky mají tendenci přecházet z prostředí se svou vyšší koncentrací do prostředí s nižší koncentrací. Přirozenou vlastností látek je, že pokud se její částice mohou pohybovat (molekuly v nehybném roztoku se pohybují na základě náhodného pohybu), tak se rozptylují do celého prostoru, kterého mohou dosáhnout, a postupně ve všech jeho částech vyrovnají svou koncentraci. Ať živé či neživé, vše podléhá teplotním fluktuacím. Látky difundují určitým směrem tak dlouho, dokud se jejich koncentrace nevyrovná. Proces je stejný např. u solí rozpuštěných ve vodě nebo třeba u kyslíku obsaženém v kompostu. Při tzv. prosté difúze látky samovolně přecházejí samovolně (náhodným pohybem) z prostředí, kde je jejich koncentrace vyšší směrem tam, kde byla dosud jejich koncentrace nižší. Tato kapitola si klade za cíl osvětlit dynamiku difúze živých i neživých systémů a metody pro její studium. Výklad a příklady uvedené v této kapitole čerpají především z monografií Berg (1993), Keener (2009) a Murray (2002). 1 Mikroskopický pohled a náhodná procházka Difúzí rozumíme proces náhodného přesunování molekul nebo malých částic vyplývající z pohybu na základě tepelné energie, přičemž trajektorie pohybu 30 jednotlivých částic lze považovat za náhodné. Jako typické příklady difúze uveďme pohyb molekul nebo jednoduchých organismů, např. bakterií. Částice s absolutní teplotou T (jednotkou je kelvin, K), má v průměru kinetickou (pohybovou) energii (jednotkou SI je joule, J) \kBT (1) podél každé z os v prostoru. Přitom fcB = 1,38 • 10~23 J K_1 je tzv. Boltzman-nova konstanta. Albert Einstein v jednom se svých nejslavnějších článků z roku 1905 dokázal, že tento vztah platí bez ohledu na velikost částice, dokonce i pro ty, které jsou pozorovatelné mikroskopem a vykazují Brownův pohyb. Částice o hmotnosti m a rychlosti v podél x-ové osy má okamžitou kinetickou energii -mv2. (2) 2 Tato sice kolísá, ale v průměru (značíme čarou nad veličinou) musí platit -m^2 = -kBT, (3) 2 2 odkud můžeme vyjádřit střední kvadratickou rychlost částice, * = m. (4) m Rovnici (4) lze použít pro odhad rychlosti malé částice, například molekuly lysozymu. Lysozym je enzym, který se vyskytuje např. ve slinách, slzách, vaječném bílku, nosním hlenu, krevní plazmě a mateřském mléce, a má silné antibakteriální účinky. Hmotnost 1 molu lysozymu, je 14 kg. Jeden mol má přitom 6,022 • 1023 molekul, toto množství je dáno tzv. Avogadrovou konstantou. Hmotnost jedné molekuly lysozymu je tedy 1,4 -104 9n m = —'-=r g = 2,32 • 10"20 g. 6,022 • 1023 6 6 Při teplotě 20 °C, tedy T = 293,15 K, dostaneme dle (4) odhad rychlosti -\ ^-^l,33-103 cms_1. (5) \ m Pokud by molekule lysozymu nebyly kladeny žádné překážky v pohybu, urazila by průměrně více než 13 metrů za 1 sekundu. Protože však protein není 31 ve vakuu, ale ve vodním prostředí, nedostane se molekula lysozymu daleko, neboť záhy narazí do molekuly vody a odrazem změní směr pohybu. Je tak nucena se pohybovat způsobem, který matematicky popisujeme jako náhodnou procházku. Pokud by větší množství takových částic bylo na počátku koncentrováno v malé oblasti, postupem času by se částice pohybovaly v různých směrech a šířily by se do prostoru. Tento jev nazýváme difúzí. Abychom dokázali matematicky popsat proces difúze, omezíme požadavky modelu difúze jen na nezbytné minimum. Budeme tak předpokládat, že částice se mohou pohybovat pouze podél jedné prostorové osy, x. Na počátku pozorování, v čase ř = 0, jsou všechny částice v bodě x = 0 a následně se jejich pohyb řídí zobecněnou náhodnou procházkou s následujícími pravidly: (i) Každá částice se pohybuje rychlostí v, a jednou za čas Ař se na ose x posune o vzdálenost Ax = ±v Ař, buď vpravo, Ax = v Ař, anebo vlevo, Ax = —v Ař. Pro jednoduchost přitom Ař a Ax považujeme za konstanty. Reálně by tyto veličiny mohly záviset na velikosti částice, druhu tekutiny (plynu či kapaliny) tvořící zkoumané prostředí a na teplotě. (ii) Pravděpodobnosti posunu vpravo i vlevo jsou v každém časovém kroku rovny 1/2. Částice se v tekutině pohybují nezávisle na historii svých trajektorií, po nárazu do molekuly tekutiny mění směr zcela náhodně. Po sobě jdoucí kroky pohybu částic jsou tedy stochasticky nezávislé. (iii) Jednotlivé částice se pohybují nezávisle na sobě, nereagují spolu. To je v praxi pravda, pokud suspenze částic je dostatečně zředěná. Tato pravidla popisují symetrickou zobecněnou náhodnou procházku, tedy s parametrem p = 1/2. Od jednoduché náhodné procházky se liší v tom, že časové a prostorové kroky kroky nemají jednotkovou velikost. Časový krok má délku Ař a částice se na reálné ose pohybují jen po násobcích délky Ax, jak je schematicky naznačeno na Obr. 1. Polohu částice v čase ř označme jako náhodnou veličinu X(t) diskrétního typu. Složením pravidel (i)-(iii) dostáváme rovnici X(t + At)=X(t)+Z(t), (6) 32 x-2Ax x-Ax x x+Ax x+2Ax Obrázek 1: Schéma symetrické zobecněné náhodné procházky. Za časový interval délky Ař se částice z bodu x přesune buď do bodu x — Ax s pravděpodobností 1/2, anebo do bodu x + Ax, taktéž s pravděpodobností 1/2. kde f i 2, z = Ax P{Z(t)=z} = h, z = -Ax (7) 0, z £ ±Ax a náhodné veličiny X(t), Z(t), X(s), Z(s), t ^ s, jsou vzájemně stochasticky nezávislé. Lehce spočítáme střední hodnotu a rozptyl náhodné veličiny Z, E[Z(ř)] = \ Ax-\ Ax = 0, (8) Var[Z(ř)] = E[Z2(ř)] = ±(Ax)2 + ±(-Ax)2 = (Ax)2. (9) Spojením (6), (8) a (9) nyní spočítáme střední hodnotu a rozptyl polohy částice po jednom časovém kroku, E[X(t + At)]=E[X(t)]+E[Z(t)}=E[X(t)}, (10) Var[X(ř + Ař)] = Var[X(t)] + Var[Z(ř)] = Var[X(ř)] + (Ax)2. (11) Vidíme, že střední poloha částice se v čase nemění, zatímco rozptyl s každým krokem narůstá o konstantní hodnotu (Ax)2. Pro jednoduchost dále předpokládejme, že částici začínáme pozorovat v čase ř = 0, a spočítejme rozptyl její polohy po n krocích, tedy v čase ř = n Ař. Pomocí rekurentního vztahu (11) lehce spočítáme (Ax)2 Var[X(t = nAt)] = n(Ax)2 = t v y , (12) tedy rozptyl polohy částice je přímo úměrný první mocnině času ř. Pokud budou časové kroky Ař dostatečně krátké, zmenší se i délka kroků Ax, a zlomek 33 na pravé straně (12) bude konvergovat k jistému kladnému číslu. Na základě tohoto pozorování se definuje tzv. difúzni koeficient (Ax)2 D= lim --—, (13) At^o+ 2 Ar pomocí něhož potom můžeme psát Var[X(t)] = 2Dt. Difuzní koeficient D charakterizuje pohybové možnosti částic daného typu v určité tekutině při určité teplotě. Obecně proto závisí na velikosti částic, struktuře tekutiny a teplotě. Rozměr veličiny D je rovný jednotce plochy dělené jednotkou času, v praxi se často uvádí v cm2 s_1. Pro molekuly malých rozměrů ve vodě při pokojové teplotě je difuzní koeficient řádově rovný D & 10~5 cm2 s_1. Při takové hodnotě difuzního koeficientu se částice posune o vzdálenost Ax = 1 /im, tj. zhruba délky bakterie, za čas (Ax)2 4 Ař =--— w 5 • 10"4 s = 0,5 ms. 2D Urazit 104krát větší vzdálenost Ax = 1 cm by jí však trvalo (Ax)2 4 Ař = --— s* 5 • 104 s, 2D tedy téměř 14 hodin, přesně 108krát delší čas. Změna polohy není úměrná času, ale jeho druhé odmocnině. Při difúzi ve vícerozměrném prostoru se pravidla (i)-(iii) použijí pro každou ze souřadnic, x, y, příp. z. Navíc přitom předpokládáme, že průměty pohybu částic do těchto os jsou stochasticky nezávislé. Ve dvourozměrném prostoru tak pro střední kvadratickou vzdálenost od počátku, kdy R2 = X2 + Y2, platí E[R2(t)} =4Dt, a ve třírozměrném prostoru, kde R2 = X2 + Y2 + Z2, pak E[R2(t)} =6Dt. Nyní nás bude zajímat rozdělení pravděpodobnosti polohy částice po n krocích, tedy v čase ř = nAř. (14) 34 Použijeme opět matematický popis pomocí (6) a (7). Ze znalosti délky časového kroku A ŕ dokážeme určit, kolik kroků pohybu částice proběhlo, označme tento počet V každém z těchto kroků se náhodná veličina Z dle (7) realizovala buď jako posun vpravo o Ax, anebo jako posun vlevo o Ax. Označme jako náhodnou veličinu N+ počet těch z celkem n realizací, které byly posuny vpravo. Podle pravidel (i), (ii) je zřejmé, že tato náhodná veličina má binomické rozdělení pravděpodobnosti s parametry n a 1/2, tedy má pravděpodobnostní funkci P{JV+ = m} = r^Vljm, m = 0,l,...,n. (15) Spočítáme střední hodnotu a rozptyl náhodné veličiny N + E(iV+) = -, Var(iV+) = -. (16) 2 4 Při modelování procesu difúze náhodnou procházkou musíme brát do úvahy, že časový krok Ař musí být dostatečně krátký. Pokud v (14) položíme Ař —> 0, počet kroků n poroste nade všechny meze. K aproximaci binomického rozdělení za této situace využijeme centrální limitní větu. Ta říká, že rozdělení pravděpodobnosti standardizované náhodné veličiny N+ pro n —> oo konverguje ke standardizovanému normálnímu rozložení, N+-E(N+) y/Var(N N 0,1 Konkrétně, po dosazení (14) a (16) pak z vlastností lineární transformace plyne, že rozdělení pravděpodobnosti náhodné veličiny N+ lze pro Ař —> 0 aproximovat normálním rozdělením N+ ~ NÍ——,——V (17) Pokud při pohybu částice v n krocích došlo k N+ posunům vpravo o Ax, plyne z předpokladů difúze, že muselo dojít k n — N+ posunům vlevo. Poloha částice po n krocích je tak dána rozdílem Ař+Ax + (n-JV+)(-Ax) = 2AxiV+-nAx, neboli, mezi X (ŕ) aN+ existuje jednoznačný lineární vztah X(t) = 2AxN+-— Ax. (18) w Ař 35 Ze znalosti rozložení (17) již nyní lehce dopočítáme asymptotické rozdělení pravděpodobnosti polohy částice X(t), X(t) ~ Nl °» Aŕ ŕ J = N (0,2Dŕ). (19) Za podmínek modelu difúze pomocí náhodné procházky tak poloha částice, X(t),v čase ř má normální rozdělení s nulovou střední hodnotou a rozptylem 2Dt, lineárně závislém na čase a na hodnotě difuzního koeficientu D. Hustota pravděpodobnosti tohoto rozdělení je rovna f(x,t) = -z^=exp\——\. (20) Podle pravidel (i), (ii) se za časový interval Ař jednoho kroku částice do místa x mohla dostat s pravděpodobností 1/2 z místa x — Ax, anebo s pravděpodobností 1/2 z místa x + Ax. Použitím formule úplné pravděpodobnosti tak můžeme psát P{X(ř + Ař) = x} = |P{X(ř) = x-Ax} + |p{X(ř) = x + Ax}. (21) Zaveďme pro tyto pravděpodobnosti označení p(x,ř) = P{X(t) =x}, xeR,t>0. (22) Rovnici (21) přepíšeme do tvaru p(x,ř + Ax) = ±[p(x-Ax,ř)+p(x + Ax,ř)]. (23) Od obou stran nyní odečteme výraz p(x, ř), zlomek na pravé straně rozšíříme výrazem (Ax)2 a celou rovnici vydělíme Ař. Po přeuspořádání výrazů na obou stranách dostaneme p(x, t + Ař)-p(x, ř) (Ax)2 p(x-Ax, ř) +p(x + Ax, ř)-2p(x, ř) Ař ~~ 2Ař (Ax)2 ' (24) (Ax)2 Provedeme limitní přechod Ař —* 0, Ax —* 0 tak, aby podíl Ař byl konečný. Na levé straně (24) tak vznikne první parciální derivace funkce p podle času ř. První zlomek na pravé straně je roven difuznímu koeficientu D. Limita druhého 36 zlomku na pravé straně je rovna druhé parciální derivaci funkce p vzhledem k prostorové proměnné x, jde o hlavní člen tzv. tříbodové formule. Vztah (24) tak uvedeným limitním přechodem přejde na rovnici kterou nazýváme rovnicí difúze. Není náhodou, že tuto parciální diferenciální rovnici splňuje právě hustota normálního rozdělení pravděpodobnosti (20). 2 Makroskopický pohled a Fickův zákon V předchozí části jsme zkoumali mikroskopický pohyb jedné konkrétní částice. Proces difúze jsme modelovali zobecněnou náhodnou procházkou, rovnici difúze a její obecné řešení jsme nalezli zkoumáním pravděpodobnostních vlastností tohoto modelu. V této části rovnici difúze odvodíme pomocí makroskopické teorie difúze. Prostředkem matematického modelování difúze bude popis koncentrace zkoumaných částic (objektů) a jejich změn v čase a prostoru. Při odvození se opět omezíme na jednorozměrný prostor, reprezentovaný dutým válcem o konstantním průřezu S umístěným podél x-ové osy. Než odvodíme tzv. Fickův zákon a rovnici difúze, popíšeme si použité fyzikální veličiny. Koncentrací C(x, ř) rozumíme počet částic na jednotkový objem v místě x a čase ř. Počet částic v určité oblasti je pak rovný součinu koncentrace a objemu oblasti, za předpokladu, že koncentrace je ve všech místech stejná, což bude splněno např. když oblast bude dostatečně malá. Počet částic v části dutého válce o průřezu S umístěného podél x-ové osy a o dostatečně malé délce Ax je tedy přibližně roven Rychlostí ryzího vzniku rozumíme veličinu Q(x, t), udávající rozdíl počtu nově vzniklých (dodaných) a zaniklých (odebraných) částic na jednotkový objem a jednotkový čas v místě x a čase t. Je-li Q(x, t) > 0, částice v dané oblasti vznikají, je-li Q(x, t) < 0, částice v dané oblasti zanikají. Rozdíl počtu vzniklých a zaniklých částic v části dutého válce o průřezu S a dostatečně malé délce Ax za časový interval At je přibližně roven součinu rychlosti ryzího vzniku, (25) C(x, ř)SAx. 37 objemu oblasti a délky časového intervalu, Q(x,ŕ)SAxAŕ. Veličina čistý tok J(x0, t) popisuje množství částic, které za jednotku času překročí plochu o jednotkovém obsahu kolmou na osu x v daném místě x0 a čase t, přičemž částice procházející touto plochou ve směru kladné x-ové poloosy se započítávají s kladným znaménkem a částice procházející ve směru opačném se započítávají se znaménkem záporným. Je-li J(x0, ř) > 0, větší množství částic se pohybuje přes plochu x = x0 v kladném směru, pokud J(x0, ř) < 0, větší množství částic se pohybuje přes plochu x = x0 ve směru záporném. Rozdíl počtu částic, které projdou plochou o obsahu S kolmou na osu x v bodě x0 za časový interval délky Ař v kladném a záporném směru, je pak rovný součinu J(x0,t)S Ař. Z uvedených vyjádření lze již nyní určit jednotky uvedených veličin. Zůstaneme u obvyklého vyjádření počtu částic pomocí molárního množství, tedy jednotkou mol. Budeme-li délku Ax vyjadřovat v centimetrech a obsah S v cm2, koncentrace C(x, ř) vyjadřující počet částic v jednotkovém objemu bude mít jednotku mol cm-3. Aby platily výše uvedené vztahy, musí pak mít rychlost ryzího vzniku Q (x, ř) jednotku mol cm-3 s-1 a čistý tok J(x, ř) musí být vyjádřen v jednotkách mol cm-2 s_1. Vezměme si dvě nepříliš vzdálená místa trubice, x0 a x0 + Ax, Ax > 0. Předpokládejme, že známe koncentrace částic v těchto místech v čase ř a klademe si otázku, kolik částic se přemístí přes jednotkovou plochu mezi těmito dvěma body za časový okamžik Ař > 0. Jaký je čistý tok J(x, ř) trubicí v prostoru mezi uvažovanými dvěma body? Sledujme schéma na Obr. 2. Pohyb částic modelujeme symetrickou zobecněnou náhodnou procházkou. Během časového intervalu [ř, ř + Ař) se polovina částic z bodu x0 přesune do bodu x0 + Ax, tedy v kladném směru osy x, těch je N+ = \S AxC(x0,ř). Naopak polovina částic původně v místě x0 + Ax se dostane do bodu x0, tedy v záporném směru osy x, jejich počet je N~ = \ S AxC(x0 + Ax,ř). 38 Obrázek 2: Schéma k odvození 1. Fickovy rovnice, (28). Uvažujeme válec o průřezu S a délce Ax umístěný podél x-ové osy, s čely v bodech x0 a x0 + Ax. V místě levého čela je S C(x0, ř) částic na jednotkovou délku válce, analogicky v místě pravého čela je S C(x0, ř) částic na jednotkovou délku válce. Za dobu Ař přejde N+ částic z levého čela doprava a N~ částic z pravého čela směrem doleva. Tok J(x, ř) částic bodem x uvnitř válce obdržíme uvážením chování při limitního zkrácení válce na nulovou délku. Čistý tok válcem mezi body x0 a x0 + Ax je dán rozdílem N+ — N~ vyděleným plochou S průřezu válce a délkou Ař časového intervalu, J(x'ř) = š~a7 -Ar]= 2JÄI [c(x°'ř) s Ax ~ c(x°+ Ax'ř) s Ax] ■ (26) Zlomky na pravé straně zkrátíme obsahem průřezu S a rozšíříme členem Ax. Po přeuspořádání dostaneme (Ax)2 C(x0 + Ax,ř)-C(x0, ř) J(x,t)=---J-^^-'—L-(27) v ; 2Ař Ax Provedeme limitní přechod Ax —* 0+, Ař —* 0+ tak, aby podíl v Ař; zůstal konečný. Tím přejde bod x uvnitř uvažovaného intervalu do bodu x0. V prvním zlomku na pravé straně (27) rozpoznáme difuzní koeficient (13), druhý zlomek přejde v parciální derivaci koncentrace podle x v bodě x0. Po formálním přeznačení x0 na obecný bod x (na volbě polohy uvažované části válce nezáleží) obdržíme tzv. 1. Fickovu rovnici , \ dC(x,t) J(x,t)=-D ; '. (28) ax Ta vyjadřuje tzv. Fickův zákon, který říká, že čistý tok v bodě x a čase ř je úměrný prostorové derivaci koncentrace v tomto bodě. přičemž konstanta úměrnosti je zápornou hodnotou difuzního koeficientu, —D. Všimněme si záporného 39 znaménka, které vyjadřuje že kladný tok směřuje do místa s nižší koncentrací. Zároveň je zřejmé zahrnutí faktoru 1 /2 do definice difuzního koeficientu D dle (13), aby zápis (28) by co nejjednodušší. Z (28) také lehce ověříme, že jednotkou D musí být cm2 s_1. Poznamenejme, že obecně může být difuzní koeficient funkcí koncentrace C(x, t), resp. polohy x a času ř. V takovém případě se často používá k aproximaci závislosti D na C(x, t) Taylorova rozvoje vhodného řádu. Pokud je parciální derivace koncentrace částic v bodě x0 podle x nulová, je v bodě x0 tok nulový, J(x0, ř) = 0. Prostorové rozložení částic v okolí bodu x0 se nemění, systém je v rovnováze. Zdůrazněme, že rozložení částic se však stále může měnit v čase, parciální derivace vhledem k ř není nijak omezena. Pokud je v pevném čase ř0 tok v bodě x0 konstantní, J(x0, t0) = a, má v tomto čase prostorové rozložení koncentrace částic v okolí bodu x0 tvar lineární funkce C(x, t0) = — ax + b. Pro další odvození budeme opět uvažovat část B válce mezi body x0 a x0 + Ax, umístěného podél x-ové osy, jak je ukázáno na Obr. 3. Obě čela (podstavy) válce mají obsah S. Objem zkoumané části válce je tak rovný S Ax. Obecný zákon rovnováhy říká, že změna počtu částic v B o objemu S Ax za dobu Ař je rovna součtu počtu nově vzniklých částic uvnitř B a čistého počtu částic, které přešly přes obě čela (hranice) dovnitř B. Obrázek 3: Schéma k odvození 2. Fickovy rovnice, (33). Uvažujeme válec B o průřezu S a délce Ax umístěný podél x-ové osy, s čely v bodech x0 a x0 + Ax. Toky částic čely válce B jsou rovny J(x0, t), resp. J(x0 + Ax, t), čistá rychlost vzniku částic (vnější příčinou) v bodě x uvnitř B je rovna Q(x, t). x xo x0+Ax 40 Počet částic v B je (obecně při nekonstantní koncentraci) rovný Xq + Ax N(t)=S J C(x,t)dx. Změna počtu částic za časový interval [t, t + At) je rovna rozdílu iV(ŕ +Aŕ)-iV(ŕ). Čistý počet nově vzniklých částic uvnitř B za dobu A ŕ je dán výrazem x0+A SAt J Q(x,ŕ)dŕ. Čistý počet částic, které přešly přes čelo x = x0 dovnitřB je rovný J(x0, t) S At. Analogicky, přes čelo x = x0 + Ax přejde —J(x0 + Ax, t) S At částic dovnitř B. Záporné znaménko je zde proto, že tok částic čelem x = x0 + Ax dovnitř B je v negativním směru osy x. Matematický zápis zákonu rovnováhy je tedy tvaru Xq + Ax Xq + Ax S J C(x,t + At)dx-S J C(x,t)dx = x0 x0 Xq + Ax = AAř J Q(x,t)dx + J(x0,t)S At-J(x0 + Ax,t)S At. (29) xo Rovnici vydělíme výrazem S At, čímž přejde do tvaru Ä7 Xq + Ax Xq + Ax C(x,ř + Ař)dx- C(x,ř)dx *0 *0 Xq + Ax Q(x, ř) dx + J(x0, ř) - J(x0 + Ax, ř). (30) xo 41 Za předpokladu spojitosti použitých funkcí použijeme ve všech třech integrálech Lagrangeovu větu o střední hodnotě. Podle ní existují tři (ne nutně různé) body Xi,x2,x3 e [x0,x0 + Ax] takové, že (30) lze napsat ve formě C(x1,t + At)Ax-C(x2,t)Ax . . Ař Celou rovnici vydělíme Ax, (31) C(x1,ř + Ař)-C(x2,ř) J(x0 + Ax,ř)-J(x0,ř) -Ät-= Q(X3'f)--A~x-' (32) a provedeme limitní přechod Ax —> 0+, Ař —> 0+. Tím body xl5 x2, x3 přejdou v x0, na levé straně dostaneme parciální derivaci koncentrace podle času a na pravé straně parciální derivaci toku podle x. Po přeznačení x0 na obecné x dostáváme tzv. 2. Fickovu rovnici 3C{x,t) dJ(x,t) -—— = Q(x,ř)----. (33) dt ax Dosazením (28) do pravé strany (33) a za podmínky Q(x, ř) = 0 obdržíme dC(x,t) d2C(x,t) \ ) = D ^ \ (34) dt dxz což je tzv. rovnice difúze. Porovnáním s (25) zjistíme, že se matematicky jedná o tu samou diferenciální rovnici. Tvar (25) jsme však odvodili pro pravděpodobnost výskytu jedné částice, zatímco (34) popisuje koncentraci všech částic v daném místě a čase. Obecný tvar řešení této parciální diferenciální rovnice je uveden v úkolu 1, více se mu pak budeme věnovat v následující kapitole. S touto rovnicí se lze setkat nejen u difúze. Pokud veličina C(x, ř) popisuje teplotu tyče v daném místě a čase, nazývá se (34) rovnicí vedení tepla. Kromě modelování biochemických a biofyzikálních procesů se rovnice difúze využívá např. také ve finanční matematice pro popis závislosti ceny akcie a opce. Z pohledu matematické analýzy se jedná o parciální diferenciální rovnici 2. řádu parabolického typu. Jako u jedné z mála parciálních diferenciálních rovnic známe exaktní postupy pro nalezení řešení za různých počátečních, koncových a okrajových podmínek. Tvary řešení některých z těchto situací si ukážeme v další kapitole tohoto výukového materiálu. 42 Ve vyšších dimenzích prostoru se difuzní rovnice odvozuje analogicky, přičemž výše uvedené principy se aplikují složky pohybu částic ve směrech jednotlivých os. Ve dvourozměrném prostoru tak má rovnice difúze pro koncentraci C(x,y, t) v bodě (x,y) a čase ř tvar dC(x,y,t) dt D 82C(x,y,t) 82C(x,y,t + (35) dx2 dy2 Podobně, ve třírozměrném prostoru se rovnice difúze pro popis koncentrace částic C(x,y,z, t) obvykle zapisuje ve tvaru 8C(x,y,z,t) dt DVzC(x,y,z,ř), kde vystupuje tzv. Laplaceův operátor V2 (nabla na druhou), , d2 d2 d2 (36) (37) dx2 dy2 dz2 Na závěr uvedeme ještě speciální, ale v praxi častý, případ difúze v třírozměrném prostoru. Jedná se o sféricky symetrický model, kdy jsou podmínky difúze stanoveny tak, aby koncentrace částic v prostoru byla vždy v konkrétním čase konstantní na celé sféře kolem počátku. Pro popis modelu pak místo tří prostorových souřadnic x,y,z stačí jen jedna proměnná, a to vzdálenost od počátku r = y/x2 + y2+z2 > 0. Rovnice (36) se klasickou transformací z kartézských do sférických souřadnic převede na tvar dC(r,t) dt D 2 dC(r,t) d2C(r,t) r d r d r2 (38) Obecný tvar řešení této parciální diferenciální rovnice je uveden v úkolu 2, vrátíme se k němu i v následující kapitole. Úkoly 1. Ověřte, že hustota normálního rozdělení pravděpodobnosti N(0,4Dt), 1 f x2\ V4nDt exp 4Dt 43 řeší pro ieR, t >0, D > 0 parciální diferenciální rovnici difúze df{x,t) Dd2f(x,t) dt dx2 2. Ověřte, že funkce f(r,t) exp v/(47rDt)3 * V 4Dt, řeší pro r > 0, t >0, D > 0 parciální diferenciální rovnici df(r,t) dt D 2df(r,t) | d2f{r,t) r dr dr2 44 Řešení rovnice difúze, 1. část Předpoklady Pro studium této výukové jednotky je potřebná znalost • základních principů řešení diferenciálních rovnic • tvaru a odvození parciální diferenciální rovnice difúze Výstupy z učení Studenti se naučí • používat Diracovu delta-funkci • porozumět řešení rovnice difúze v ID 1 Úvod Náš výklad řešení rovnice difúze pro jednorozměrný prostor začneme nejjed-nodušším případem, kdy předpokládáme, že chceme popsat časový vývoj koncentrace částic v jednotlivých místech nekonečně dlouhé trubice. Ačkoliv se na první pohled může zdát, že tato situace je ryze teoretická a pro praxi nepoužitelná, ukazuje se, že je i dobrou aproximací situace, kdy trubice je sice konečné délky, ale ta část trubice, která je zájmem modelování, je natolik vzdálena od konců trubice, že vliv tzv. okrajových podmínek můžeme zanedbat. Matematicky přesné řešení rovnice difúze ponecháme až zcela na závěr, ukáže se tak totiž, jak principiálně jiné a na výpočet mnohem náročnější je přesné řešení s uvažovanými okrajovými podmínkami. Připomeňme parciální diferenciální rovnici 2. řádu (DlF.34), resp. (DlF.25). V této difúze, nazývané též rovnicí vedení tepla, vystupuje jedna prostorová proměnná, x, a časová proměnná, ř > 0. Konstanta D > 0 značí difuzní parametr, který kvantifikuje, jak rychle se částice v uvažovaném objektu (trubici) šíří. ř>0 (1) 45 V našem výkladu se nebudeme zabývat řešením uvedené diferenciální rovnice, které je z velké míry založeno na Laplaceově a Fourierově transformaci a jisté dávce praktických zkušeností s využitím těchto integrálních transformací. Zvídavého čtenáře odkazujeme na některou z monografií či skripta o parciálních diferenciálních rovnicích. Naopak se více zaměříme na interpretaci řešení a vliv parametrů pro několik praktických počátečních, resp. okrajových, úloh. Dále si všimneme, že v řešení většiny úloh určitým způsobem figuruje hustota nebo distribuční funkce normálního (Gaussova) rozdělení pravděpodobnosti. Předpokládáme, že s těmito pojmy se čtenář dosud setkával výhradně v kurzech pravděpodobnosti a matematické statistiky, tedy ve spojení s náhodnými veličinami. Zde však žádná náhoda není, tvar rovnice difúze i počáteční, příp. podmínky, jsou jednoznačně zadány. I přesto hraje normální (Gaussovo) rozdělení velkou roli v řešení rovnice difúze. Základ tohoto možná poněkud překvapivého propojení hledejme v závěrečné části předchozí kapitoly. Právě tam jsme si ukázali dva matematicky ekvivalentní pohledy na difúzi: mikroskopický, založený na náhodné procházce, resp. Brownově pohybu, jednotlivých částic, a makroskopický, popisující koncentraci částic v určitém místě a čase. K matematickému popisu některých podmínek budeme potřebovat pojem tzv. Diracovy delta-funkce (či jen delta-funkce) 5(x). Jednou z možností je definovat tuto funkci jako limitu hustoty pravděpodobnosti rovnoměrného rozdělení pravděpodobnosti na intervalu [a — h, a + h], a e IR, h > 0, pro h —> 0+, Alternativní definice Diracovy delta-funkce je ve tvaru hustoty pravděpodobnosti singulárního normálního rozdělení pravděpodobnosti, tedy jako limitní případ normálního rozdělení se střední hodnotou a e IR a s rozptylem cr2 konvergujícím k nule, Z obou definic je patrné, že Diracova delta-fukce není funkcí v klasickém pojetí, (2) (3) 46 neboť pro x = a nabývá nekonečné hodnoty (není reálným číslem), 5(x-a) = { ' ' (4) + 00, x = a 0, x ^ a Navzdory nekonečné hodnotě integrandu v jednom bodě však v důsledku definic založených na hustotách pravděpodobnosti pro Lebesgueův integrál delta-funkce platí 5(x — a) dx = 1. (5) Při řešení rovnice difúze za různých podmínek budeme dále využívat následující vlastnost delta-funkce, nazývanou sampling property. Pro spojitou funkci f (x) platí f (x) 5(x — a) dx = f (a), (6) tedy integrál součinu funkce f(x) s delta-funkcí 5(x — a) je rovný hodnotě f (a) této funkce v bodě a. Pro přehlednost a zjednodušení zápisu budeme v celé této kapitole dále používat značení f(x, t) pro hustotu centrovaného normálního rozdělení pravděpodobnosti s rozptylem 2Dt, „2 f(x,t) V4nDt a F(x, t) pro odpovídají distribuční funkci, exp 4Dt F(x,ť. f(u, t) du. (7) (8) 2 Nekonečná trubice Uvažujeme nekonečně dlouhou trubici, přičemž počáteční koncentrace částic v ní je popsána funkcí g(x), C(x,0)=g(x), X6R. (9) Difuzní rovnice (1) pro x e IR s počáteční podmínkou (9) má řešení C(x,t) V4nDt exp [x — u) 4Dt g(u)du. (10) 47 12 h = 1/20 10 h= 1/8 8 h = 1/4 h = 1/2 6 4 2 0 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 x Obrázek 1: Schéma konvergence hustoty rovnoměrného rozdělení k delta-funkci. Obsah obdélníka je stále rovný jedné, při zmenšující se šířce (2fr) tak jeho výška (l/2fr) konverguje k nekonečnu. Všimněme si, že člen s odmocninou a exponenciální funkce má tvar hustoty normálního rozdělení pravděpodobnosti se střední hodnotou u a rozptylem 2Dt. V souladu se zavedeným označením (7) tedy lze psát Všimněme si rozměrů zde vystupujících veličin. Délkové míry budeme uvádět v centimetrech a čas v sekundách. Koncentrace C(x,t) a g(x) mají rozměr mol cm-3, což odpovídá chemické jednotce molární koncentrace 1000 M. Veličina (7) vystupující v (11) pak musí mít rozměr cm-1. Zdůrazněme, že hustota pravděpodobnosti není bezrozměrnou veličinou, v praxi má vždy rozměr rovný převrácené hodnotě jednotky veličiny v argumentu hustoty, v našem případě délky. Teprve její integrál, distribuční funkce (8), je bezrozměrnou veličinou nabývající hodnot z intervalu [0,1]. (ID 48 Uvažujme dále speciální případ počáteční podmínky — 5(x-a; (12) kdy, v jinak prázdné trubici s plochou průřezu S, je v čase ř = 0 vypuštěno N0 částic v místě a trubice. Množství částic N0 předpokládáme uváděné v molech, plochu průřezu trubice v cm2, tedy delta-funkce musí mít rozměr cm-1. Je to v souladu s definicí delta-funkce jako limitního případu hustoty pravděpodobnosti. Delta-funkce v počáteční podmínce určuje prostorovou polohu zdroje na x-ové ose, má tedy rozměr převrácené hodnoty jednotky délky. Dosadíme (12) do (11) a obdržíme f °° N C(x,ř)= f (x— u,t) — 5 (x— a) du, J-oo S (13) což s využitím (6) dává řešení C(x,t) ■a, t SV4nDt exp (x — a: 4Dt (14) 49 ve tvaru N0/S -násobku hustoty normálního rozdělení se střední hodnotou a a rozptylem 2Dt. Graf řešení je je zobrazen v Obr. 3. Rozložení koncentrace v trubici v několika časových okamžicích je zachyceno v Obr. 4. Koncentrace má v každém čase přesně tvar hustoty normálního rozdělení se střední hodnotou a a rozptylem 2Dt lineárně rostoucím v čase. Časové vývoje koncentrace v několika místech trubice jsou zachyceny v Obr. 5. V místě zdroje, a, v trubici koncentrace exponenciálně klesá, v ostatních místech nejdříve roste a později klesá. Nejvyšší koncentrace je vždy v místě zdroje, a, směrem od tohoto místa symetricky klesá. Obrázek 3: Koncentrace C(x, ř) v mM = 10~6 mol cm-3 dle (14) v závislosti na prostorové souřadnici x e IR v cm a na čase ř > 0 v sekundách. Na počátku je v místě a = 0,5 cm vypuštěno NQ = 10~6 mol částic, jejichž difuzní koeficient je roven D = 10~5 cm2 s_1. Trubice má konstantní průřez o ploše S = 0,1 cm2. Uvažujme další příklad, se dvěma počátečními zdroji, s částicemi v bodě 50 0.0 0.2 0.4 0.6 0.8 1.0 Obrázek 4: Rozložení C (x, t) v mM z Obr. 3 v závislosti na x (v cm) dle (14) v několika časových okamžicích. 0 2000 4000 6000 8000 t Obrázek 5: Vývoj C(x, t) v mM z Obr. 3 v čase (v sekundách) dle (14) v čase t v několika místech (v centimetrech) trubice. 51 a: a s Aí2 částicemi v bodě a2 trubice. Řešíme tedy (1) na x e R s počáteční podmínkou g(x) = y5(x-a1) + y5(x-a2). (15) Parciální diferenciální rovnice difúze je lineární, tzn. pro řešení platí tzv. princip superpozice. Při lineární kombinaci počátečních podmínek získáme řešení stejnou kombinací řešení, které odpovídají jednotlivým počátečním podmínkám. V tomto případě tedy zkombinujeme dvě řešení tvaru (14) a obdržíme řešení úlohy (1), (15) ve tvaru iVi Nn C(x,t) = -±f(x-ai,t) + -Zf(x-a2,t). (16) Jeho průběh je zobrazen na grafu v Obr. 6. Rozložení koncentrace v trubici v několika časových okamžicích je zachyceno v Obr. 7. Řešení (16) je směsí násobků hustot normálních rozdělení se středními hodnotami a: a a2 a s rozptyly lineárně rostoucími v čase. Prostorový průběh koncentrace se postupně vyhlazuje, po určitém čase (v závislosti na vzdálenosti zdrojů) dojde ke změně z bimodálního tvaru na tvar unimodální. Předchozí úlohu se dvěma zdroji lze principem superpozice zobecnit na případ nejvýše spočetně mnoha zdrojů, které na počátku pozorování vypustí požadované koncentrace částic v daných místech trubice. Řešení (16) pak bude tvořeno součtem hodnot hustot N(cij,2Dt) rozdělení pravděpodobnosti násobených plošnými koncentracemi N^/S. 3 Nekonečná trubice s podélným zdrojem Jak se řešení (16) změní, když budou zdroje rozmístěny v nespočetně mnoha bodech trubice? Představme si, že na počátku bude zajištěna koncentrace C0 částic v celém úseku trubice mezi dvěma body a a b, —oo < a < b < oo. Hledáme tedy řešení difuzní rovnice (1) pro x e IR s počáteční podmínkou a < x < b C{x,0) = g{x) = { . (17) x < a nebo x > b I v tomto případě lze díky linearitě difuzní rovnice využít princip superpozice. Nyní však součet hodnot hustot v nespočetně mnoha bodech přejde v integrál 52 0 1.0 Obrázek 6: Koncentrace C (x, t) v mM dle (16) v závislosti na prostorové souřadnici x e IR v cm a na čase ŕ > 0 v sekundách. Na počátku je v místě a: = 0 cm vypuštěno N1 = 10-6 mol částic a v místě a2 = 0,5 cm množství N2 = 0,5 • 10-6 mol částic, jejichž difúzni koeficient je roven D = 10-5 cm2 s-1. Trubice má konstantní průřez o ploše S = 0,1 cm2. hustoty f (x —u, t) rozdělení N(u, 2Dt) v uvedených mezích, r a r b C(x,ŕ) = C0 f(x-u,t)du = -C0\ /(x-u,ŕ)du. (18) Jb Ja Z vlastností spojitých náhodných veličin víme, že určitý integrál hustoty pravděpodobnosti je roven rozdílu hodnot distribuční funkce, řešení difúzni rovnice 53 o 100 40 20 60 80 0 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 x Obrázek 7: Rozložení C(x, t) v mM z Obr. 16 v závislosti na x (v cm) v několika časových okamžicích. s podélným počátečním zdrojem tak lze zapsat také ve tvaru s distribuční funkcí F podle (8). Opět připomeňme rozměry veličin. Koncentrace C(x, t) i C0 mají jednotky mol cm-3, distribuční funkce F je bezrozměrnou veličinou. Pokud bychom počáteční koncentraci chtěli specifikovat pomocí množstvím částic, podobně jako u bodového zdroje, můžeme psát Tento zápis však v praxi není obvyklý, zejména proto, že umožňuje pracovat jen se zdrojem konečné délky (b — a). Průběh řešení je zachycen v grafu v Obr. 8. Rozložení koncentrace v trubici v několika časových okamžicích je pak zobrazeno v Obr. 9. V čase ř = 0 je C(x, t) rovno počáteční podmínce g(x) ve tvaru funkce se skoky v bodech a a b. Pro ř > Oje již rozložení koncentrace hladkou funkcí ve tvaru rozdílu hodnot distribuční funkce regulárního normálního rozdělení pravděpodobnosti. S rostoucím časem se tvar řešení dále více vyhlazuje, neboť roste rozptyl 2Dt. (19) o — (b-a)S' 54 40000 Obrázek 8: Koncentrace C(x, t) v mM dle (19) v závislosti na prostorové souřadnici x e IR v cm a na čase t > 0 v sekundách. Na počátku je mezi body a = —0,5 cm a b = 0,5 cm koncentrace částic rovna C0 = 10~5 mol cm-3. Difuzní koeficient je roven D = 10~5 cm2 s_1. Počáteční úlohu nyní zobecníme na nekonečně dlouhý počáteční zdroj, který je umístěn podél trubice od daného bodu a e IR. Uvažujeme tedy počáteční podmínku ve tvaru Tato skoková funkce je (pro c0 = 1) známá pod názvem Heavisideova funkce (zprava spojitá). Řešení počáteční úlohy (1), (20) lehce obdržíme z řešení (18) limitním přechodem b —> oo, Vidíme, že řešením je C0-násobek hodnoty distribuční funkce F normálního rozdělení pravděpodobnosti se střední hodnotou a a rozptylem 2Dt. (20) (21) 55 10 u -1.5 — t = Os — t = 1 min t = 30 min — t = 1 h -1.0 -0.5 0.0 x 0.5 1.0 1.5 Obrázek 9: Rozložení C (x, ř) v mM z Obr. 8 v závislosti na x (v cm) v několika časových okamžicích. Průběh řešení je pro volbu D = 0,05, c0 = 1, a = 0 zachycen v grafu v Obr. 10. Rozložení koncentrace v trubici v několika časových okamžicích je pak zobrazeno v Obr. 11. Vidíme tvary distribuční funkce normálního rozdělení s rozptylem lineárně rostoucím v čase. Počáteční podmínce t = 0 odpovídající Heavisideova funkce je také distribuční funkcí, a to tzv. singulárního normálního rozdělení pravděpodobnosti (jemu odpovídající hustota je 5-funkce). Již nás nepřekvapí, že s rostoucím časem se počáteční rozložení koncentrace částic v trubici vyhlazuje. 4 Trubice s nulovou koncentrací na okraji Přejdeme nyní k úlohám o trubici s jedním okrajem, který pro zjednodušení zápisu obvykle umisťován do bodu 0. Podobně jako u nekonečně dlouhé trubice lze i tento model s určitou nepřesností použít v případě trubice konečné délky, kdy druhý okraj je buď dostatečně vzdálený nebo koncentrace částic na něm není vnějším prostředím ovlivňována. 56 -1.0 40000 Obrázek 10: Koncentrace C(x, t) v mM dle (21) v závislosti na prostorové souřadnici xeRv cm a na čase t > 0 v sekundách. Parametry mají hodnoty a = 0 cm, C0 = 10"5 mol cm"3, D = 10"5 cm2 s"1. Uvažujeme tedy trubici s jedním okrajem, přičemž počáteční koncentrace částic v ní je popsána funkcí g(x), Na okraji trubice je koncentrace částic stále udržována na nulové hodnotě, tedy tzv. okrajová podmínka je tvaru x > 0. (22) ř > 0. (23) Řešení difuzní rovnice (1) pro x > 0 s počáteční podmínkou (22) a okrajovou podmínkou (23) je obvykle uváděno jako (24) 57 10 u - — t = Os - — t = 1 min - — t = 30 min — t = 1 h 4 - 0 - -1.0 -0.5 0.0 x 0.5 1.0 Obrázek 11: Rozložení C (x, ř) v mM z Obr. 10 v závislosti na x (v cm) v několika časových okamžicích. kde g je tzv. liché rozšíření funkce g, o, —g(—x) x > 0 x = 0 . x < 0 (25) Můžeme si virtuálně představit nekonečnou trubici, v níž je v bodě a zdroj počáteční koncentrace a v bodě —a, tedy symetricky vzhledem k okraji skutečné trubice, je místo, z něhož jsou částice stejnou rychlostí odčerpávány. Tato konfigurace zajistí nulovou koncentraci částic v bodě 0, tedy splnění okrajové podmínky, díváme-li se jen na polovinu trubice. Tvar řešení je principiálně stejný jako řešení (11), ale vystupuje v něm liché rozšíření počáteční podmínky. Aby byl vliv lichého rozšíření na řešení více patrný, rozdělíme (24) na dva integrály a využijeme (25), C(x,ť -f(x — u,t)g(—u)du+ f(x — u,t)g(u)du. 58 V prvním integrálu následně provedeme substituci — u >-> u a zaměníme meze. Tím budou oba integrály přes interval [0, oo) a lze řešení lze přepsat do tvaru f oo C(x,t)=\ [f(x-u,t)-f(x + u,t)]g(u)du. (26) Jo Vidíme, že v integrandu vystupuje rozdíl hodnot dvou hustot normálního rozložení pravděpodobnosti: reálné N(u,2Dt) a zdánlivé, umístěné symetricky vzhledem k okraji trubice, N(—u,2Dt). Právě tento princip zajistí nulovou koncentraci částic na okraji trubice, tedy platnost okrajové podmínky. Uvažujme dále speciální případ Nn g(x) = -±5(x-a), (27) kdy je v čase ř = 0 vypuštěno množství N0 částic v místě a > 0 trubice o ploše průřezu S. Liché rozšíření této počáteční podmínky je zřejmě (např. si načrtněte graf) Nn g(x) = -g-[5(x-a)-5(x + a)] (28) a řešení pak dle (24) či (26) s pomocí (6) nalézáme ve tvaru násobku rozdílu hodnot hustoty normálního rozdělení se střední hodnotou a (poloha reálného zdroje) a rozptylem 2Dt a hustoty normálního rozdělení se střední hodnotou —a (poloha zdánlivého zdroje) a stejným rozptylem, Nn C(x, t) = -±\f(x- a, t) -f (x + a,t)}. (29) Řešení je zobrazeno na grafu v Obr. 12. Rozložení koncentrace v trubici v několika časových okamžicích je zachyceno v Obr. 13. Průběh řešení je podobný řešení pro nekonečnou trubici, zde však máme řešení jen na poloprostoru a navíc je reflektována okrajová podmínka nulové koncentrace. 5 Trubice s podélným zdrojem a nulovou koncentrací na okraji Předchozí počáteční podmínku nyní rozšíříme na podélný zdroj o koncentraci C0 v intervalu mezi body a a b, 0 < a < b < oo, podél trubice. Máme tedy trubici s jedním okrajem, přičemž počáteční koncentrace částic v ní je popsána 59 0.6 0.4 x 0.2 0.0 Obrázek 12: Koncentrace C (x, t) v mM dle (29) v závislosti na prostorové souřadnici x > 0 v cm a na čase ŕ > 0 v sekundách. Hodnoty parametrů jsou a = 0,25 cm, N0 = 10"6 mol, S = 0,1 cm2, D = 10"5 cm2 s-1. funkcí C(x,0) = g(x) CQ, a b Liché rozšíření počáteční podmínky je přitom zřejmě (30) x = < CQ, a < x < b 0, x < — b nebo —ab. (31) -CQ, -b 0 s počáteční podmínkou (30) a okrajovou podmínkou (23) obdržíme dosazením do (24), r —a r b C(x,t)=-C0 f(x-u,t)du + C0\ f(x-u,t)du. J-b Ja 60 V integrálech provedeme substituci v = x — u a oba integrály poté s využitím (8) přepíšeme jako rozdíly distribučních funkcí. Obdržíme výsledek C(x,t) = C0 [-F(x + b,t) + F(x + a,t) + F(x-a,t) -F(x-b,t)] (32) ve tvaru kombinace hodnot distribučních funkcí normálních rozdělení se středními hodnotami —b,—a,a,b a stejným rozptylem 2Dt. Řešení je zobrazeno na grafu v Obr. 14. Rozložení koncentrace v trubici v několika časových okamžicích je zachyceno v Obr. 15. Je patrné vyhlazování prostorového rozložení s rostoucím časem, současně je reflektována okrajová podmínka nulové koncentrace. 61 Obrázek 14: Koncentrace C (x, t) v mM dle (32) v závislosti na prostorové souřadnici x > 0 v cm a na čase ŕ > 0 v sekundách. Podélný počáteční zdroj částic o koncentraci C0 = 10~5 mol cm-3 je umístěn mezi body a = 0,5 cm a b = 1 cm. Difúzni koeficient je roven D = 10~5 cm2 s-1. 62 0.0 0.5 1.0 1.5 2.0 x Obrázek 15: Rozložení C(x, t) v mM z Obr. 14 v závislosti na x > 0 (v cm) v několika časových okamžicích. 63 Řešení rovnice difúze, 2. část Předpoklady Pro studium této výukové jednotky je potřebná znalost • základních principů řešení diferenciálních rovnic • tvaru a odvození parciální diferenciální rovnice difúze Výstupy z učení Studenti se naučí • používat Diracovu delta-funkci • porozumět řešení rovnice difúze v ID • porozumět řešení rovnice difúze při radiální difúzi 1 Trubice s nulovým tokem na okraji V této a následující části prozkoumáme řešení rovnice difúze dC(x,t) d2C(x,t) dt oxÁ za podmínky nulového toku na okraji trubice. Počáteční podmínku přitom stále předpokládáme ve tvaru C(x,0)=g(x), x>0. (2) Nulový tok znamená, že první derivace koncentrace částic vzhledem k prostorové souřadnici je na okraji trubice nulová. V každém čase je tedy koncentrace částic v nějakém malém, ale celém, okolí okraje konstantní, čímž zde nedochází k difúzi. To však neznamená, že koncentrace je také konstantní v čase. Nikoliv, během času se může měnit, ale tak, aby se v celém, malém, okolí okraje měnila stejným způsobem. Okrajová podmínka je tedy matematicky zapsána ve tvaru d C — (0,t)=J(0,t)=0, t>0. (3) dx Nulovost toku na okraji trubice, tzn. platnost (3) lze, podobně tomu bylo jako u nulové koncentrace na okraji, zajistit trikem. Symetricky vzhledem k okraji, 64 tzn. bodu 0, umístíme zdánlivý zdroj(e) stejné koncentrace, jako je zadáno počáteční podmínkou. Počítáme tedy nyní s tzv. sudým rozšířením g funkce g, g(x) = g(\x\), (4) pomocí níž lze řešení difuzní rovnice (1) pro x > 0 s počáteční podmínkou (2) a okrajovou podmínkou (3) zapsat jako f oo C(x,t)= f(x-u,t)g(u)du. (5) J —oo Abychom prozkoumali vliv sudého rozšíření, rozdělíme (5) na dva integrály a využijeme (4), r 0 r oo C(x,t)=\ f(x-u,t)g(\u\)du+\ f(x-u,t)g(u)du. J-oo Jo V prvním integrálu provedeme substituci — u <-> u a zaměníme meze. Tím budou oba integrály přes interval [0, oo) a lze řešení lze přepsat do tvaru f oo C(x,t)= [f(x-u,t)+f(x + u,t)]g(u)du. (6) Jo Vidíme, že v integrandu vystupuje součet hodnot dvou hustot normálního rozložení pravděpodobnosti: reálné N(u,2Dt) a zdánlivé, umístěné symetricky vzhledem k okraji trubice, N(—u, 2D t). Tento princip zajistí nulový tok částic na okraji trubice, tedy platnost okrajové podmínky. Stejně jako úlohách z předchozí kapitoly dále prozkoumáme speciální počáteční podmínku g(x) = -±5(x-a). (7) Sudé rozšíření této počáteční podmínky je zřejmě Nn S(*) = y [5(x-a) + 5(x + a)} (8) a řešení pak dle (5) či (6) s pomocí (6) nalézáme ve tvaru N0 /S-násobku součtu hodnot hustoty normálního rozdělení se střední hodnotou a (poloha reálného zdroje) a rozptylem 2Dt a hustoty normálního rozdělení se střední hodnotou —a (poloha zdánlivého zdroje) a stejným rozptylem, Nn C(x,t) = ^ [f(x-a,t)+f(x + a,t)}. (9) 65 Řešení je zobrazeno v grafu na Obr. 1. Rozložení koncentrace v trubici v několika časových okamžicích je pak zachyceno v Obr. 2. Na průbězích z Obr. 2 si zejména všimněme, že koncentrace na okraji, tzn. v bodě 0, není nulová, dokonce se v čase mění. Jako funkce prostorové proměnné, x, však koncentrace na okraji dosahuje lokálního extrému, první parciální derivace koncentrace vzhledem k x je zde skutečně nulová. Je zde tedy nulový tok, jak je předepsáno okrajovou podmínkou (3). x Obrázek 1: Koncentrace C(x, t) v mM dle (9) v závislosti na prostorové souřadnici x > 0 v cm a na čase ř > 0 v sekundách. Okrajovou podmínkou je stanoven nulový tok na okraji, x = 0, na počátku je v bodě a = 0,25 cm vloženo množství NQ = 10~6 mol částic. Průřez trubice má obsah S = 0,1 cm2, difuzní koeficient má hodnotu D = 10~5 cm2 s_1. 66 100 80 60 - 40 - 20 -■ 0 -■ 0.0 0.2 0.4 x — t = 1 min — t = 5 min t = 30 min — t= lh 0.6 Obrázek 2: Rozložení C (x, t) v mM z Obr. 1 v závislosti na x > 0 (v cm) v několika časových okamžicích. Koncentrace na okraji není obecně nulová, ale vzhledem k prostorové proměnné zde má nulovou derivaci. 2 Trubice s podélným zdrojem a nulovým tokem na okraji Počáteční podmínku bodového zdroje na trubici s nulovým tokem na okraji opět rozšíříme na podélný zdroj o koncentraci C0 v intervalu mezi body a a b, 0 < a < b < oo. Sudé rozšíření této počáteční podmínky, C(x,0) = g(x) C0, a < x < b 0, 0 < x < a nebo x > b (10) je rovne C0, —b < x < —a nebo a < x < b 0, x < —b nebo — a < x < a nebo x > b (H) Po prostudování předchozích kapitoly již čtenář patrně intuitivně očekává tvar řešení difuzní rovnice (1) pro x > 0 s počáteční podmínkou (10) a okrajovou podmínkou (3), C(x,t) = C0[F(x + b,t)-F(x + a,t)+F(x-a,t)-F(x-b,t)], (12) 67 ve tvaru kombinace hodnot distribučních funkcí normálních rozdělení se středními hodnotami —b, —a, a, b a stejným rozptylem 2D t. Rozdíl od (ŘrdI.32) je jen v opačných znaménkách zdánlivého zdroje. Formální výpočet spočívá v dosazení (11) do (5), r —a r b C(x, t) = C0 f(x — u,t)du + C0 f (x — u,t)du, J-b Ja následné substituci v = x — u v obou integrálech a využití (Řrd1.8). Detaily tohoto výpočtu a úpravu na tvar (12) již ponecháváme čtenáři jako cvičení. Řešení C(x, t) je zobrazeno na grafu v Obr. 3, rozložení koncentrace v trubici v několika časových okamžicích je pak zachyceno v Obr. 4. Je vidět postupné vyhlazování prostorového rozložení s rostoucím časem. Současně je patrné dosažení lokálního extrému prostorových průběhů koncentrací v bodě x = 0, tedy nulovost toku na okraji. 3 Trubice s trvalým zdrojem na okraji Dalším, v praxi častým, typem úlohy je difúze v trubici, na jejímž okraji je po celou dobu experimentu udržována předepsaná konstantní koncentrace C0 částic, které jsou na okraji technickým zařízením dodávány, nebo naopak odčerpávány. Řešíme tedy difuzní rovnici (1) pro x > 0 s počáteční podmínkou C(x,0) = 0, x>0, (13) a okrajovou podmínkou C(0,ř) = C0, ř>0. (14) Difuzní rovnice má za těchto podmínek řešení C(x,ř) = 2C0[l-F(x,ř)]=2C0 • X ,2 1 \ * . , 1--, exp--du (15) přičemž F označuje distribuční funkci (Řrd1.8). Výraz v hranatých závorkách v (15) je tedy funkcí přežití a grafem závislosti koncentrace C(x, ř) na x je 2 C0-násobek pravé poloviny grafu funkce přežití normálního rozdělení. 68 i_,_í_._,_f_I_._i_._I_I_I_>_i_1_1_■_1_>_O 2.0 1.5 1.0 0.5 0.0 x Obrázek 3: Koncentrace C (x, t) v mM dle (12) v závislosti na prostorové souřadnici x > 0 v cm a na čase t > 0 v sekundách. Okrajovou podmínkou je stanoven nulový tok na okraji, x = 0, na počátku jsou do části trubice mezi body a = 0,5 cm a b = 1 cm vloženy částice o koncentraci C0 = 10~5 mol cm-3. Průřez trubice má obsah S = 0,1 cm2, difúzni koeficient má hodnotu D = 10-5 cm2 s"1. Průběh koncentrace částic odpovídající řešení této úlohy je zobrazeno na grafu v Obr. 5. Pozorujeme, že na okraji je, v souladu s okrajovou podmínkou, koncentrace částic konstantní v čase. Rozložení koncentrace v trubici v několika časových okamžicích je zachyceno v Obr. 6. Prostorová rozložení koncentrace mají tvar konvexních funkcí klesajících z hodnoty C0 asymptoticky k nule. 69 10 t= Os t = 1 min 8 t = 30 min t= lh 6 4 2 0 0.0 0.5 1.0 1.5 2.0 x Obrázek 4: Rozložení C(x,t) v mM z Obr. 3 v závislosti na x > 0 (v cm) v několika časových okamžicích. Koncentrace na okraji není obecně nulová, ale vzhledem k prostorové proměnné zde má nulovou derivaci. 4 Trubice s bodovým zdrojem a nulovou koncentrací na okra- Jako poslední z úloh na řešení jednodimenzionální rovnice difúze si ukážeme řešení pro trubic konečné délky a tedy se specifikovanými dvěma okrajovými podmínkami. Tato úloha, ač možná z prvního pohledu prakticky nejčastější, je z matematického hlediska mnohem složitější, než všechny předchozí. Současná specifikace obou okrajových podmínek, tedy hodnot koncentrace ve dvou bodech, je značně omezující a řešení takové úlohy je tak možné pouze ve tvaru Fourierovy řady, tedy nekonečné řady goniometrických funkcí. Uvažujeme trubici délky L, 0 < L < oo, jejíž jeden konec pro jednoduchost umístíme do počátku souřadnicového systému. Prostorovou souřadnici tedy uvažujeme v intervalu x e [0,L]. Předpokládáme, že počáteční koncentrace částic v trubici je popsána funkcí g(x), a že koncentrace částic na obou koncích trubice je trvale udržována na nulové jích (16) 70 o Obrázek 5: Koncentrace C(x, t) v mM dle (15) v závislosti na prostorové souřadnici x > 0 v cm a na čase ř > 0 v sekundách. Na okraji (x = 0) je stále udržována koncentrace částic na hodnotě C0 = 10~5 mol cm-3. Difuzní koeficient má hodnotu D = 10~5 cm2 s_1. hodnotě, C(0,ř) = C(L,ř) = 0, ř>0. (17) Difuzní rovnice (1) s počáteční podmínkou (16) a okrajovými podmínkami (17) má řešení . . 2 fknx\ C (x, ř) = - 2^ srn J exp ve tvaru Fourierovy řady, nekonečné řady funkcí, kde jsou za bázové funkce voleny goniometrické funkce (zde konkrétně sinus) s proměnlivými periodami. Pro jednoduchost uvažujme opět speciální případ, kdy na počátku jsou vypuštěno množství N0 částic v místě uprostřed trubice, g(x) = ^ô(x-^, xg[0,L}. (19) 71 k2n2Dt f , . fknu\ , ---— g(u)sin — )du (18) Integrály v (18) jsou pak rovny — o x--sin knu au = — srn knL ~2L { No S > 0, JVo S ' k k k 1,5,9,... 2,4,6,... 3,7,11,... a řešení této speciální úlohy je ve tvaru nekonečné řady -IV'-1 srn C x,ř "(2j-l)7ix" (2j-l)27i2Dř exp L L2 (20) Graf řešení (20) je zobrazen v Obr. 7. Při numerických výpočtech je nutno vzít dostatečně členů nekonečné řady, řádově několik set, aby na grafu koncentrace nebyly patrné artefakty vlnění o vysoké frekvenci v důsledku skládání řešení z goniometrických funkcí. Rozložení koncentrace v trubici v několika časových okamžicích je zachyceno v Obr. 8. Krátce po začátku je rozložení podobné Gaussově křivce s chvosty rovnými nule na okrajích trubice. S narůstajícím časem se tvar postupně mění na tvar vlny jakožto části průběhu goniometrických funkcí. Časové vývoje koncentrace v několika místech trubice jsou 72 zachyceny v Obr. 9. Uprostřed trubice, tzn. v bodě x = L/2 = 0,5, koncentrace exponenciálně klesá, v jiných částech trubice nejdříve roste, později také pozvolna klesá. Nejvyšší koncentrace je v každém časovém okamžiku uprostřed trubice, směrem k okrajům trubice klesá, na okrajích je v souladu s okrajovými podmínkami (17) udržována nulová hodnota. Obrázek 7: Koncentrace C(x, t) v mM dle (20) v závislosti na prostorové souřadnici x e [0,1] v cm a na čase ř > 0 v sekundách v trubici jednotkové délky, L = 1 cm. Okrajové podmínky stanovují nulovou koncentraci na obou okrajích, počáteční podmínkou je dáno množství N0 částic vložených v čase ř = 0 do středu trubice, tedy místa x = L/2 = 0,5 cm. Trubice má průřez o obsahu S = 0,1 cm2, difúzni koeficient má hodnotu D = 10~5 cm2 s_1. 73 0.0 0.2 0.4 0.6 0.8 1.0 x Obrázek 8: Rozložení C (x, ŕ) vmMzObr. 7 v závislosti na x e [0,1] (vcentime-trech) v několika časových okamžicích. Na počátku má rozložení koncentrace přibližně tvar Gaussovy funkce, s rostoucím časem se postupně mění na tvar půlvlny grafy goniometrických funkcí, avšak vždy s nulovými hodnotami na obou okrajích. 5 Radiální difúze z mikropipety Na závěr si ukážeme řešení triviálního případu difúze v třírozměrném prostoru. Mikropipeta je naplněna částicemi fluorescenční látky a její konec je vložen přibližně do středy nádoby s vodou. V čase ř = 0 do vodního prostředí vstříknuto množství N0 fluorescenční látky a ta se difúzí šíří rovnoměrně všemi směry. Sledujeme vývoj koncentrace v čase v okolí tohoto místa. Vzhledem k tomu, že vlastnosti vodního prostředí jsou ve všech směrech stejné, jedná se o sféricky symetrickou úlohu (DlF.38). Řešíme tedy rovnici dC(r,t) 8t D 2 dC(r,t) d2C(r,t) r d r d r2 (21) pro čas ŕ > 0 a proměnnou r > 0, popisující vzdálenost od místa vstříknutí fluorescenční látky. Počáteční podmínka je matematicky zapsána ve tvaru C(r,0)=N053 (r) (22) 74 Obrázek 9: Vývoj koncentrace C(x, t) v mM z Obr. 7 v čase, ř, (v sekundách) v několika místech trubice, x e [0,1] (v centimetrech). V pravé části trubice jsou průběhy symetrické vzhledem ke středu trubice. Na obou okrajích, tzn pro x = 0 cm a x = L = 1 cm, je koncentrace nulová. kde 53 je třírozměrná delta-funkce. Ta je odlišná od jednorozměrné delta-funkce, neboť je definována jako limitní případ třírozměrného rovnoměrného, resp. třírozměrného normálního rozdělení pravděpodobnosti. Pro popis třírozměrné delta-funkce nám však bude stačit následující vyjádření pomocí jednorozměrné delta-funkce, z7rrz z něhož je patrné, že rozměr funkce 5 3 je rovný převrácené hodnotě třetí mocniny jednotky veličiny r. Difuzní rovnice (21) s počáteční podmínkou (22) má řešení tvaru C(r, ř) = N° exp (-—). (24) Jedná se vlastně o hustotu pravděpodobnosti centrovaného třírozměrného normálního rozdělení pravděpodobnosti s diagonální kovarianční maticí a homogenním rozptylem, kde místo proměnných x, y, z vystupuje pouze veličina 75 r = ^/ x2 + y2 + z2. Časový vývoj koncentrace na sférách v několika vzdálenostech od počátku je zobrazen v grafu na Obr. 10. Pozorujeme obecně podobný tvar jako v Obr. 9, ve třírozměrném prostoru je však pokles koncentrace po dosažení maxima rychlejší. o — r = 0,10 — r = 0,15 — r = 0,20 — r = 0,25 50 100 150 200 250 300 350 Obrázek 10: Průběh koncentrace C (r, t) (v M) v třírozměrném prostoru podle (24) v závislosti na čase (v sekundách) na sférách o vzdálenosti r (v centimetrech) od místa, kde bylo na počátku experimentu vstříknuto N0 = 0,1 mol fluorescenční látky o difuzním koeficientu D = 5 • 10~5 cm2 s-1. Jak se průběh koncentrace změní, když látka nebude vstříknuta najednou, ale celé množství N0 bude z mikropipety uvolňováno rovnoměrně, rychlostí N0/10 během intervalu délky ř0? Řešíme tak difuzní rovnici (21) s podmínkou na rychlost l(x, t) přidávání částic v bodě x = 0 a čase ř, 0 < ř < ř0 í(0, ŕ) = \t0 . (25) t>t0 Závislost koncentrace na čase a vzdálenosti od počátku získáme integrováním 76 (24) přes časovou proměnnou v intervalu [0, t0]. Obdržíme výsledek N, C(r,t) = { 2nDrt *—[l-F(r,t)], 0 < t < tr o N, ^ 2 n D r t0 [F(r,t-t0)-F(r,t)], t > t0 (26) v němž vystupují hodnoty distribuční funkce F podle (Řrd1.8). Závislosti C(r, ř na čase jsou vykresleny v Obr. 11 pro několik různých hodnot difuzního koeficientu D. Pro velké hodnoty D koncentrace C(r, ř) velmi rychle reaguje na ukončení vstřikování látky v čase ř0. Pro malé hodnoty D je difúze pomalá, a dosažení čas dosažení maximální koncentrace je několikanásobně delší, než doba ř0. o 1000 2000 t 3000 4000 Obrázek 11: Průběh koncentrace C(r, ř) (v M) v třírozměrném prostoru podle (24) v závislosti na čase (v sekundách) na sféře o vzdálenosti r = 1 mm od místa, kde bylo po dobu ř0 = 400 sekund od počátku experimentu rovnoměrnou rychlostí vstřikováno celkové množství N0 = 0,1 mmol fluorescenční látky. Závislosti jsou uvedeny pro různé hodnoty difuzního koeficientu D (v cm2 s_1). 77 Vybraná témata z teorie informace Předpoklady Pro studium této výukové jednotky je potřebná znalost • základních pojmů z teorie pravděpodobnosti Výstupy z učení Studenti se naučí • rozumět a interpretovat pojem entropie • počítat entropii diskrétních a spojitých rozdělení pravděpodobnosti 1 Úvod do teorie informace Zájemcům o hlubší pochopení principů teorie informace doporučujeme následující standardní literaturu, Ash (1965), Cover (1991), Gallager (1968) a Sveš-nikov (1971). Uvedené zdroje zároveň sloužily jako zdroj textu této kapitoly. Základní problematika teorie informace byla vymezena americkým matematikem Claude Elwood Shannonem (1916-2001) v jeho práci A mathemati-cal theory oj communication z roku 1948, ve které dokázal klíčové teorémy o přenosu a reprezentaci informace. Shannonovou motivací byla snaha exaktně popsat komunikaci, tedy způsob, jak ve zvoleném místě reprodukovat (přesně nebo přibližně) sdělení či signál vyslaný z jiného místa. Teorie informace ukazuje, že existují omezení na efektivitu reprezentace a spolehlivé komunikace sdělení: (a) Reprezentace informace: Sdělením, či signálem, mohou být hodnoty nějaké fyzikální veličiny (elektrického pole, frekvence zvuku, atp.) či symbolické abstrakce (písmena abecedy). Shannon ukázal, že sdělení lze popsat jako náhodný proces a přiřadit mu informační obsah v jednotkách bit (z anglického binary digit). V nejjednodušším případě je sdělení modelováno jako posloupnost nezávislých pozorování diskrétní náhodné veličiny U s konečnou množinou 78 stavů {ui,u2,...,uM} a pravděpodobnostní funkcí, p(ui)=P{U = ui}, i = l,2,...,M. (1) Průměrný informační obsah (v bitech) najeden symbol takového sdělení je definován jako entropie veličiny U, M Význam entropie a zdůvodnění pojmu informační obsah spočívá vtom, že sdělení realizované veličinou U nelze popsat méně než H (U) bity na symbol. Entropie tedy mimo jiné určuje dolní mez na bezeztrátovou kompresi dat, při nižší velikosti dat by již komprese musela být ztrátová. (b) Spolehlivost komunikace: Sdělení je přenášeno mezi zdrojem a cňem (příjemcem informace) tzv. informačním kanálem. Informační kanál je reprezentován např. kovovým drátem v případě elektronické komunikace, atmosférou v případě zvuku, atp. Vlastní přenos informace kanálem je typicky ovlivněn šumem, kdy výstup z kanálu má stochastický (náhodný) charakter. V nejjednoduš-ším případě je diskrétní informační kanál popsán svou vstupní abecedou X = {x1,x2,...,xn}, výstupní abecedou <8f = {yi,y2, ■ ■ -,ym} a podmíněnou pravděpodobností f (YíM = V {Y = yt\X = xk], i = l,...,m; k = l,...,n, (3) tedy pravděpodobností, že yř byl přijatý signál za podmínky, že x^ byl signál vyslaný. Kritická je zejména situace, kdy dvě různá vstupní sdělení nelze po přenosu rozlišit. Před Shannonovou teorií bylo široce přijímáno, že stochastickým (a tedy nespolehlivým) informačním kanálem nelze přenášet informaci s libovolně malou tolerancí na chybu přenosu, tedy, že současně s snižováním hranice na pravděpodobnost chyby při přenosu bude klesat i informace najeden symbol („jedno užití" kanálu). Shannon však dokázal, že existuje způsob jak informačním kanálem přenášet nenulovou informaci na symbol i když pravděpodobnost chyby je (2) 79 libovolně malá. Toto množství informace se nazývá kapacita kanálu, C, a je definováno jako C = maxVVp(xfc)/(yř|xfc)log: f(yiM p(yd (4) kde p(yi) = 2fc=iP(xfc)/(yilxfc) a maximum je bráno přes všechna možná pravděpodobnostní rozdělení na vstupní abecedě 3C'. Zdroj, U —Zdrojový en kodér Kanálový cnkodcr x = (xu, . .,xN) Kanál, X -> Y y = (/!.....yN) Cíl, V Zdrojový dekodér Kanálový dekodér Obrázek 1: Schéma přenosu informace dle C. E. Shannona. Celý systém komunikující sdělení ze zdroje do cňe je znázorněn na Obr. 1. Sdělení realizované náhodnou veličinou U je nejprve zdrojovým enkodérem komprimováno, v praxi často do podoby binárního řetězce. Tento binární řetězec je kanálovým enkodérem převeden na vstup do informačního kanálu. Sekvence N vstupních symbolů {xl5..., Xjv} (tzv. blokový kanálový kód) je kanálem zobrazena na N výstupních symbolů, přičemž úkolem kanálového dekodéru je určit vstupní sekvenci s minimální pravděpodobností chyby. Výstup je v cíli dekomprimován zdrojovým dekodérem. Teoreticky lze v případě H(U) < C takovým systémem přenášet H(U) bitů na symbol a přitom zajistit, aby pravděpodobnost že zdrojové a cílové sdělení se liší byla libovolně malá. V praxi nebývá teoretických limitů dosaženo z důvodů přílišné výpočetní složitosti, zejména kanálového dekódování. 80 2 Axiomatické odvození Shannonovy entropie Vzorec pro Shannonovu entropii (2) lze odvodit na základě minimálně omezujících a intuitivních axiomů. Entropie zde hraje roli míry nejistoty, náhodnosti či predikovatelnosti pozorování generovaných diskrétní náhodnou veličinou X. Poznamenejme, že entropie je definována rovnicí (2) i pokud M = oo. Předpokládejme, že X je diskrétní náhodná veličina s konečnou množinou stavů {x1,x2, ■ ■ ■, xM} a příslušnými pravděpodobnostmi {p\,p2, ■ ■ ■ ,Pm)-Předpokládejme, že pro funkci H (X) = H(p1,.. .,Pm) platí následující axiomy: (i) f (M) = H {—,...,— ] je rostoucí funkcí vzhledem k M, v ' \M M J (ii) pro pozorování {yl5..., yL} náhodné veličiny Y nezávislé na X platí f(ML)=f(M)+f(L), (iii) tzv. seskupovací axiom: pro r = 1,..., M platí H(pl,...,pM) =H(p1-\-----Vpr,pr+l^-----hpM) + + (p1 + ---+Pr)H[-^,...,-p—] + \Zji=iPi Zjí=iPíJ + (Pr+1 + ---+PM)H(-^,...,-^\, (5) V2ji=r+1 Pi Zji=r+1 Pi. (iv) H(p,l—p) je funkce spojitá v proměnné p,pe [0,1]. Pak jediná funkce splňující všechny čtyři výše uvedené axiomy je tvaru M H(X) =~cJ]Pi l°SaPu (6) i=l kde c > 0 a pro základ logaritmu platí a > 1. Ještě než přistoupíme k vlastnímu důkazu, všimněme si podmínek (i)-(iv). Heuristicky se často (i) označuje jako extensivita, (ii) + (iii) jako aditivita a (iv) 81 jako spojitost entropie. Uvažujme náhodný experiment se dvěma nezávislými náhodnými veličinami, X = {xl5...,xM} aY = {y1,... ,yL}, Současné pozorování realizace (X,Y) pak nabývá jedné z ML stejně pravděpodobných možností, tedy průměrná nejistota v předpovědi výsledku sdruženého experimentu je f (ML). V případě že známe výsledek pozorováními, celá „nejistota" či náhodnost v experimentuje dána pouze výsledkem pozorování Y. Z důvodu nezávislosti X a Y pak očekáváme, že průměrná nejistota v predikci výsledku sdruženého experimentu (X,Y) mínus nejistota získaná ze znalosti výsledku X, musí dát nejistotu experimentu zahrnujícího jen pozorování Y. Z tohoto pohledu je axiom (ii) zcela přirozený. Nyní k vlastnímu důkazu tvrzení. Z (ii) plyne / (1 • 1) = /(l)+/(l), takže musí platit / (1) =0. Tento výsledek intuitivně koresponduje s představou, že míra náhodnosti v experimentu s diskrétní náhodnou veličinou s rozdělením Pi = V{X = Xi} = 1 s jediným možným výsledkem je nulová. Zároveň máme / (Mfc) = k f (M), takže intuitivní volbou pro / je logaritmická funkce, což dokážeme následujícím způsobem. Nechť M > 1 je celé číslo a r je libovolné přirozené číslo. Pak musí existovat číslo k e N0 takové, že Mk<2r 1, dostáváme fc< logn(2) <(fc + l) r loga(M) r Abychom ukázali že /(x) = loga(x), všimneme si, že oba podíly v (8) a (9) leží ve stejných mezích a tedy pro jejich rozdíl platí /(2) loga(2 /(M) loga(M; < -. (10) r 82 Jelikož r je libovolné, můžeme položit r —> oo, a tedy oba dva členy na levé straně nerovnosti (10) jsou si rovny, čímž dostáváme /(M) = clogaM, (11) kde c = / (2)/ loga 2 > 0, protože / (l) = 0a / má být rostoucí funkce. Dále, nechť p = r /s je racionální číslo dané nějakými přirozenými čísly r < s. Z axiomu (iii) plyne f(s) H[-, — ) + -f(r) + —f(s-r). (12) \ s s I s s Poznamenejme, že platí m -.....- 1 = m r/s r/s HÍ> ■.....' M- (13) r s — r Při označení p = - můžeme psát - = 1 —p, a podle (ii) a (11), lze vý- s s raz (12) zapsat ve tvaru c logas =H(p,l-p) +cplogar + c(l-p)loga(s-r), (14) z čehož plyne H(p,l~p) =-c[plogar-logas + (l-p)loga(s-r)]. (15) Odečtením a přičtením výrazu p logas k pravé straně (15) obdržíme H (p, 1 -p) = -c [p loga p + (1 -p) loga(l -p)]. (16) Rovnice (16) pak platí pro všechna p e [0,1] podle axiomu (iv). Rovnici (6) pak dokážeme matematickou indukcí. Její platnost pro pro M = 1 a M = 2 jsme si již ukázali výše. Pro M > 2 využijeme axiom (iii), pro rozdělení na pravděpodobnosti q = '^L^ Pí a pM, H(p1,...,pM)=H(q,pM)+qH(^-,...,^-)+pMH(l). (17) 83 Předpokládejme, že rovnice (6) platí pro všechna přirozená čísla menší nebo rovna M — 1, a dokážeme její platnost pro M. Z (17) máme M-l H(pi,...,pM) = -c[qlogaq+pMlogapM]-cq V — loga —+pM-0, (18) éíg g čímž je platnost (6) indukcí dokázána, neboť M-l M-l Q Y\ — loga" = Y! Pí l°ga Pí" 9 l°ga T/*- (19) éíg g éí Pro základ logaritmu se často, zejména ve spojení s informačními technologiemi, volí a = 2 a c = 1, jednotkou entropie je pak bit. Lze se setkat i se základem a = 10, kdy je jednotkou entropie tzv. dit. V matematice má výsadní postavení přirozený logaritmus loge = ln, tedy logaritmus o základu a = e, kdy se jednotka entropie nazývá nat. Všimněme si ještě, že axiom (ii) plyne z (iii), protože H(l/ml,..., 1/ml) můžeme seskupit do m segmentů délky l, takže 2f=i Pí = a následně /(ml) =H( ) = jy ' \ml mlj f 1 n vil ( m m \ p, . p, . = H +> —H i-- =/ m)+/ l . (20) Vm mj p*m Vml ml) jy ' jy ' Přesněji, axiomy (ii) + (iii) lze nahradit jediným axiomem H(Pi, ■ ■ -,Pri> Prj+lj • ■■>Pr2' Pru + lj • • • >Prk) = (21) pro k skupin, každou o rk prvcích (položme r0 = 0), přičemž součet pravděpodobností v i-té skupině označíme qř = 2jLr- j+i^j- Na závěr poznamenejme, že funkce (6) není definována pro události xk s nulovou pravděpodobností, V{X = xk} = 0. Prostým dosazením dostáváme neurčitý výraz „0 • oo". Funkci však můžeme snadno dodefinovat limitou, neboť z 1'Hospitalova pravidla plyne lim plogap = 0, (22) p->0+ tedy události s nulovou pravděpodobností k entropii nepřispívají. 84 3 Kvantizace a entropie spojité náhodné veličiny Shannonova entropie ve vzorci (2) je definována pouze pro diskrétní náhodné veličiny. Uvažujme nyní spojitou náhodnou veličinu X s Riemannovsky inter-grabilní hustotou pravděpodobnosti f(x), a najděme entropii její jednoduché kvantizace. Obor hodnot náhodné veličiny X rozdělíme na ekvidistantní intervaly [(i-l)A, i A), přičemž předpokládáme, že hustota / (x) je na těchto intervalech spojitá. Z Lagran-geovy věty o střední hodnotě plyne, že v každém dělicím intervalu existuje takové xř, že r i A f(xi)A= /(x)dx. (23) J(i-1)A Definujeme kvantizovanou diskrétní náhodnou veličinu XA tak, že XA = xt pro (i-l)A 0, diferenciální entropie h(X) může nabývat i hodnot záporných a tedy není ekvivalentem Shannonovy entropie. Pro A —> 0+ však můžeme položit h(X) & H(XA) + log2 A. Zvolíme-li A = 2~n, pak lze tvrdit, že entropie n-bitové kvantizace X je přibližně h(X) + n. Úkoly 1. Mějme alternativní rozložení s pravděpodobnostní funkcí P {X = 0} = p a P{I = 1} = 1 — p, kde 0 < p < 1. Určete Shannonovu entropii jako funkci parametru p. 2. Mějme náhodnou veličinu X s Gaussovým (normálním) rozdělením s nulovou střední hodnotou a rozptylem a2, tzn. s hustotou Určete diferenciální entropii h(X) náhodné veličinyX jakožto funkci směrodatné odchylky cr > 0. 86 3. Střelec střílí do dvou terčů, do velkého terče dvakrát a do malého terče třikrát. Pravděpodobnosti zásahu v každém pokusu jsou rovny 1/2 při střelbě do velkého a 1/3 při střelbě do malého terče. Do kterého z terčů je výsledek střelby, daný počtem zásahů, méně neurčitý? Který terč má menší směrodatnou odchylku počtu zásahů? Spočítejte i střední hodnoty obou veličin. 87 88 4. Pokračujme v příkladu 3 se střelbou do dvou terčů. Za každý zásah do malého terče se vyplácí odměna 100 Kč, za každý zásah do velkého terče pak 50 Kč. Náhodná veličina U udává celkovou výhru střelce na malém terči, náhodná veličina U pak celkovou výhru na velkém terči. Spočítejte střední hodnoty obou náhodných veličin. Která z veličin U, V je méně neurčitá? Která z veličin U, V má menší směrodatnou odchylku? 89 5. Spočítejte entropii diskrétního rovnoměrného rozdělení pravděpodobnosti X na množině {xl5 ...,xN}. 6. Spočítejte diferenciální entropii exponenciálního rozdělení pravděpodobnosti X jako funkci intenzity X > 0. 7. Spočítejte diferenciální entropii rovnoměrného spojitého rozdělení pravděpodobnosti X na intervalu [a, b], oo < a < b < oo. 90 91 Matematické modely neuronu Předpoklady Pro studium této výukové jednotky je potřebná znalost • základních matematicko-logických funkcí Výstupy z učení Studenti se naučí • vysvětlit pojem logického neuronu • rozumět jednoduchým matematickým modelům neuronu • popsat dynamiku modelu reálného neuronu 1 Logický neuron Nejznámějším matematickým modelem neuronu je tzv. logický neuron nazývaný též McCullochův-Pittsův logický neuron. Americký neurofyziolog a kybernetik Warren Sturgis McCulloch (1898-1969) a Walter Harry Pitts Jr. (1923-1969), zabývající se matematickou logikou a početními neurovědami, tento model poprvé popsali roku 1943 v článku A logical calculus of the ideas immanent in nervous activity v odborném časopise Bulletin of Mathematical Biophysics. Obrázek 1: Logický neuron s prahem 6, se čtyřmi excitačními vstupy a jedním inhibičním vstupem. McCullochův-Pittsův logický neuron bývá schematicky znázorňován tak, jak je ukázáno na Obr. 1. Ve své nejjednodušší podobě můžeme takový neuron chápat jako zařízení s několika vstupy a jedním výstupem. Vstupy jsou kresleny 92 jako spojnice s šipkami směřujícími do neuronu, výstup je reprezentován jednou spojnicí vycházející ze symbolu neuronu. Výstup se však může dělit a může být následně napojen do většího počtu (jiných) neuronů. Vstupy mohou být buď excitační, kreslené klasickou šipkou, anebo inhibiční, u nichž se ve schématu doplňuje kolmá čárka. Na výstupu se přitom objeví signál právě tehdy, když se současně signály objeví na určitém počtu vstupů. Tento minimální počet aktivních vstupů je nazýván práh a obvykle značen 9. V nejjednodušším modelu může být práh 9 rovný jakémukoliv nezápornému číslu. Logický neuron je idealizací reálného neuronu a zachovává tak nejdůle-žitější charakteristické vlastnosti skutečného neuronu. Předpokládáme, že logický neuron může změnit svůj stav jen v diskrétních ekvidistantně vzdálených časových okamžicích. Výstup z jednoho neuronu se na vstupu napojených neuronů objeví v bezprostředně následujícím časovém kroku. Síť logických neuronů se chová synchronizované, okamžiky možné změny stavu jsou pro všechny neurony v síti totožné. Biologové a neurofyziologové mohou tento model neuronu kritizovat pro přílišně až nereálné zjednodušení, zvlášť kvůli výše popsané časové závislosti. Na druhou stranu, obrovská výhoda logického neuronu je právě jeho jednoduchost. Ta často umožňuje analýzu očekávaného chování neuronové sítě, která by při zahrnutí jiných fyziologických parametrů neuronů byla příliš složitá. Logický neuron zachovává podstatné rysy chování reálného neuronu. Daní za jednoduchost modelu je zanedbání méně podstatných charakteristik reálného neuronu. Jedná se o obvyklý postup matematického modelování reálného světa. Často si musíme realitu zjednodušovat, abychom dokázali matematický model schůdně analyzovat a modelovat. Na provedená zjednodušení však přitom nezapomínáme, a jejich možný vliv na chování modelu diskutujeme. Nyní podáme formální definici logického neuronu. (i) Logický neuron se může nacházet v jednom ze dvou stavů, které nazveme aktivní a neaktivní. (ii) Logický neuron má jeden výstup, který může být současně napojen do jednoho nebo více jiných neuronů nebo také zpět do stejného neuronu. V každém okamžiku je na všech těchto napojeních stejný výstup. (iii) Logický neuron má celkem ne + n{ vstupů, z nichž je ne excitačních a n{ 93 inhibičních. Počty excitačních i inhibičních vstupů mohou nabývat hodnot jakéhokoli nezáporných celých čísel. (iv) Logický neuron má práh 9. V nejjednodušší podobě požadujeme 9 e N0, obecně stačí 9 > 0. (v) Čas je kvantizován, stav logického neuronu se může změnit pouze v diskrétních časových okamžicích (krocích), ř = k A, k = 0,1,____Stav logického neuronu zůstává nezměněn v časovém intervalu [fc A, (fc + 1)A). Délka časového kroku A je stejná pro všechny neurony v síti, neuronová síť je synchronizovaná. (vi) Daný vstup je aktivní v čase (fc + 1) A právě tehdy, když neuron, jehož výstupem uvažované spojení je, byl aktivní v čase k A. Označujeme ATe(ř) počet aktivních excitačních vstupů a iVj(t) počet aktivních inhibičních vstupů v čase ř. Zřejmě platí ATe(ř) < ne, iVj(t) < /i; v každém čase ř. (vii) Logický neuron je aktivní v čase ř = k A, tedy přesněji ve všech časech ř e [k A, (fc + 1)A), právě když Ne(kAt)-Ni(kAt)>B. (1) Konstanta

0 charakterizuje požadovaný poměr mezi počtem aktivních excitačních a aktivních inhibičních vstupů pro dosažení aktivního výstupu neuronu. Často je 1, tj. pokud je aktivní alespoň jeden z jeho (excitačních) vstupů. Pokud by tento neuron měl prahovou hodnotu 9 = 0, může být aktivní, přestože by nebyl aktivní žádný z jeho vstupů. Takový neuron je někdy nazýván spontánně aktivní. Obrázek 3: Logický neuron s prahem 1, jedním excitačním vstupem a jedním inhibičním vstupem. Na Obr. 3 vidíme neuron s ne = 1 excitačním a ri; = 1 inhibičním vstupem a prahovou hodnotou 9 = 1. Předpokládejme = 1, jedním externím vstupem (in) a výstupem (out). Chování této sítě v závislosti na časovou posloupnost aktivity vstupu uvádí Tabulka 2. Je zřejmé, že výstup je aktivní pouze tehdy, když v předchozích (alespoň) dvou po sobě jdoucích krocích byl vstup aktivní. Tato jednoduchá neuronová síť tak může představovat velice primitivní mozek, jehož funkcí je reagovat na bezprostředně opakovanou aktivitu na vstupu. 95 Ni Ne-> 0 0 0 0 1 -1 1 0 1 1 1 0 Tabulka 1: Funkce logického neuronu z Obr. 3 pro = 1. Jaký bude vývoj aktivity této sítě, když ji začneme pozorovat v některém z možných stavů? Využijeme zápisu pomocí binárního zápisu. Ciferný zápis čísla ABC ve dvojkové soustavě kóduje po řadě stav neuronu A, B, a C. Např. 101 značí stav, kdy neurony A a C jsou aktivní, neuron B je neaktivní. Uvedená síť se může nacházet v kterémkoliv z 23 = 8 možných stavů. Přechodový diagram (ověření přenecháváme čtenáři jako cvičení) 101-► 111-► 110-► 10 0 000 o 001 popisuje vývoj aktivity sítě v čase, každá šipka znázorňuje změnu stavu za krok délky Ař. Vidíme, že pokud síť začne pracovat v kterémkoliv stavu jiném než 000 (všechny neurony neaktivní), po několika krocích dospěje k oscilačnímu chování, kdy se pravidelně budou měnit stavy 001 a 100, tedy aktivní bude střídavě vždy pouze neuron A a C, neuron B bude neaktivní. 97 2 Logický neuron a logické funkce McCulloch a Pitts ve svém článku z roku 1943 upozornili také na to, že logické neurony definované dle (i)-(vii) mohou jednoznačně reprezentovat funkce matematické logiky. Tento poznatek tehdy vzbudil značnou pozornost a předpokládalo se, že neurony i mozek fungují právě na principu binárních logických funkcí. Početní neurovědy na současném stupni poznání již toto chování považují za nereálné, především z důvodu obrovského počtu neuronů a spojení a s tím spojené těžkopádnosti takového mozku. Spojení logického neuronu a funkcí matematické logiky nicméně zůstává zajímavou částí matematického modelování. Aktivní vstup či výstup reprezentuje logickou hodnotu true (pravda) označovanou číslem 1, neaktivní vstup či výstup znamená logickou hodnotu falše (nepravda) označovanou číslem 0. Aktivita výstupu logického neuronu je dána pomocí logické funkce vstupů neuronu v bezprostředně předcházejícím časovém kroku. Ukážeme alespoň pár příkladů reprezentace základních logických funkcí, složitější funkce k procvičení čtenář nalezne v úkolech na konci kapitoly. X y xVy x A y -■x x =>y 0 0 0 0 1 1 0 i 1 0 1 1 1 0 1 0 0 0 1 i 1 1 0 1 Tabulka 3: Tabulka logických operací, zleva doprava logický součet, logický součin, logická negace a implikace. Logický součet (or, nebo) je funkce, která nabývá hodnoty 1 právě když alespoň jeden ze dvou vstupů má hodnotu 1. Logický součin (and, a) je funkce, která nabývá hodnoty 1 právě když oba vstupy mají současně hodnotu 1. Logická negace (noř, ne) je unární funkce, která mění vstupní hodnotu 1 na výstupní 0 a naopak. Implikace je binární logickou funkcí funkcí, na jejímž výstupu je hodnota 0 právě tehdy, když první vstup má hodnotu 0 a druhý hodnotu 1. Definice těchto funkcí shrnuje Tabulka 3. Není těžké ověřit, že uvedené logické funkce lze modelovat pomocí logických neuronů podle následujících schémat: logický součet na Obr. 7, logický součin na Obr. 8, negaci podle 98 Obr. 9 a implikaci podle Obr. 10. Obrázek 7: Logický součet (or, nebo). Obrázek 8: Logický součin (and, a). x —K -■x Obrázek 9: Logická negace (noř, ne). Obrázek 10: Logická implikace. 3 Model reálného neuronu Logický neuron má řadu vlastností reálného neuronu. Jednou z těch pro reálné neurony typických, které však logický neuron postrádá, je spojitost času, tedy to, že neuron může svým výstupem reagovat na vstup v kterémkoliv časovém okamžiku, nejen v časových krocích délky Ař. V této kapitole se přidržíme mj. popisu modelu neuronu z úvodu monografie Gerstner & Kistler (2002). U reálného neuronu rozlišujeme tři funkčně odlišné části, analogické logickému neuronu. Vstupy jsou označovány dendrity, tělo neuronu jako soma a výstup je nazýván axon. Soma je považována za výpočetní jednotku, která na 99 základě určité, obecně nelineární, funkce převádí signál z dendritů na axon. Pokud celkový součet signálů na dendritech překročí prahovou hodnotu 9, na axonu je vygenerována výstupní signál. Spojení mezi neurony je nazýváno synapse. Neuron, který signál vysílá, je nazýván presynaptický, neuron přijímací signál je označován jako postsynap-tický. Typický neuron v mozkové kůře obratlovců je napojen na více než 104 postsynaptických neuronů. Neuronální signál je tvořen krátkými elektrickými impulsy. Každý takový impuls je nazýván akční potenciál {action potential), častěji spike. Posloupnost špiků tvořící signál je označována jako spike train. Špiky mají typicky amplitudu cca 100 mV a trvání 1 až 2 ms. Špiky generované stejným neuronem mají vždy stejný tvar a tento tvar zůstává neměnný i při šíření axonem. Protože všechny špiky generované neuronem vypadají stejně, nenesou informaci. Předpokládá se tak, že jednotkou informace přenášené mezi neurony je jeden akční potenciál (spike) a signál je nějakým způsobem zakódována ve spike trainu, konkrétně v rozložení po sobě jdoucích špiků. Akční potenciály ve spike trainu jsou obvykle jasně oddělené. I při velmi silném vstupu neuron nedokáže generovat další akční potenciál ve chvíli generování prvního špiku a po určitou dobu po něm. Minimální časový odstup mezi dvěma po sobě následujícími akčními potenciály je nazýván refraktorní periodou. Vliv akčního potenciálu na postsynaptický neuron je experimentálně zjišťován pomocí elektrody umístěné blízko soma nebo axonu neuronu. Tou se měří potenciál V, tedy rozdíl elektrického napětí mezi vnitřkem neuronu a jeho blízkým okolím; přesněji se nazývá membránový potenciál. Bez přítomnosti akčních potenciálů na vstupech neuronu je membránový potenciál udržován na konstantní, tzv. klidové hodnotě UT = —70 mV. Při příchodu akčního potenciálu na vstup neuronu se membránový potenciál změní a poté se vrátí na klidovou hodnotu. Přitom akční potenciál na excitačním vstupu vyvolává kladnou změnu membránového potenciálu, akční potenciál na inhibičním vstupu vyvolává zápornou změnu membránového potenciálu. Dynamika neuronu je popsána pomocí závislosti membránového potenciálu V(t) na čase, ř, a jeho reakce na příchod akčního potenciálu na některý ze 100 vstupů. Tzv. spike response model (SRM) předpokládá časovou závislost membránového potenciálu ve tvaru V(r) = Vr + ]T]reř(ř-ř, t-tu) + rj(t-ř). O) i j Konstanta VT je klidový potenciál. Funkce eř je změna vyvolaná příchozími akčními potenciály na vstupech od presynaptického neuronu i, tzv. postsynap-tický potenciál (PSP). Pokud je vstup od presynaptického neuronu i excitační, jedná se o excitační postsynaptický potenciál (EPSP), pokud je vstup inhibiční, označuje se jako inhibiční postsynaptický potenciál (IPSP) a příslušné eř má zápornou hodnotu. Ve vyjádření závislosti V(t) vystupuje součet postsynaptic-kých potenciálů ze všech neuronů, jejichž axony jsou synapsemi napojeny na dendrity sledovaného neuronu. Čas řřj označuje okamžik j-tého akčního potenciálu na axonu neuronu i, čas ř pak představuje okamžik posledního akčního potenciálu generovaného sledovaným neuronem do času ř, ř = max{řj ; ij < ŕ}. (4) Zde ij značí okamžik j-tého akčního potenciálu sledovaného neuronu. Jako poslední člen ve V(t) vystupuje funkce 17, která modeluje nelineární chování membránového potenciálu, vyvolané překročením prahové hodnoty 9 a generováním akčního potenciálu. Excitační postsynaptický potenciál eř(ř — i,s) je jako funkce proměnné s obvykle uvažován ve tvaru unimodální funkce, která v okamžiku generování akčního potenciálu presynaptického neuronu roste z nulové hodnoty, nabývá maxima a poté asymptoticky klesá k nulové hodnotě. Analogicky to platí pro IPSP —eř. Tvar závislosti PSP je přitom pro daný neuron neměnný. První proměnná, ř — ř je čas od posledního akčního potenciálu, PSP bývá neklesající funkcí této proměnné, pro velmi malé hodnoty t— i PSP nabývá malých hodnot (v absolutní hodnotě). Pokud je příchozích akčních potenciálů málo, je chování membránového potenciálu lineární, V(t) je součet PSP všech presynaptických neuronů. Pokud počet akčních potenciálů na dendritech sledovaného neuronu překročí jistou mez, chování V(ŕ) se změní na nelineární. Pokud hodnota membránového potenciálu při svém růstu překročí prahovou hodnotu 9 v čase ř = i, tedy za 101 podmínky V (t ) = 9 a —r^>0, (5) dř dojde k vygenerování akčního potenciálu. Chování V(t) se přitom zásadně změní, hlavní roli převezme poslední člen v (3). Dojde k rychlému nárůstu potenciálu, který je následován poklesem až pod klidovou hodnotu VT s následným exponenciálním přiblížením k VT. Často je tvar 17 popisován funkcí 0 < ŕ < Ař V(t) = \~ . (6) -r)0exp(--J, ř>Ař Zde vystupující konstanty jsou dané typem neuronu, Ař určuje délku intervalu, během níž je hodnota 17 a tím i V(t) velmi vysoká, větší než 9, n0 udává minimální potenciál po vygenerování akčního potenciálu, t pak určuje rychlost přiblížení ke klidovému potenciálu. Výška skoku l/Ar se často uvádí kolem 100 mV. V jednoduchých modelech se často klade Ař —> 0+, čímž dojde k vytvoření krátkého impulsu nekonečné hodnoty, který matematicky popisujeme Diracovou delta-funkcí, viz (Rrd.2), rj(ř) = 5(ř)-rj0 exp^-^j. (7) Po dobu Ař a také po jistou dobu, kdy je 17 (ř — ř < 9, je vliv PSP epsilorii na potenciál V (ř) velmi malý čímž je je modelována refraktorní perioda neuronu. V reálných neuronech je bývá délka refraktorní periody rovna 2 až 4 ms. Prahová hodnota neuronů bývá asi o 20 mV vyšší než klidový potenciál VT, tzn. okolo hodnoty 9 = —50 mV. Amplituda PSP (v absolutní hodnotě) je však typicky jen 1 až 2 mV. Ke splnění podmínky (5) a vygenerování akčního potenciálů je tak potřeba příchod 20 až 50 akčních potenciálů na vstupy sledovaného neuronu v dostatečně krátkém časovém intervalu. Více či méně komplexním matematických modelů reálného neuronů je více. Pro jejich složitost jim zde nebudeme věnovat pozornost, zmiňme však alespoň jejich názvy. Tzv. leaky integrate-and-fire model je založen na popisu pomocí elektrického obvodu a diferenciálních rovnic pro membránový potenciál a proud. Jiným modelem je Hodgkinův-Huxleyův model, který k popisu membránového potenciálů využívá soustavu dvou diferenciálních rovnic. 102 Tuto kapitolu ukončíme velmi jednoduchým modelem, který je postaven na vlastnostech závislosti V(t). Tento model předpokládá, že chování membránového potenciálu je v každém čase dáno vždy právě jednou z následujících možností: (a) Pokud presynaptický neuron generoval akční potenciál v čase ř, potom se hodnota potenciálu V v čase ř + Ař zvýší o konstantu 17. Konstanta 17 zde představuje amplitudu EPSP (pokud je kladná), resp. amplitudu IPSP (je-li záporná). Časový posun Ař reprezentuje zpoždění mezi okamžikem generování akčního potenciálu presynaptickým neuronem a okamžikem příchodu maxima EPSP, resp. minima IPSP, na dendrit sledovaného neuronu. (b) V čase i, pro nějž platí (5), dojde k vygenerování akčního potenciálu. Po dobu délky R > 0 po jeho vygenerování, tedy v intervalu t G [i, i + Ŕ), klademe V(t) = 00. To odpovídá volbě (7) v SRM modelu. Dále klademe V (i = Ŕ) =VT. Konstanta R tedy odpovídá refraktorní periodě neuronu. (c) Ve všech ostatních případech nepopsaných pravidly (a) a (b) membránový potenciál exponenciálně klesá, přesněji jeho časový vývoj je popsán diferenciální rovnicí kde a > Oje průměrná časové konstanta rychlosti konvergence PSP k nule. Uvádí se, že l/a má hodnotu jednotek milisekund. 1. Jaká logická funkce vstupů x,y,z je na výstupu w reprezentována logickým neuronem na Obr. 11? Uvažujte