Parameter estimation je metoda využívaná často při syntéze (vytváření) modelu. Pokud máme model s neznámými parametry (podle formalismu, ve kterém je zapsaný - například v modelu daném diferenciálními rovnicemi mohou být některé konstanty neznámé), umožňuje nám na základě chování (experimentu) určit jejich hodnotu (s určitou přesností).
V biologických modelech může být tato metoda důležitá. Například modely metabolismu vychází z reakčních schémat. Pro daný zkoumaný jev (například signální dráhu) je možné reakční schemata získat (například z literatury nebo výzkumu). Problematické jsou ale reakční konstanty, které určují rychlost reakce, a tím výrazně ovlivňují chování modelu. Přitom neexistují žádné metody jak z organismu získat reakční konstanty vhodné pro daný model. Proto se používá paremeter estimation.
Tento článek zkoumá chování parameter estimation na modelu fotosyntézy popsaném v (Lazár, 2009a) který je zároveň veřejně dostupný v rámci projektu e-photosynthesis. Využívá (a zároveň testuje) přitom nástroj Copasi. Hlavním účelem tohoto článku je ukázat fungování parameter estimation na složitém reálném modelu.
Pojem parameter estimation může označovat celou řadu podobných metod. Pro potřeby tohoto článku je použita definice podle nástroje Copasi.
Vstupem pro parameter estimation je model s neznámými parametry a experiment ve formě časové řady. Výstupem pak je přiřazení hodnot těmto parametrům tak, že simulace modelu se nejvíce blíží experimentu. Zřejmě se bere v úvahu pouze jedna počáteční podmínka.
Pro parameter estimation je zásadní pojem účelové funkce (objective function). Je to funkce, která pro každé přiřazení hodnot parametrů vyjadřuje, jak moc se liší simulace modelu s danými hodnotami parametrů liší od experimentu. Výstupem experimentu je časová řada obsahující měření nějaké veličiny popisující chování systému (například teplota nebo koncentrace nějakého proteinu). Podobně výstup simulace pro dané hodnoty parametrů je časová řada s hodnotami této veličiny. Účelová funkce pak je statistická veličina daná rozdíly hodnot veličiny z experimentu a ze simulace v daných časových bodech (přesná definice: manuál Copasi, sekce 3.1). Pokud experiment sleduje hodnoty více veličin, musí se z účelových funkcí pro každou z hodnot vytvořit (vážený) průměr.
Čím nižší má účelová funkce hodnotu, tím více se simulace pro dané hodnoty parametrů blíží experimentu. V optimálním případě, kdy je její hodnota nulová, simulace experimentu odpovídá (což ale v realitě nemusí platit, protože experiment většinou obsahuje pouze několik málo měření hodnoty veličiny). Parameter estimation pak spočívá v hledání optimální hodnoty účelové funkce.
Pokud máme n neznámých parametrů, je účelová funkce funkcí z n rozměrů. Na tomto n-rozměrném prostoru vytváří funkce řadu "vyboulenin" a "kopců", daných jejími lokálními minimy. Parameter estimation hledá globální minimum. Tento parametrický prostor se často nazývá evoluční krajina (evolutionary landscape).
|
| Příklad evoluční krajiny pro dva neznámé parametry |
Problém hledání globálního minima na evoluční krajině je známý jako numerická optimalizace. Existuje celá řada heuristik pro jeho řešení. Zde je třeba zdůraznit slovo heuristika. Úplné prohledání parametrického prostoru je časově exponenciální vůči počtu parametrů, takže, s výjimkou malých prostorů o málo rozměrech, není uskutečnitelné v reálném čase. Heuristiky pro numerickou optimalizaci, a tím i parameter estimation, mají mnoho společného s magií. Studie ukazují, že na praxi fungují, ale není řečeno za jakých okolností a to, že výsledná hodnota bude skutečně optimální, není zaručeno.
Většina heuristik pro numerickou optimalizaci je založena na lokálním prohledávání. Při něm se začne v nějakém bodě parametrického prostoru a pak se po něm pohybuje tak, aby se snižovala hodnota účelové funkce. Protože se pohybujeme pouze o malou vzdálenost, je hledání časově zvládnutelné. Problémem ale může být to, že se heuristika zastaví v lokálním minimu (ze kterého už nemůže dál), které však není globální. Různé heuristiky tento problém řeší různě. Na FI MU se některými těmito heuristikami okrajově zabývá předmět Algoritmika pro těžké problémy.
Copasi implementuje větší množství heuristik, v následujícím testování budeme používat tyto (manuál Copasi):
| Simulated annealing (simulované žíhání) | Simulované žíhání řeší problém uváznutí v lokálním minimu tím, že umožňuje s určitou pravděpodobností dovoluje posunout se do bodu s větší hodnotou účelové funkce. Tato pravděpodobnost se postupně snižuje. |
| Particle Swarm (optimalizace hejnem částic) | V parametrickém prostoru se zároveň pohybuje více částic; jejich rychlost a směr pohybu se mění podle úspěšnosti minulého pohybu. |
| Evolutionary Programming (Evoluční programování) | Evoluční algoritmus pracující ve více iteracích. V parametrickém prostoru nachází určitá populace jedinců. Při každé iteraci se vybere polovina s vyšší hodnotou účelové funkce, z každého vybraného jedince se vytvoří kopie a ta "zmutuje", tedy se náhodně pozmění. |
| Newtonova metoda | Využívá při pohybu v parametrickém prostoru derivaci funkce. Je určena spíše pro jednoduché funkce a je zde použita jen jednou, pro srovnání. |
Je nutné poznamenat, že těchto heuristik vykazuje dobré výsledky pro "pěkné" evoluční krajiny, tedy ty, které mají málo lokálních minim a jejich globální minimum je výrazné. Naopak problematické jsou "hrbolaté" krajiny, které mají hodně lokálních minim s podobnou hodnotou. Ukazuje se, že většina biologických modelů vytváří právě složité, "hrbolaté" krajiny.
Protože většina výše uvedených heuristik je založena na lokálním prohledávání, zřejmě bude důležité, ve kterém bodě prohledávání začíná.
V dalším textu se zaměříme především na evoluční programování, protože je v Copasi přednastavené jako metoda pro parameter estimation.
Model fotosyntézy byl zvolen proto, že se jedná o relativně velký model reálného jevu, který už je odladěný (tzn. nejsou v něm žádné neznámé parametry), což znamená, že můžeme porovnávat výsledky parameter estimation se známými hodnotami parametrů. Obsahuje celkem 45 metabolitů, 73 reakcí (některých reverzibilních, tedy celkem 112 diferenciálních rovnic) a 30 konstant, které lze brát jako parametry. Model byl převeden do formátu Copasi (ke stažení).
Při reálných experimentech se sleduje fluorescenční vzestup (FLR), jehož křivka má velmi charakteristický tvar (tak zvaný O-J-I-P přechod). Více informací o měření FLR je možné najít v (Lazár, 2009b). V daném modelu je veličina určující fluorescenční vzestup označena jako globální hodnota Funq. Protože experimentální data nebyla k dispozici, jako základ experimentu byla použita simulace už odladěného modelu (ke stažení). Posléze byly označený některé hodnoty parametrů za neznámé a provedena na nich parameter estimation. Tento přístup je specifický tím, že známe hodnoty, kterých by měly parametry nabývat a výsledek parameter estimation s nimi můžeme srovnat.
 |
| Fluorescenční vzestup (simulace) - O-J-I-P přechod |
Použitý model má svá specifika. Konkrétně, skládá se z velkého množství diferenciálních rovnic (viz výše), tedy je pravděpodobné, že výsledná evoluční krajina bude velmi "hrbolatá". Dále, O-J-I-P přechod lze sledovat pouze v logaritmickém časovém měřítku, zatímco měření (a simulace) je v čase lineárním. Z těchto důvodů to nemusí být úplně vhodný model pro aplikaci parameter estimation. Na druhou stranu, jedná se o příklad reálného modelu, pro jehož vytváření může být parameter estimation nutné použít.
Při aplikaci parameter estimation na výše uvedený model nám šlo především o praktickou použitelnost této metody. Konkrétně, otázkou je, jestli dává rozumně správné výsledky v rozumném čase. Konkrétně, měřítky jsou:
Co se trvání výpočtu týče, všechny výpočty byly spouštěny v 64b linuxovém systému běžícím na 1.8GHz dvoujádře s 2.5 GB RAM.
Kromě parametrů, jejichž hodnoty se mají určit a vstupního experimentu umožňuje Copasi pro každý parametr nastavit omezení jeho hodnoty, tvořené intervalem. To může být výhodné, pokud známe přibližnou hodnotu parametru a chceme ji určit přesně. Dále je možné nastavit hodnotu, z které se začíná průzkum parametrického prostoru. Za ni je rovněž možné vzít náhodnou hodnotu, čehož budeme využívat.
Každá z heuristiky implementovaných v Copasi má rovněž svá nastavení. Ty by měly ovlivnit kvalitu výsledků, potažmo i délku výpočtu. Není však zřejmé, jaká nastavení použít pro dosažení požadovaných výsledků, proto budeme používat pouze základní nastavení. Výjimkou je pouze počet iterací evolučního programování, který (společně s velikostí populace, tedy počtu částic -- viz výše) jasně udává délku výpočtu. Podobně metoda particle swarm má omezení na maximální délku výpočtu.
Jak nastavit úlohu parameter estimation v Copasi je popsáno v přiloženém návodu.
Na začátku zkusíme hodně naivní přístup. Zkusíme parameter estimation na všech parametrech a náhodných počátečních hodnotách. Výsledky jsou většinou nevalné. U metod evolučního programování a particle swarm to může být způsobeno na počátku zadaným omezením doby výpočtu, ale pokud dobu výpočtu prodlužujeme, velmi často nám vyjde výsledek, kdy většina parametrů je blízká nule a průběh FLR se podobá původnímu pouze v tom, že má vzestupnou tendenci. Podivné je pak chování v případě, že počáteční hodnoty parametrů jsou nastaveny na jejich původní hodnoty. Dalo by se předpokládat, že v těchto hodnotách by měla být hodnota účelové funkce nulová a heuristika by se tak hned na začátku měla zastavit. Není tomu ale tak. Nabízí se dvě možnosti, proč by tak tomu mohlo být:
Zřejmě bude třeba hodnoty parametrů nějakým způsobem omezit. Pro účely následujících testů byla provedena parameter estimation opět na všech parametrech, přičemž jejich hodnoty byly omezeny na interval ±50 % původní hodnoty. Nastavení metod bylo ponecháno na přednastaveném, a počáteční hodnota byla opět zvolena náhodně. Parameter estimation byla spuštěna za použití každé z výše uvedených metod. Zvláštní pozornost byla věnována evolučnímu programování, kde byla spuštěna s velikostmi populace 20, 200, 2000 a 20000 jako zkouška toho, jak velikost populace ovlivní výsledek.
Poznámka: Průběh parameter estimation pomocí Newtonovy metody byl přerušen poté, co se po téměř 150 hodinách výpočtu nepodařilo snížit hodnotu účelové funkce pod 0,1. Je tedy vysoce pravděpodobné, že Newtonova metoda je pro reálné modely nevhodná.
| metoda | Simulated annealing | Particle swarm | EP (20) | EP (200) | EP (2000) | EP (20000) |
|---|---|---|---|---|---|---|
| trvání výpočtu | 44 h | 1 h 30 m | 4 m | 37 m | 6 h 20 m | 61 h |
| účelová funkce | 8,66e-8 | 1,66e-7 | 1,58223e-4 | 1,20e-5 | 6,91e-6 | 6,24e-7 |
| Trvání výpočtu paramater estimation na všech parametrech a omezení na ±50 % a výsledné hodnoty účelových funkcí | ||||||
Z výsledků je těžké něco vyčíst. U evolučního programování zřejmě existuje přímá úměra mezi velikostí populace a délkou chování. Což je ovšem důležitější, se zvyšující populací se snižuje hodnota účelové funkce, což by mělo znamenat lepší výsledek. Zajímavá je také účinnost heuristiky particle swarm, která dosáhla poměrně malé hodnoty účelové funkce za velmi rozumný čas. Ostatní heuristiky pro podobnou hodnotu potřebovaly až o dva řády více času. Na obrázku (viz níže) je srovnání FLR původního modelu (tmavší) s výsledky parameter estimation.
| Původní | Simulated annealing | Particle swarm | EP (20) | EP (200) | EP (2000) | EP (20000) |
| Porovnání FLR původního modelu a výsledků parameter estimation | ||||||
Většinou má křivka získaná pomocí parameter estimation velmi podobný tvar. Výjimkami jsou výsledky metody particle swarm a evolučního programování s populací 200, které mají o jeden inflexní bod navíc. Zvlášť markantní je rozdíl mezi evolučním programováním s populacemi 20 a 200. Zřejmě menší hodnota účelové funkce nemusí nutně znamenat "lepší" tvar křivky. To, že mezi původním modelem a výsledky parameter estimation jsou rozdíly je zřetelnější na simulacích některých vnitřních proměnných modelu (viz níže).
| Původní | Simulated annealing | Particle swarm | EP (20) | EP (200) | EP (2000) | EP (20000) |
| Porovnání vývoje vnitřních proměnných modelu a výsledků parameter estimation | ||||||
Pokud se podíváme na hodnoty parametrů, zjistíme, že jsou (v rámci omezení omezení ±50 %) poměrně náhodně rozmístěné a vůbec neodpovídají hodnotám v původním modelu. To však také může znamenat, že použitý model je robustní.
 |
| Hodnoty parametrů (relativně k hodnotám v původnímu modelu) |
Pro další testování byl počet neznámých parametrů snížen na šest (kf_ct, kb_ct, kf_(HC), kb_(HC), kf_F a kb_F), aby se snížil počet rozměrů parametrického prostoru, čímž by se měla úloha zjednodušit. Testovala se pouze metoda evolučního programování, s nastavením populace na 250. Hodnoty parametrů byly omezeny, podobně jako v předchozí úloze. Vzniklo tak pět úloh s následujícími omezeními hodnot parametrů vůči původní hodnotě: ±25 %, ±50 %, ±100 %, ±200 % a bez omezení. Každá z těchto pěti úloh byla provedena desetkrát, přičemž pokaždé se začínalo z jiné náhodné hodnoty parametrů.
| omezení | ±25 % | ±50 % | ±100 % | ±200 % | žádné |
|---|---|---|---|---|---|
| účelová funkce | 3,7e-4 ±0,1 ‰ | 2,6e-4 ±0,8 ‰ | 1,4e-4 * | 1,1e-4 ±2,2 % | 5.6e-4 ±128 % |
| doba trvání | 41 m | 44 m | 46 m | 1 h | 1 h 25 m |
| Průměrná doba trvání a hodnota účelové funkce s relativní standardní odchylkou (* značí zanedbatelnou odchylku). | |||||
Na výše uvedených hodnotách je obzvláště zajímavé postupné snižování účelové funkce omezení na ±200 % původní hodnoty. To, že kvalita heuristiky degraduje, pokud se neomezí interval na hodnotách parametrů, lze očekávat. Za pozornost také stojí rozptyl hodnot parametrů získaných z jednotlivých výpočtů. Až do omezení na ±200 se vyskytují na kraji intervalu daného omezením (viz níže).
 |  |  |
| Relativní rozptyl parametrů pro omezení ±25 % | Relativní rozptyl parametrů pro omezení ±50 % | Relativní rozptyl parametrů pro omezení ±100 % |
 |  | Poznámka: Záporné hodnoty parametrů nemají smysl. |
| Relativní rozptyl parametrů pro omezení ±200 % | Relativní rozptyl parametrů bez omezení |
Zřejmě heuristika vnímá evoluční krajinu jako klesající k okrajům daným omezením hodnoty parametrů a nedokáže najít předpokládané minimum v bodě jejich původních hodnot (pokud není parametrický prostor nějak zkreslený). Zajímavá je také fakt, že zatímco pro omezení na ±50 % vykazují první dva parametry rozptyl, při omezení na ±100 % jsou shodně "natlačeny" k okraji, což by odpovídalo evoluční krajině klesající s tím, jak se vzdaluje od původních hodnot parametrů.
Přitom platí, že rozptyl na omezeních na ±25 %, ±50 % a ±100 % nevytváří žádnou viditelnou změnu na průběhu FLR a vnitřních proměnných (které se opticky shodují s jejich průběhy u původního modelu). U omezení na ± 200 % už jsou patrné změny průběhů pro různé výsledné hodnoty parametrů, podobně je tomu pokud se rozsah hodnot nijak neomezí.
Pro následující testování byl počet parametrů omezen na jeden (k_L1). Testování bylo přitom provedeno stejně jako minule: desetkrát pro každé z pěti různých druhů omezení. Zajímavé bylo, že nezávisle na omezení a počáteční hodnotě byl výsledek vždy přesně stejný, ačkoli jiný než původní hodnota (představoval jejích 75 %). Se zvětšujícím se omezujícím intervalem se dala vysledovat pouze mírně zvyšující se doba trvání (která se ale vždy pohybovala mezi 40 a 45 minutami). Křivka FLR přitom nelze od původní okem rozlišit.
Zřejmě omezení na jeden parametr zásadním způsobem zjednodušuje úlohu. Na druhou stranu, zdá se, že vybraná vstupní data určitým způsobem zkreslují výsledek.
Dále byla testována metoda particle swarm na šesti parametrech (kf_ct, kb_ct, kf_(HC), kb_(HC), kf_F, kb_F). Z časových důvodů byl maximální počet iterací nastaven na 500. Podobně jako u evolučního programování bylo spuštěno deset úloh pro každé z následujících omezení hodnot parametrů (relativně vůči původním hodnotám): ±50 %, ±100 % a ±200 %. Podle manuálu Copasi má kromě časového omezení metoda particle swarm ještě nastavení přesnosti hodnoty účelové funkce, při které se výpočet ukončí. Toho však nikdy nebylo dosaženo (a to ani s přednastaveným počtem iterací 2000). Když nebylo nastaveno žádné omezení, parameter estimation vždy skončilo chybou.
| omezení | ±50 % | ±100 % | ±200 % |
|---|---|---|---|
| účelová funkce | 2.6e-4 | 1.5e-4 | 1.1e-4 |
| doba trvání | 23 m | 20 m | 29 m |
| Průměrná doba trvání a hodnota účelové funkce. | |||
Z výsledků je patrná především snižující se hodnota účelové funkce se zvětšujícím se intervalem, na které jsou hodnoty parametrů omezeny, podobně jako v případě testování evolučního programování na stejných parametrech (viz výše). Zásadní je také, že pro různé úlohy při stejném omezení vyšla vždy stejná hodnota účelové funkce. Jak je vidět na relativním rozptylu hodnot parametrů (viz níže), výsledky úlohy byly stejné. Významnou výjimku tvoří úlohy na omezení ±100 %. Ty zároveň skončily před vyčerpáním počtu iterací, a to aniž by dosáhly nastavené přesnosti. Lze předpokládat, že v průběhu došlo k chybě. Jinak se metoda chová deterministicky.
 |  |  |
| Relativní rozptyl parametrů pro omezení ±50 % | Relativní rozptyl parametrů pro omezení ±100 % | Relativní rozptyl parametrů pro omezení ±200 % |
Z rozptylu hodnot parametrů lze vyvodit podobné závěry jako v případě testování evolučního programování na stejných parametrech (viz výše), tedy že z pohledu particle swarm má evoluční krajina směrem od globálního minima sestupnou [sic!] tendenci. Když se podíváme na rozdíl mezi FLR křivkou skutečného modelu a výsledku (viz níže), zjistíme, že s rozšiřujícím se intervalem omezení se rozdíl opticky zvětšuje, a to přesto, že účelová funkce se snížila.
| Původní | ±50 % | ±100 % | ±200 % |
| Porovnání FLR původního modelu a výsledků parameter estimation | |||
Z předešlých testů vyplývá jedna zásadní skutečnost pro parameter estimation: Závislost mezi hodnotou účelové funkce, podobností křivek a hodnotami parametrů není triviální. Nižší hodnota účelové funkce nemusí nutně znamenat větší podobnost křivek. Podobně je zajímavá skutečnost, že přestože experimentální data měla zdroj ve stejném modelu s určitými hodnotami parametrů, nepodařilo se těchto hodnot pomocí parametrů dokázat. To samozřejmě není důležité, pokud "správné" hodnoty parametrů předem neznáme, ale přesto to stojí za pozornost. Ve výsledku je pro výsledek parameter estimation zřejmý velký význam experimentálních dat, tedy toho, jak se účelová funkce počítá. Účelová funkce definována pomocí vzdáleností hodnot veličiny experimentu a simulace nemusí tedy příliš vypovídat o tvaru křivky. Může tak být zajímavé zaobírat se i derivací (první a druhou) veličiny a nebo přístupem na bázi temporálních logik. Je také možné, že tato skutečnost je způsobena "statistickou kompresí", kdy z velkého množství vzdáleností je složeno jedno číslo.
Další zajímavou skutečností je nutnost omezit rozsah možných hodnot parametrů, aby parameter estimation došla k vhodnému výsledku. To předpokládá nějakou další znalost o hodnotách parametrů. Míra této nutnosti přitom závisí na počtu parametrů: pokud je parametrů větší počet, s velkou pravděpodobností bude třeba omezit jejich hodnoty na menší interval, zatímco pro jeden parametr vykazuje parameter estimation dobrý výsledek. Určit, jak přesně hodnoty parametrů omezit nemusí být jednoduché, především proto, že parameter estimation vždy vrací nějaký výsledek, u kterého ale není jasné, jak moc dobrý je (ačkoli hodnota účelové funkce může být určitým ukazatelem). Úplně přímočaré není ani nastavení parametrů heuristik používaných k parameter estimation. Například u časově omezených heuristik nemůže být přesně určena kvalita výsledku. Problémem parameter estimation na parametrickém prostoru o větším počtu rozměrů je pak časová složitost, která efektivně znemožňuje interaktivní práci s nastavením vstupu.
Parameter estimation je tedy použitelná i na tak komplikovaných modelech, jako je model fotosyntézy, nemůže být však používána zcela bez rozmyslu a bez dalších informací o hodnotách parametrů. Její výsledky přitom můžou mít různou kvalitu. Délka výpočtu přitom souvisí s požadovanou kvalitou výpočtu a v případě velkého množství neznámých parametrů může být poměrně velká.
