Kapitola 1 Úvod

1.1 Exponenciální funkce

 

 

(1.1)

     

                 

Derivace je úměrná funkční hodnotě.

Řešením je exponenciální funkce:

 

(1.2)

 

(1.3)

 

tedy

(1.4)

 

protože:

 

pak

(1.5)

 

tato rovnice bude platit pokut hodnota e se v bodě 0 bude limitně blížit:

(1.6)

 

což lze snadno dokázat dosazením:

 

Proto tedy při h=1/10 a h=1/100

(1.7)

 

 

 

 

 

1.2 Taylorovy řady

 

Princip spočívá v tom, že známe-li hodnotu funkce f(t) v bodě t=a můžeme její funkční hodnotu v bodě t=a+h aproximovat pomocí vztahu:

(1.8)

 

Princip této aproximace spočívá v tom, že funkci f(t) nahradíme polynomiální aproximací: k hodnotě funkce v bodě a (tj. f(a))[1] budeme přičítat polynomiální aproximační funkci:

(1.9)

 

kde : h=t-a

 

tedy:

 

uvažujme pro jednodnoduchost aproximaci polynomem třetího stupně:

 

 

Pak po dosazení

 

Tedy:

 

 

Předpokládejme, že aproximační funkce g(h)=g(t-a) bude mít v bodě t=a stejné derivace jako funkce f(t). Porovnáním derivací obou funkcí pak vypočteme hodnoty koeficientů b a c.

 

 První derivace:

 

 

Při t=a:

(1.9)

 

Druhá derivace:

 

Při t=a:

(1.10)

Proto tedy:

Po dosazení koeficientů:

do (viz 1.9):

 

 

dostaneme:

(1.12)

V (1.12) uvažujeme jen dva členy, v Taylorově rozvoji lze obecně n-tý člen vyjádřit vztahem:

 

(1.13)

Obecné řešení můžeme pak vyjádřit:

(1.13a)

Pro exponenciální funkci f(t)=et můžeme fodnotu vypočítat pro bod t=0 kde e0=1 . Pro okolí bodu t=0 můžeme aproximovat:

 

Protože u exponenciální funkce jsou její derivace rovny funkční hodnotě:

Protože e0=1 :

(1.19)

Při t=1 dostaneme podle 1.14 e=2.7083, což je jen 0.5% správné hodnoty.

 

1.3 Imaginární exponenty

Komplexní číslo:

z=x+yi

je vhodné vyjádřit též v Gausově tvaru:

z=r(cos(i)+i sin(i))

kde r je modul komplexního čísla z resp. jeho absolutní hodnota (dá se vyjádřit jako délka vektoru v Gausove rovině). Eulerův vzorec vyjadřuje obecnou exponenciální funkci:

(1.15)

Podle Taylorova vzorce můžeme vyjádřit aproximaci.

 

Vzpomeňme, že:

 

Aproximace reálné části dle Taylorova rozvoje:

Aproximace imaginární části:

¨

 

Tedy výsledná aproximace:

(1.15)

 

1.4 Náhrada ekvivalentními rovnicemi

 

Příklad:

(1.17)

můžeme nahradit soustavou rovnic:

(1.17a)

 

druhou rovnici převedeme na diferenciální:

 

Protože

 

 

 

a tedy výsledná soustava rovnic bude:

(1.18)

 

 

Počáteční podmínka pro y(0): podle 1.17 při t=0 y=1.

 

Řešení soustavy rovnic v Simulinku (viz Example1_8.mdl) ukazuje následující obrázek:

 

Obrázek 1.1 Řešení rovnice 1.17 v Simulinku.

Kapitola 2 Diferenciální rovnice prvného stupně

2.1 Základní diferenciální rovnice prvního stupně

 

Obrázek 2.1 Kompartmentové schéma exponenciálního odtoku.

 

Nejjednodušší diferenciální rovnice popisuje exponenciální  odtok nějakých entit z kompartmentu, případně jeho exponenciální růst (proč exponenciální bude uvedeno za chvíli). Základním principem je, že výtok (případně vtok) F je úměrný množství entit v kompartmentu (velikosti kompartmentu).Označíme-li velikost kompartmentu pomocí proměnné x, pak:

(2.0a)

 

Rychlost změny velikosti kompartmentu (tj. množství entit v kompartmentu) je rozdíl mezi všemi přítoky do kompartmentu a všemi oddtoky z kompartmentu. Máme li pouze přítok (exponenciální růst), pak:

(2.0b)

V případě oddtoku (viz obr. 2.1) je znaménko toku záporné

((2.0c)

Dosadíme-li 2.0c do 2.0a, pak:

(2.1)

Řešením je exponenciální rovnice:

(2.1a)

O správnosti se můžeme snadno přesvědčit. Protože derivace této funkce (používám zde symbolický manipulátor Matlabu - dřinu strojům!):

 

syms('A','a','t')

diff(A*exp(a*t))

 

ans =

A*a*exp(a*t))

 

Je tedy:

(2.1b)

Substitucí 2.1a a 2.1b do 2.1 dostaneme:

(2.1c)

Po vykrácení dostaneme:

(2.1d)

Po substituci do 2.1a:¨

(2.1e)

 

Konstanta A může mít libovolnou hodnotu. Abychom zabezpečili jedinečné řešení, zvolíme ji jako počáteční podmínku při t=0 bude – význam bude možno interpretovat např.jako množství enttit v komparmentu (tj. množství látky apod.) v čase 0.

 

Dosádíme-li  do 2.1e dostaneme výsledné řešení:

(2.2)

Homogenní diferenciální rovnice obsahují na pravé straně neznámou funkci x(t) (např. 2.1)[2]. Nehomogenní rovnice obsahují na pravé straně ještě další člen, nezávislý na x(t), např.:

(2.3)

V kompartmentovém vyjádření tomu odpovídá přidání přítoku, závislého na externí funkci S(t) a časové konstante  - viz následující obrázek:

 

 

Obrázek 2.2 Kompartmentové zobrazení kinetiky popsané nehomogenní diferenciální rovnicí 2.3.

 

Máme tedy odtok z kompartmentu závislý na velikosti kompartmentu x a časové konstantě  a zároveň přítok do kompartmentu  závislý na externí funkci S(t) a časové konstantě . Protože jak přítok tak odtok jsou závislé na stejné časové konstantě , můžeme kompartmentové schéma nahradit ekvivalentním schématem s jedním  tokem .

 

Obrázek 2.3 Ekvivalentní kompartmenové schéma odpovídající obrázku 2.2.

 

Tok F směřuje do kompartmentu – bude-li mít kladnou hodnotu bude přítokem, v případě, že bude mít zápornou hodnotu bude de fakto odttokem. Numerické řešení příslušné diferenciální rovnice bude záviset na počáteční hodnotě kompartmentu v čase t=0. Označme si počáteční hodnotu x(0)=A. Kompartmentové schéma na předchozích dvou brázcích můžeme dekomponovat na schéma s dvěma kompartmenty: , jeden má v počátečním stavu obecně nenulovou hodnotu  a druhý je čase t=0  prázdný, tj. . Chování prvního kompartmentu můžeme popsat klasickou rovnicí:

za předpokladu, že v nulovém čase může mít obecně nenulovou hodnotu A. Chování druhého kompartmentu, který je v čase nula prázdný popíšeme:

 

 

Obrázek 2.4 Kompartment x z předchozího obrázku 2.2 je možné dekomponovat jako součet dvou kompartmentů x1 a x2, kompartment x1 zobrazuje kinetiku oddtoku počáteční hodnoty kompartmentu x a kompartment x2, který byl v čase t=0 prázdný mění svoji hodnotu v závislosti na externí funkci S(t). Jeho kinetika je závislá na celé historii funkce S(t) od času 0 a na časové konstantě , nikoli však na počáteční hodnotě kompartmentu x.

 

 

 

 

Obrázek 2.5 Kompartmentové schéma ekvivalentní předchozímu obrázku, toky Fin a Fout jsou soustředěny do jednoho sumárního toku Fx2=Fin-Fout.

 

Protože  kompartment x je součtem kompartmentů x1 a x2 pak I rychlost změn jeho velikosti bude součtem derivací kompartmentů x1 a x2:

 

Pro kompartment x1 platí již známá exponenciální kinetika:

 

A je počáteční hodnota kompartmentu x1, resp. i celého kompartmentu x - kompartment x2 je na počátku prázdný, proto platí x(0)=x1(0). Kompartment x2, který byl v čase t=0 prázdný, mění svoji hodnotu v závislosti na časové konstantě  a na průběhu externí funkce S(t). Jak si záhy dokážeme, kinetika kompartmentu x2 může být vyjádřena pomocí vztahu:

Na rozdíl od kinetiky kompartmenty x1 je zde místo počáteční hodnoty funkce H(t), která reprezentuje historii průběhu funkce S(t) od počátečního času t=0 do aktuálního okamžiku t.

Protože x=x1+x2, řešení rovnice 2.3 (jak se záhy dosazením přesvědčíme) bude mít formu:

(2.4)

Do rovnice 2.3:

 

dosadíme rovnici 2.4 a pak se pokusíme se vyjádřit funkci H(t). Dosadíme si nejdříve x z rovnice 2.4 do pravé strany rovnice 2.3, pak dosadíme x do levé strany rovnice 2.3 a prvoedeme derivaci (derivujme rovnici 2.4 podle času). Nakonec porovnáme obě strany a z výsledného výrazu vyjádříme derivaci funkce H(t) podél času. Využíváme opět symbolický manipulátor Matlabu:

 

clear

syms ('A','t','tau','H','x','levaStrana','pravaStrana','vyraz')

 pravaStrana=simple(subs('1/tau*(-x+S(t))','x','A*exp(-t/tau)+H(t)*exp(-t/tau)'))

levaStrana=simple(diff('A*exp(-t/tau)+H(t)*exp(-t/tau)'))

vyraz=simple(pravaStrana-levaStrana)

pretty(vyraz)   

 

pravaStrana =

-(A*exp(-t/tau)+H(t)*exp(-t/tau)-S(t))/tau

levaStrana =

-exp(-t/tau)*(A-diff(H(t),t)*tau+H(t))/tau

vyraz =

1/tau*S(t)-diff(H(t),t)/exp(t/tau)

 

                                       d

                                       -- H(t)

                                S(t)   dt

                                ---- - --------

                                tau         t

                                       exp(---)

                                           tau  

 

Pravá strana rovnice 2.4 po dosazení:

Levá strana rovnice 2.4 po dosazení a derivování:

Rozdíl mezi pravou a levou stranou musí být nulový, vypočteme a zjednodušíme výraz:

Vyjádříme derivaci funkce H(t):

Vyjádříme funkci H(t) integrací. V tomto integrálu funkci integrujeme podle pomocné proměnné t' od t'=0 do t'=t:

(2.5)

Výsledný výraz pro funkci H(t) dosadíme do původní rovnice 2.4 a dostaneme:

(2.5a)

Exponentu můžeme "zahnat" do integrálu a x(t) vyjádřit ekvivalentním výrazem:

(2.5b)

Dostali jsme tzv. konvoluční integrál. Konvoluce zachycuje vliv minulých hodnot funkce S(t) na aktuální hodnotu x v čase t. Tento vliv exponencielně klesá vzhledem k minulosti: výraz  v rovnici 2.5b při bude .

Máme tedy odvozenu následující větu:

 

 

Je vhodné poznamenat, že místo zdlouhavého odvozování je možno dostat výsledek přímým řešením diferenciální rovnice  symbolickým manitpulátorem Matlabu (dřinu strojům!):

clear

syms ('A','t','tau','x','S')

dsolve('Dx=1/tau*(-x+S(t))','x(0)=A')

  

 

ans =

(Int(1/tau*S(z1)*exp(1/tau*z1),z1 = 0 .. t)+A)*exp(-1/tau*t)  

 

Příklad na využití věty 1:

Pro čas v milisekundách vyřešte rovnici:

pro počáteční podmínky x(0)=0.

 

Podle věty 1 můžeme psát:

Vypočteme integrál a nakreslíme graf funkce:

 

clear

syms ('A','t','h','x')

A=0;

x=A*exp(-t/10)+(1/10*int(exp(-(t-h)/10)*50,h,0,t))

ezplot(x,[0,40])

grid

  

 

x =

50-50*exp(-1/10*t)

 

  

Výsledek je tedy:

Pro počáteční podmínku x(0)=70 dostaneme:

clear

syms ('A','t','h','x')

A=70;

x=A*exp(-t/10)+(1/10*int(exp(-(t-h)/10)*50,h,0,t))

ezplot(x,[0,40])

grid  

 

x =

20*exp(-1/10*t)+50

 

Při x(0)=A=70 je výsledek:

 

Přímé řešení v MATLABU:

Řešení rovnice:

Pro  počáteční podmínku x(0)=0:

 

 

clear

syms('t','h','x');

x=dsolve('Dx=1/10*(-x+50)','x(0)=0')

ezplot(x,[0,40]);

grid on

  

 

x =

50-50*exp(-1/10*t)

  

 

 

 

Pro  počáteční podmínku x(0)=70:

 

 

clear

syms('t','h','x');

x=dsolve('Dx=1/10*(-x+50)','x(0)=70')

ezplot(x,[0,40]);

grid on;

 

 

x =

50+20*exp(-1/10*t)

  

 

 

 

 

Důležité pro zapamatování je to, že rovnice podobného typu je možné řešit analyticky (pomocí symbolického manipulátoru Matlabu)

 

 

 

 

Další příklad:

Pro počáteční podmínky x(0)=0 vyřešte rovnici:

Využíváme-li symbolický manipulátor Matlabu vůbec nemusíme ručně dosazovat do rovnice podle věty 1. Celou práci je možno svěřit stroji a diferenciální rovnici, kterou chceme vyřešit můžeme v zadání přímo formulovat:

clear

syms ('t','x')

x=dsolve('Dx=1/20*(-x+40*exp(-t/20))','x(0)=0')

ezplot(x,[0,140])

grid 

 

x =

2*exp(-1/20*t)*t

 

Výsledek:

Při počátečních podmínce x(0)=40 dostaneme:

 

clear

syms ('t','x')

x=dsolve('Dx=1/20*(-x+40*exp(-t/20))','x(0)=40')

ezplot(x,[0,140])

grid 

 

x =

2*exp(-1/20*t)*t+40*exp(-1/20*t)

 

Výsledek:

2.2 Kaskády diferenciálních rovnic prvního řádu

Příklad kaskády:¨

(2.6)

Kaskádu je možné si představit jako propojené kompartmenty (viz následující obrázek 2.6).

 

Obrázek 2.6 Kaskáda propojených kompartmentů – používá se pro modelování postupně se šířících signálů s časovým zpožděním. V daném případě chování kaskády závisí na koeficientu zesílení k a časové konstantě , které jsou v daném případě společné pro všechny kompartmenty.

 

Pro řešení kaskádových rovnic je možné využít větu 1. Ukažme si řešení pro počáteční podmínky x(0)=1, y(0)=y(0)=0.

Řešení první rovnice je triviální:

Po dosazení x do druhé rovnice kaskády dostaneme:

(2.7)

Pro řešení použijeme větu 1 a dostaneme:

(2.8)

Dosadíme y do třetí rovnice pro z:

 

Nakonec použijeme opět větu 1 a vypočteme výraz pro finkci z:

(2.9)

 

Stejné řešení ale dostaneme bez námahy, použijeme-li Matlab:

clear

syms ('t','x','y','z','tau','k','S');

S=dsolve('Dx=-1/tau*x','Dy=1/tau*(-y+k*x)','Dz=1/tau*(-z+k*y)','x(0)=1','y(0)=0','z(0)=0');

z=S.z;

pretty(z)

 

 

 

                                        t    2  2

                                 exp(- ---) t  k

                                       tau

                             1/2 ----------------

                                          2

                                       tau  

 

Pro hodnoty k=20,  =10 např. dostaneme:

 

z=S.z;

z=subs(z,{'k','tau'},{'20','10'});

pretty(z);

ezplot(z,[0,100]);grid;

  

 

 

                                               2

                              2 exp(- 1/10 t) t

 

 

 

Pro hodnoty k=20,  =20  pak dostaneme: 

z=S.z;

z=subs(z,{'k','tau'},{'20','20'});

pretty(z);

ezplot(z,[0,100]);grid  

 

 

                                                2

                             1/2 exp(- 1/20 t) t

 

 

Zatím jsme uvažovali tři propojené kompartmenty. Je možno odvodir obecné řešení rovnice 2.9 pro N kompartmentů:

(2.10)

2.3 Odpovědi jednoduchého modelu neuronu

Rovnice známá jako Naka.Rushonova (1966) uvádí vztah mezi intenzitou posdráždění P myšlenou jako čistý postsynaptický potenciál dostihující místo generování impulzů, jejichž frekvenci S(P) je možno vyjádřit:

(2.11)

V této rovnici (vyjadřující obdobný vztaj jako kinetika Michaelis-Mentenové) znamená M maximální frekvenci vzruchů a určuje bod, kdy S(P) dostihuje poloviny svého maxima.

Tak např. hodnoty M a  se pohybují kolem 100 a 50 (M=100, =50), zatímco hodnoty N bývají různé, tak např pro laterální genikulární neurony je N=1.4 pro visuální kortikální neurony N=2.4  a pro neurony střední temporální kůry se hodnoty pohybují kolem N=3.0 (Sclar 1990). Albrecht a Hamilton (1982) udávají střední hodnoty N=3.4 a M=120.

Pro matematické "pohodlí" uvažujme zatím M=100 a N=2.

 

 clear

syms ('M','P','ro','N','S');

     ro=50;M=100;N=2;

S=M*P^N/(ro^N+P^N);

ezplot(S,[0,140])

grid 

 

 

 

 

Při N=4 je křivka strmější:

 

 clear

syms ('M','P','ro','N','S');

     ro=50;M=100;N=4;

S=M*P^N/(ro^N+P^N);

ezplot(S,[0,140])

grid 

 

 

 

Uvažujme odpověď jednoho neuronu na arbitrální stimulus. Dané úrovni stimulu P bude odpovídat určitá frekvence stimulů S(P) [vzruchy/sek] podle Naka-Rushonovi rovnice (2.11) . Ta ovšem bude nabíhat s určitým zpožděním, popsatelným diferenciální rovnicí. Odpověď jednoho neuronu R (vzruchy/sec) na daný stimulus P, kde S(P) je požadovaná frekvence podle rovnice (2.11) Wilson popisuje jako:

(2.12)

podle věty 1 můžeme psát:

(2.13)

Pro jednoduchost řešení uvažujme že daná intenzita stimulu P (a tedy I odpovídající S(P) je celou dobu konstantní. Na počátku je odpověď R nulová – tedy počáteční podmínka: A=R(0)=0;

pak prostou integrací dostaneme:

(2.14)

Kontrola pomocí Matlabu nám ukáže, že jsme rovnici vyřešili dobře:

clear

syms ('M','P','ro','N','S','s','R','r','tau');

 

R=dsolve('DR=1/tau*(-R+S(P))','R(0)=0');

pretty(R)  

 

 

                                          t

                            S(P) - exp(- ---) S(P)

                                         tau  

 

 

Vypočtěme si řešení pro hodnoty P=40, P=80 a P=120 a =15 msec. Funkční hodnotu S(Pp spočítáme podle rovnice (2.11). Předpokládejme že v rovnici (2.11) bude M=100, N=2 a =50

 

clear

syms ('M','P','ro','N','S','s','R','r','tau');

R=dsolve('DR=1/tau*(-R+S(P))','R(0)=0');

ro=50;M=100;N=2;tau=15;

P=40

S=M*P^N/(ro^N+P^N)

r=subs(R,{'S','tau'},{S,'15'})

ezplot(r,[0,100])

grid

 

 

P =

    40

S =

   39.0244

r =

1600/41-1600/41*exp(-1/15*t)

 

 

 

clear

syms ('M','P','ro','N','S','s','R','r','tau');

R=dsolve('DR=1/tau*(-R+S(P))','R(0)=0');

ro=50;M=100;N=2;tau=15;

P=80

S=M*P^N/(ro^N+P^N)

r=subs(R,{'S','tau'},{S,'15'});

ezplot(r,[0,100])

grid  

 

P =

    80

S =

   71.9101

 

 

 

clear

syms ('M','P','ro','N','S','s','R','r','tau');

R=dsolve('DR=1/tau*(-R+S(P))','R(0)=0');

ro=50;M=100;N=2;tau=15;

P=120

S=M*P^N/(ro^N+P^N)

r=subs(R,{'S','tau'},{S,'15'});

ezplot(r,[0,100])

grid 

 

P =

   120

S =

   85.2071

 

 

Pokud P(t) je konstantní, tj. S(P(t)) = S = konst, systém se po čase dostane do rovnovážného stavu kdy dR/dt=0. Jeho hodnotu můžeme z rovnice (2.12) anadno spočítat:

(2.15)

Rovnovážný stav je takový stav, kdy se proměnné systému v čase nemění. Z rovnice (2.12) a (2.14) vyplývá, že systém se asymptoticky blíží rovnovážnému stavu při , hovoříme, že je asymptoticky stabilní.

2.3 Excitační a inhibiční postsynaptické potenciály

 

 

Proud přes iontové kanály je určen Ohmovým zákonem:

 

 

kde I je proud, g je vodivost, V je membránový potenciál a E  je rovnovážný (ekvilibrační) potenciál daný Nerstovou rovnicí:

(2.16)

 

 

Při 20°C  RT/F=25mv. Pokud V=E (s opačným znaménkem), jedná se o reverzní potenciál, při kterém neteče žádný proud. 

 

(2.17)

 

Kde -  je vodivost a ekvilibrační potenciál oblasti pasivího vytékání,  je vodivost a ekvilibrační potenciál v oblasti kanálků pro excitační postsynaptický potenciál (EPSP), a  je vodivost a ekvilibrační potenciál oblasti kanálků pro inhibičního postsynaptického potenciálu (IPSP).

Na obr. 2.6 ve Wilsonovi jsou sledovány odpovědi na presynaptické podráždění v oblasti EPSP při různé úrovni postsynaptického potenciálu  - pokud je postsynaptický potenciál blízky nule, je nulová odpověď – neteče žádný proud a tudíž ekvilibrační potenciál pro EPSP kanálky je nulový.

 

Obdobně (viz obr. 2.6 ve Wilsonovi) je ekvilibrační potenciál pro kanálky v oblasti IPSP kolem -75 mv (GABA reversal).

 

Pro lidské neurony pak platí (při , :

(2.18)

 

Pokud není žádný mediátor jak v excitační tak v inhibiční synapsi (příslušné vodivosti jsou nulové), pak v ekvilibriu je V=-75 mV. Předpokládejme, že je aktivní excitační synapse po dobu 1 ms (příslušná vodivost =2 nS, a vodivost inhibiční synapse zůstává nulová ():

(2.19)

 

To umíme řešit podle věty 1.1 (pří počátečních podmínkách A=-75, a pro):

(2.20)

Lépe je ale přenechat dřinu strojům a nechat si to spočítat symbolickým manipulátorem v Matlabu:

clear

syms ('t','V','V1','V0', 'ge','gi','Eq');

ge=2;

gi=0;

V0=-75;

V=dsolve('DV=-1/12.5*((V+75)+ge*V+gi*(V+75))','V(0)=V0');

V1=eval(V)

ezplot(V1,[0,1])

grid

  

 

V1 =

-25-50*exp(-6/25*t)

  

 

Po jedné milisekundě se kanálky EPSP uzavřou ( ), hodnota potenciálu bude V=-64.3 mV.

pak pro t>1:

(2.21)

Opět podle věty 1 dostáváme (kde za čas t dosadíme (t-1), abychom posunuly počáteční podmínky):

(2.21)

 

Správnost výsledku můžeme zkontrolovat Matlabem:

clear

syms ('t','V','V2');

V=dsolve('DV=-1/12.5*(V+75)','V(1)=-64.3');

V=simple(V)  

 

V =

-75+107/10*exp(-2/25*t+2/25)  

 

Můžeme si to i graficky znázornit – řešíme první rovnici pro ge=2, gi=0, a počáteční podmínky V(0)=V0=75, vykreslíme graf pro první rovnici V1, nastavíme možnost překreslování (hold on), spočítáme V0 pro čas 1 msec, a spočítáme druhou rovnici při V(1)=V0, gi=ge=0 a vykreslíme graf pro funkce V2:

 

 

clear

syms ('t','V','V1','V2','V0', 'ge','gi','Eq');

hold off

ge=2;

gi=0;

V0=-75;

V=dsolve('DV=-1/12.5*((V+75)+ge*V+gi*(V+75))','V(0)=V0');

V1=eval(V)

ezplot(V1,[0,1])

hold on

grid on

 

V0=eval(subs(V,t,'1'))

ge=0;

V=dsolve('DV=-1/12.5*((V+75)+ge*V+gi*(V+75))','V(1)=V0');

 

V2=eval(V)

ezplot(V2,[0,30])

grid on

hold off

 

 

V1 =

-25-50*exp(-6/25*t)

V0 =

  -64.3314

V2 =

-75+96094268539191552/8314692867212933*exp(-2/25*t)

  

 

Rovnice (2.18) při otevřených kanálcích na excitačních I inhibičních synapsích:

(2.21)

 

 

Podle věty 1 dostaneme:

(2.24)

Ověříme v Matlabu:

 

clear

syms ('t','V','V3');

V=dsolve('DV=-1/12.5*(15*V+975)','V(0)=-75')

ezplot(V,[0,1])

grid

 

 

V =

-65-10*exp(-6/5*t)

  

Na rozdíl od situace, kdy jsou otevřenys pouze kanálky pro excitační synapsi, v tomto případě je vzestup napětí jen na -68 mV (místo -64.3 mV).

 


Kapitola 3 Dvoudimenzionální systémy a stavový prostor

 

3.1 Diferenciální rovnice druhého stupně

 

Obrázek 3.1 Vektory sil a vektory rychlosti a zrychlení u hmotného bodu hmotnosti "m" na ideální pružině délky  "x".

 

Z jednoduchého příkladu z mechaniky vyplývá:

tj.

(3.1)

 

Pro řešení této diferenciální rovnice druhého stupně zkusíme dosadit exponenciální funkci :

(3.2)

rovnice bude mít řešení když:

(3.3)

tedy

 

(3.4)

máme pak dva kořeny při nichž po dosazení  nebo  platí:

(3.4 a)

 

(3.4 b)

 

 

obecným řešením je lineární kombinace:

(3.5)

jak se snadno přesvědčíme když dosadíme (3.5)  do (3.1), kdy dostaneme součet rovnic (3.4a) a (3.4b):

 

 

Konkrétní řešení závisí na koeficientech a a b, jejichž hodnoty určíme z počátečních podmínek v čase .

 

Počáteční hodnoty proměnné x a její první derivace v čase t=0 označíme jako  a :

 

 

 

Při bude , takže dosadíme-li do tohoto výrazu (3.5) dostaneme soustavu dvou rovnic:

 

 

a z počátečních hodnot  a  určíme koeficienty a a b.

Diferenciální rovnici druhého stupně (3.1) můžeme přepsat do tvaru soustavy dvou diferenciálních rovnic prvního stupně:

(3.6)

Hovoříme o vyjádření diferenciálních rovnic vyšších stupňů v tzv. normální formě. Tak např. lineární diferenciální rovnici třetího stupně:

(3.7)

můžeme vyjádřit jako soustavu tří obyčejných diferenciálních rovnic prvního stupně:

(3.8)

 

Do normální formy se dají převést lineární diferenciální rovnice, nebo diferenciální rovnice alespoň takové, které jsou lineární v jejich nejvyšší derivaci – pak hovoříme o tzv. kvazilineárních diferenciálních rovnicích. Uvažujme známou van der Polovu rovnici (1926), formalizující srdeční frekvenci:

 

přepíšeme do tvaru soustavy dvou diferenciálních rovnic prvního stupně:

(3.9)

Mírně modifikovaná verze této rovnice (FitzHugh-Nagumova rovnice) je využívána jako zjednodušená aproximace Hodgkin-Huxleyho vztahu. Kvazilineární diferenciální rovnice se běžně převádějí do normální formy vyjádřené soustavou diferenciálních rovnic prvního stupně. Důležité je, že ne vždy se soustava nelineárních diferenciálních rovnic dá převézt z normální formy zpět do jedné rovnice vyššího řádu.

 

Lineární a kvazilineární soustava diferenciálních rovnic se používá často při formalizaci biologických dějů. Tak nař. uvažujme jednoduchý neuronový inhibičně-excitační obvod dvou neuronů:

 

Obrázek 3.2 Soustava excitačního neuronu E a inhibičního neuronu I, excitační neuron stimuluje inhibiční neuron, inhibiční neuron tlumí excitační neuron.

 

Zpětnou vazbu můžeme vyjádřit jako soustavu dvou rovnic:

(3.10)

kde funkce S je Naka-Rushtonova funkce (2.11) závislosti rychlosti vzruchů na intenzitě stimulu. Konstanty a a b reprezentují váhy synapsí z a  a p je intenzita vstupu do excitačního neuronu.

3.2 Řešení diferenciálních rovnic druhého stupně

 

Soustavy diferenciálních rovnic druhého stupně, zapsané ve své normální formě:

 

můžeme vyjádřit v maticovém zápise:

nebo elegantněji (oboustranná šipka nad symbolem označuje čtvercovou matici, jednostranná šipka označuje sloupcový vektor):

(3.11)

Při zkoumání soustavy rovnic nejprve hledáme hodnoty ustáleného stavu tj. takové hodnoty proměnných x a y při kterých je systém v rovnováze, tj. kdy se hodnoty proměnných x a  y nemění v čase. Jinak řečeno, jejich hodnoty jsou konstantní a tudíž jejich derivace jsou nulové.

 

Pokud tento stav existuje, pak v tomto stavu platí že:

a tedy z (3.11) plyne, že:

(3.11a)

Chceme-li tedy určit hodnoty ustáleného stavu, řešíme soustavu lineárních algebraických rovnic. V maticovém zápise pak toto řešení bude vypadat:

(3.12)

kde  je inverzní matice

 

Konkrétní příklad:

(3.13)

Její řešení je v Matlabu snadné :

 

A=[-9,-5;1,-3];

B=[7;1];

Xeq=inv(A)*B

  

 

Xeq =

   -0.5000

   -0.5000  

tedy:

Soustava má ustálený stav při x=0,5 a z=0,5. Pro řešení tohoto příkladu je snadné také sestavit simulinkové schéma. Při numerickém řešení v Simulinku vidíme, že systém se z počátečních hodnot x=0 a y=0 rychle dostane rychle do ustáleného stavu s hodnotami x=0,5 a y=0,5.

Obrázek 3.3 Simulinkové schéma řešení příkladu (3.13)

 Má-li systém jedinečný ustálený stav (jako výše uvedený příklad), můžeme rovnici (3.11) přepsat tak, aby počátek souřadnic vycházel z ustáleného stavu.

(3.13a)

protože  je konstanta, pak derivace se rovná derivaci :

z (3.11) pak vyplývá:

podel (3.13a) dosadíme za :

protože (podle 13.11a):

pak se v nových souřadnicích zbavíme sloupcového vektoru :

(3.14)

nebo v podrobnějším zápisu:

 

Nyní dosadíme do proměnných  a exponenciální funkce:

tj.:

 

první derivace pak budou:

Ze (3.14) pak vyplývá že:

(3.15)

Tedy:

(3.15a)

Protože:

kde I je jednotková matice:

po dosazení do (3.15a):

 

neboli:

(3.16)

Výraz v závorkách  je čtvercová matice (proto ten trik s jednotkovou maticí). Rovnice (3.16) má

triviální řešení při , které je nezajímavé. Bude však jediným řešením, pokud bude existovat inverzní matice  . Můžeme se o tom přesvědčit, když vynásobíme-li obě strany rovnice (3.16) touto inverzní maticí:

 

Výsledkem součinu matice se svou inverzní maticí je jednotková matice, na levé straně rovnice však opět dostaneme nuly:

Tedy jediné řešení je pak pouze:

Netriviální řešení (při ) existuje jen tehdy, pokud inverzní matice neexistuje. Inverzní matice k matici  neexistuje, pokud je její determinant nulový:

(3.17)

Požadavek nulovosti determinantu je požadavkem charakteristickou rovnici. Tato rovnice patří k soustavě diferenciálních rovnic (3.14), a odtud pak určit hodnotu  :

V konkrétním výše uvedeném příkladu (3.13) dostaneme:

(3.18)

Řešením kvadratické rovnice dostaneme hodnoty . Tyto hodnoty jsou vlastní čísla matice .

 

Při dvou kořenech a bude obecným řešením lineární kombinace:

protože podle (3.13a) platí, že, pak:

 

Terminologická poznámka: vlastní čísla matice

 

Problematika vlastních čísel matic souvisí s netriviálním řešením rovnice:

(3.18a)

kde je čtvercová matice velikosti n-n, a  je sloupcový vektor s délkou sloupce n (netriviálnost řešení znamená, že se hledá ještě jiné řešení, než když všechny prvky sloupcového vektoru jsou nuly), n různých hodnot , které jsou řešením výše uvedené rovnice se nazývají vlastní čísla matice A (eigenvalues) a příslušné hodnoty vektoru se nazývají pravé vlastní vektory (right eigenvectors).

Mimochodem, tato rovnice odpovídá rovnici (3.15).

Sloupcový vektor všech hodnot :

V Matlabu tento sloupcový vektor z čtvercové matice můžeme získat velice snadno pomocí funkce eig:

 

lambda=eig(A)

 

 

A=[-9,-5;1,-3]

lambda=eig(A)  

 

A =

    -9    -5

     1    -3

lambda =

    -8

    -4  

 

Pokud vyjádříme hodnoty jako hlavní diagonálu ve čtvercové matici  a sloupcové vektory  , odpovídající jednotlivým hodnotám seskupíme jako sloupce čtvercové matice :

 

pak z rovnice (3.18a) vyplývá, že:

 

Matice a získáme z matice v Matlabu snadno pomocí funkce eig:

 

[X,D]=eig(A)

 

A=[-9,-5;1,-3]

[X,D]=eig(A)

%ověření platnosti A*X=D*X:

AX=A*X

DX=D*X  

 

A =

    -9    -5

     1    -3

X =

   -0.9806    0.7071

    0.1961   -0.7071

D =

    -8     0

     0    -4

AX =

    7.8446   -2.8284

   -1.5689    2.8284

DX =

    7.8446   -5.6569

   -0.7845    2.8284  

 

Pro n-tou hodnotu jako prvek na n-tém řádku a n-tém sloupci matice :

a n-tý sloupec matice , označený jako:

musí platit původní rovnice (13.18a):

V Matlabu je možno se o tom snadno přesvědčit:

A=[-9,-5;1,-3]

[X,D]=eig(A)

X1=X(:,1)

lambda1=D(1,1)

AX1=A*X1

lambda1X1=lambda1*X1

X2=X(:,2)

lambda2=D(2,2)

AX2=A*X2

lambda2X2=lambda2*X2

  

 

A =

    -9    -5

     1    -3

X =

   -0.9806    0.7071

    0.1961   -0.7071

D =

    -8     0

     0    -4

X1 =

   -0.9806

    0.1961

lambda1 =

    -8

AX1 =

    7.8446

   -1.5689

lambda1X1 =

    7.8446

   -1.5689

X2 =

    0.7071

   -0.7071

lambda2 =

    -4

AX2 =

   -2.8284

    2.8284

lambda2X2 =

   -2.8284

    2.8284  

 

Vraťme se nyní k problematice řešení diferenciálních rovnic druhého stupně podle věty 2.

Řešení našeho konkrétního příkladu:

(3.19)

začneme s určením hodnot rovnovážného stavu a vlastních hodnot matice :

 

A=[-9,-5;1,-3];

B=[7;1];

Xeq=-inv(A)*B

lambda=eig(A)

[X,D]=eig(A)  

 

Xeq =

    0.5000

    0.5000

lambda =

    -8

    -4

X =

   -0.9806    0.7071

    0.1961   -0.7071

D =

    -8     0

     0    -4  

 

Dostali jsme:

 

 

Nyní musíme z hodnot počátečního stavu (tj. z hodnot  ) určit hodnoty koeficientů .

Podle (3.15):

pro , ,  pak platí:

 

po roznásobení dostaneme:

 

clear;

syms('a1','b1','t');

[-9,-5;1,-3]*[a1*exp(-8*t);b1*exp(-8*t)]  

 

ans =

[ -9*a1*exp(-8*t)-5*b1*exp(-8*t)]

[    a1*exp(-8*t)-3*b1*exp(-8*t)]  

Dostaneme dvě rovnice:

 

clear;

syms('a1','b1','t','A1');

A1=solve('-9*a1*exp(-8*t)-5*b1*exp(-8*t)=-8*a1*exp(-8*t)','a1')

A1=solve(' a1*exp(-8*t)-3*b1*exp(-8*t)=-8*b1*exp(-8*t)','a1')

 

  

 

A1 =

-5*b1

A1 =

-5*b1  

  

 

Řešením obou rovnic vzhledem k proměnné je vztah:

Vztahu při ,  tedy vyhovuje vektor:

(3.19a)

(hovoří se o vlastním vektoru matice , který odpovídá vlastnímu číslu  matice .

 

Totéž provedeme pro druhé vlastní číslo matice , :

 

 

Po roznásobení dostaneme:

 

 

Z rovnic:

Vyjádříme vztah :

 

clear;

syms('a2','b2','t','A2');

A2=solve('-9*a2*exp(-4*t)-5*b2*exp(-4*t)=-4*a2*exp(-4*t)','a2')  

A2=solve(' a2*exp(-4*t)-3*b2*exp(-4*t)=-4*b2*exp(-4*t)','a2')  

 

A2 =

-b2

A2 =

-b2  

Řešením obou rovnic vzhledem k proměnné je vztah:

Vztahu při ,  vyhovuje vektor:

(3.19b)

 

Dosadíme-li do původní rovnice:

obdržené vztahy mezi koeficienty  a :

pak dostaneme:

(3.19c)

Při lze ze známých počátečních hodnot  (resp. z a ) vypočítat koeficienty a ze vztahu:

 

 

Tak např. při vyřešíme soustavu rovnic:

 

clear;

syms('b1','b2','S');

S=solve('-5*b1-b2+0.5=0','b1+b2+0.5=0','b1','b2');

b1=S.b1

b2=S.b2

  

 

b1 =

.25000000000000000000000000000000

b2 =

-.75000000000000000000000000000000  

 

Dostaneme řešení:

 

 

Upravíme-li (3.19c) a dosadíme hodnoty a  , dostaneme konečný výsledek:

(3.19d)

Wilson ve své knize pro rychlé řešení soustavy diferenciálních rovnic druhého stupně podle věty 2 programy LinearOrder2.m a Equlibrium.m.

 

Použití Wilsonova programu LinearOrder2.m

 

Nejprve musíme znát hodnoty ekvilibria – lépe je používat přímo Matlab, (Wilsonův program Equilibrium.m je natvrdo nastaven na počítání jediného příkladu):

 

A=[-9,-5;1,-3];B=[7;1];

%ekvilibrium je při A*X+B=0

Xeq=-inv(A)*B  

 

Xeq =

    0.5000

    0.5000  

V programu LinearOrder2 řešíme rovnice posunuté do počátku souřadnic, proto si dále musíme přepočítat počáteční podmínky: , tj v našem případě:

Program se nás ptá na hodnoty prvků matice  označené jako vstupy A1, A2, A3, A4 a dále na přepočítané hodnoty počátečních podmínek (označené jako vstupy B1 a B2).

V našem případě tedy po vyvolání LinearOrder2 interaktivně zadáme:

A1=-9

A2=-5

A3=1

A4=-3

B1=-0.5

B2=-0.5

 

Obrázek 3.4 Výstup z Wilsonova programu LienarOrder2.m

 

 

Použití programu LinearOrder2 pro analytické řešení diferenciálních rovnic druhého stupně není však nutné, pro rychlé přímé řešení stačí napsat pár řádek v symbolickém manipulátoru Matlabu:

 

clear;

syms('x','y','t','S');

S=dsolve('Dx=-9*x-5*y+7','Dy=x-3*y+1','x(0)=0','y(0)=0');

x=vpa((S.x),5)

y=vpa((S.y),5)

  

 

x =

-1.2500*exp(-8.*t)+.75000*exp(-4.*t)+.50000

y =

-.75000*exp(-4.*t)+.25000*exp(-8.*t)+.50000  

 

3.3 Negativní zpětná vazba

 

Čípky sítnice (C) stimulují horizonální buňky (H), které mají inhibiční vazbu na čípky (analogie obrázku 3.2). Kinetiku můžeme popsat:¨

(3.20)

 

Konkrétní hodnoty koeficientů , a , počáteční podmínky . Po dosazení dostaneme:

(3.21)

 

V rovnovážném stavu budou derivace nulové (předpokládejme konstantní hodnotu L). Tedy:

 

 

Vidíme že:

Zpětná vazba při konstantním L tlumí odpověď na pětinu, při bude

Takže:

V maticovém vyjádření (při ):

(3.20a)

 

 

Posuneme koordináty do ekvilibria:

 

Do (3.20a) dosadíme:

 

Dostaneme:

 

Po roznásobení:

 

Protože poslední část výrazu jsou podmínky pro ekvilibrium a tedy:

 

Pak se v nových souřadnicích dostaneme:

(3.20b)

Nalezení analytického řešení původní diferenciální rovnice 2. stupně pak spočívá v řešení rovnice (3.20b). Musíme si jen pamatovat, že počáteční podmínky v nových souřadnicích jsou posunuty o hodnoty ekvilibria:

Vlastní řešení pak je možné naléz pomocí Wilsonova programu LinearOrder2 kam zadáme členy matice (A1=-40, A2=-160, A3=12.5, A4=-12.5) a počáteční podmínky v nových souřadnicích (B1=-2, B2=-2). Výsledek zobrazuje obrázek 3.5. Trigonometrické členy, které se objevují v analytickém řešení vyplývají z toho, že vlastní čísla matice  jsou komplexní:

A=[-40,-160;12.5,-12.5];

lambda=eig(A)  

 

lambda =

 -26.2500 +42.5551i

 -26.2500 -42.5551i  

 

 

Podle Eulerova vzorce:

můžeme komplexní exponenty, které dosazujeme do rovnic (3.20b) vyjádřit pomocí sinů a kosinů:

 

 

 

 

Obrázek 3.5 Výstup Wilsonova programu LinearOrder2 při řešení rovnice negativní zpětné vazby v sítnici.

Celý příklad ovšem můžeme velmi snadno vyřešit přímo pomocí symbolického manipulátoru Matlabu:

 

clear;

syms('x','y','t','S');

syms('c','h','S','t','Z')

S=dsolve('Dc=1/0.025*(-c-4*h+10)','Dh=1/0.08*(-h+c)','c(0)=0','h(0)=0');

c=simplify(S.c);

h=simplify(S.h);

%vypočtení a zaokrouhlení konstant

c=vpa((c),5)

h=vpa((h),5)

%tisk grafu

i=0;