Identyfikacja funkcji

STANISŁAW KOWALIK i ZBIGNIEW MITEK (Gliwice)

Identyfikacja funkcji ^przejścia kory nadnercza

(Praca ^przyjęta do druku 26.11.1979)

W pracy przedstawiono zastosowanie metod identyfikacji dla wyznaczenia funkcji ^przejścia kory nadnercza. Problem sprowadzono do wyznaczenia nielinio- wej funkcji regresji. Zamieszczono polepszone parametry modelu kory nadnercza.

Wstęp. Badanie zasad funkcjonowania żywych organizmów będące przedmiotem

badań nauk medycznych, a pod kątem zastosowań technicznych również bioniki, napotyka wiele trudności ze względu na specyfikę budowy i własności układów

biologicznych [I], [2], [3], [6], [10].

Wobec braku ogólnych zasad interpretacji wyników ^badań prowadzonych na

układach biologicznych najczęściej stosowaną metodą ich analizy jest modelowanie.

Podstawową trudnością, która się narzuca, jest problem uzyskania zależności między zmiennymi występującymi w ^układzie. Dla rozwiązania tych trudności,

modelowanie funkcji ^przejścia kory nadnercza oparto na metodzie identyfikacji nieliniowej funkcji regresji [5].

1. Sf~rmułowanie problemu. Dobrze znanym ^układem regulacji hormonalnej jest układ przysadka mózgowa - kora nadnercza. Zjawisko samoczynnej regulacji poziomu jednego z wielu hormonów kory nadnercza, kortikosteronu, przedstawia rysunek 1 [2]. Poziom kortikosteronu regulowany jest w zależności od wartości

zadanej zmieniającej się z zakłóceniami, na które ^narażony jest organizm (stresy).

Produkcją kortikosteronu przez nadnercze steruje hormon adrenokortikotropowy (ACTH), którego wydzielanie przez przysadkę mózgową jest regulowane przez CRF - peptyd produkowany przez podwzgórze. Wydzielanie CRF w ^dużym stopniu ^zależy od poziomu ^krążącego kortikosteronu, co zamyka ^pętlę [2], [9].

Sprawdzenie zasad współdziałania poszczególnych elementów układu biolo- gicznego oraz funkcjonowanie ^układu regulacji hormonalnej jako całości możliwe jest po określeniu parametrów każdego elementu układu regulacji i określeniu postaci transmitancji. Przy badaniach ilościowych najczęściej stosowano metody analogowe wspomagane ^analizą przebiegów czasowych [2], [3]. Coraz ^częściej stosuje ^się metody numeryczne [I], [6], [10], co znacznie ^ułatwia wyznaczanie para- metrów modelu.

[85]

2. Doświadczenie identyfikacyjne na korze nadnercza. W laboratorium endokry- nologicznym w Laval w Quebec w Kanadzie, dla uchwycenia zależności ilościowych charakteryzujących korę nadnercza, przeprowadzono elementarne doświadczenie

identyfikacyjne, które pod ^względem medycznym ^było jednak bardzo ^złożone [2].

Na samcach szczurów Spraque-Dewley (300 osobników) przeprowadzono ekspe- ryment o następującym przebiegu:

- ^wycięto przysadki mózgowe (hipofizektomia)- odpowiada to przecięciu pętli sprzężenia zwrotnego ^układu biologicznego z rys. 1,

- wyczekano 4 godziny, tj. do momentu, gdy poziom kortikosteronu ^spadł do zera - odpowiada to ustaleniu zerowych warunków początkowych,

- pobudzono ^korę nadnercza zastrzykiem hormonu ACTH w trzech ^różnych dawkach C 0 w mili uncjach: 0,5 mU, 1,0 mU, 2,0 mU.

Wartość zadano

Ośrodkowy CRF

układ

nerwowy

Przysadko ACTH Koro mózgowo nadnercza

kortikosteron

wydzielony kortikosteron plazmatyczny

Obciążenie

Rys. 1 Schemat blokowy ^układu regulacji hormonalnej przysadka mózgowa - kora nadnercza Odpowiedź układu biologicznego, tj. poziom kortikosteronu we krwi, próbko- wano po upływie 2,5, 5, 10, 20, 30, 60, 120 minut ^metodą Silbera zmodyfikowaną

przez Guillemina. Wyniki ^badań na szczurach po uśrednieniu przedstawia rysunek 2.

Zastrzyk hormonu ACTH można traktować jako impuls Diraca [2], wobec

powyższego zatem krzywe doświadczalne z rysunku 2 można interpretować jako

µg g

70 60 50 I.O 30 20 10

o ^{10 20 30 I.O 50 min}

Rys. 2. Charakterystyki impulsowe „doświadczalne" kory

nadnercza dla trzech ^różnych dawek ACTH 2

charakterystyki (odpowiedzi) impulsowe g(t) [4] kory nadnercza, gdzie

(1) g(t) = 2- ¹ {G(8)},

G(s) jest transmitancją operatorową, a 2- 1 { • } oznacza odwrotne przekształcenie

1 Laplace' a.

Dane doświadczalne opracowano ilościowo [2].

Przedstawienie wyników jest o tyle interesujące, że pozwoli sprawdzić skutecz-

ność stosowanych metod analogowych.

Transmitancję kory nadnercza na podstawie kształtu uzyskanych charakterystyk impulsowych (rys. 2) przyjęto w postaci [2]:

(2) _G( 8 ) = (1 +8T ₁ )(1 _+8T2 ^K )(I _+8T3 i.

Wartości stałych czasowych i wzmocnienie dla trzech różnych pobudzeń C 0

wyznaczono na maszynie analogowej. Wyniki ujmuje tabela 1.

K T1 T2 T3

a2

TABELA 1

Zestawienie wyników uzyskanych metodami analogowymi

Co= 0,5 mU C 0 = 1,0 mU

180 560

1 min 10 s 1 min 10 s

1 min 35 s 2 min

2 min 30 s 2 min 30 s

6,445 127,354

Co= 2,0 mU 675 1 min 10 s 4 min 4 min 30 s

445,680

Miarą błędu modeli w tabeli 1 jest <J ² (6). Mimo że transmitancja (2) jest bardzo prosta, to jednak wyznaczone ^stałe czasowe ^dały charakterystyki impulsowe ^obcią­

żone zna<;znym ^błędem (tabela 1). Wobec powyższego poszukiwanie modelu w po- staci bardziej ^złożonej metodami analogowymi wydaje ^{się być} bezcelowe. Ten fakt skłania do stosowania metod numerycznych.

3. Identyfikacja charakterystyki impulsowej. Transmitancja G(s) ^określona wzorem (2) ma bieguny jednokrotne i można ją napisać w postaci [4]:

(3) G(8) = _!!_.!___ + _!?.l:_ + ~.

8-81 8-82 8-83

Charakterystyka impulsowa g(t) = 2- ¹ {G(s)} jest równa g(t) = b1es1t+b2es2t+b3es3t,

bieguny transmitancji są jednokrotne, rzeczywiste i ujemne (T _{1 ,} T _{2 ,} T ₃ > O) (rys.

2), tak więc:

(4) g(t) = b 1 e-tfT1+b2 e-tfT2+b3 e-tfTs,

gdzie T 1 , T 2 , T 3 to stałe czasowe transmitancji (2).

Ogólnie, dla rozpatrywanej klasy modeli identyfikacja charakterystyki impul- sowej odpowiada wyznaczeniu nieliniowej funkcji regresji „najlepszej" w sensie najmniejszej sumy kwadratów o postaci [5]:

(5) g(t) = bo+bif1(t)+bif2(t)+ ··· + ^bmfm(t),

gdzie.fi(!)= e-t/T!· W przyjętym modelu matematycznym w postaci nieliniowej funkcji regresji (5) założono częściową znajomość charakterystyki obiektu, pole-

gającą na znajomości funkcji fi(t) (i= I, ... , m). Za miarę błędu modelu ^przyjęto

średnią wartość z kwadratów ^błędów

N

(6) a2 = -k L ^(gi-gi)2,

i=l

gdzie gi = g(t;) oznaczają wyjścia obiektu, ag; = g(ti) wyjścia modelu.

Przyjęcie a ² dla porównania uzyskanych modeli uznano za wystarczające.

Strukturę modelu ^określa równanie (2), wyznaczenie modelu sprowadza ^się do identyfikacji parametrycznej nieliniowego równania (4), tzn. do optymalnego doboru ^stałych czasowych Ti (i = 1 , 2, 3) oraz parametrów b 1 , b 2 , b 3 , minimali-

zujących błąd średniokwadratowy modelu.

W ^związku z trudnościami stosowania metod analitycznych dobór ^stałych cza- sowych prowadzono na podstawie metody typu Monte Carlo [7].

W ten sposób zadanie numeryczne zastępuje się zadaniem z rachunku prawdo-

podobieństwa, które z kolei rozwiązuje się na drodze eksperymentu statystycznego.

Eksperyment statystyczny prowadzono w następujący sposób:

(a) wybrano arbitralnie wartości początkowe TiO ^stałych czasowych;

(b) wykorzystując analizę regresji utworzono tzw. macierz ^wejść uogólnionych [5]

o postaci

(7) [

fo(t1) f1(t1) ··· fm(t1) - X=, ~0(~2). ~1(~2~ :„. ^~m(~2).

fo(tN) !1 (tN) ^„. fm(tN)

gdzie/ 0 (t) = 1,f;(t) = e-t/Ti (i= 1, .. „ m) (uwag a: na początku jest Ti =

= Tio (i = 1, ^„ „ m));

(c) wybór nowych trzech kolejnych liczb pseudolosowych L1i (i= 1, 2, 3) z ^ciągu liczb pseudolosowych o rozkładzie normalnym N(O, 1 ),

(d) wychodząc z wartości początkowych T; 0 obliczono wartości stałych czasowych z zależności

(8) ^~ = ~ 0 ^+L'.ł 1 (i= I, 2, 3);

(e) jeśli kolejny obliczony model był „lepszy" od modelu wyznaczonego ^wcześniej (mniejsza wariancja resztkowa), to obliczone stałe czasowe T; (i = 1, 2, 3) wykorzystywano dla dalszych ^obliczeń jako nowe TiO (i= 1, 2, 3);

(f) powrót do punktu (b).

Obliczenia prowadzi ^się do momentu uzyskania zadowalających wyników - minimalnego ^błędu modelu przy dopuszczalnym czasie ^obliczeń.

Przedstawiony algorytm ma charakter heurystyczny (tzn. nie ma dowodu jego

zbieżności do ^wartości optymalnych ^stałych czasowych), jednak na podstawie przerachowanego w pracy ^przykładu numerycznego stwierdzono dobrą zbieżność

algorytmu. Czas ^obliczeń dla ^każdego modelu (tabela 2) nie przekraczał IO minut.

Uzyskane tymi metodami wyniki należy traktować ostrożnie. Interpretacja fizy- kalna rezultatów ma sens o tyle, o ile część zmiennych można przyjąć ze źródeł

pozastatystycznych [5] (np. bezpośredni pomiar stałej czasowej zanikania hormonu ACTH we krwi, [2], [8]).

4. Wyniki identyfikacji funkcji przejścia kory nadnercza. Zgodnie z sugestiami uczynionymi w § 3, jedną ze stałych czasowych, mianowicie zanikania hormonu ACTH we krwi, przyjęto z danych literaturowych [2]

T 1 = 1,66666 min (1 minuta 10 sekund).

Pozostałe stałe czasowe T 2 , T 3 dobrano wykorzystując metodę Monte Carlo. Wyniki ujmuje tabela 2.

K T1 T~ Ta

(]2

TABELA 2

Zestawienie wyników uzyskanych ^metodą typu Monte Carlo

Co= 0,5 mU C 0 = 1,0 mU

141,5 580

1,166666 1,166666

2,156 4,042

2,5335 5,374

5,1508 10,7618

Co= 2,0 mU 1628

1,166666 8,12 41,1538 8,8

Rezultaty z tabeli 2 po porównaniu z wynikami z tabeli 1 pozwalają stwierdzić, że uzyskano ponad dziesięciokrotne polepszenie dokładności modelu dla pobu-

dzeń C 0 , wynoszących 1 mU i 2 mU.

Należy zwrócić uwagę, że przeprowadzone obliczenia miały na celu głównie

pokazanie skuteczności proponowanej metody. ^Błąd dla ^pobudzeń C 0 , wynoszą­

cych 1 mU i 2 mU, jest nadal znaczny. Ze względu na fakt, że wykorzystywane dane doświadczalne opracowano w latach pięćdziesiątych [2], a przy współczesnym postępie nauk medycznych dane te zdążyły się zdeaktualizować, „udokładnienie"

modelu wydaje ^{się być} bezcelowe.

Po zamodelowaniu transmitancji z parametrami z tabeli 2 uzyskano charakte- rystyki impulsowe g(t) jak na rysunku 3.

Z powyższych wyników można wyciągnąć następujące wnioski:

(a) Odpowiedź gruczołu nie jest wprost proporcjonalna do amplitudy dawki ACTH

pobudzającej gruczoł. Punkty maksimum z rysunku 3 nie występują jednocze-

śnie - kora nadnercza nie jest więc elementem liniowym.

PROGR/IH PL()T X-T --W~KRES NR 1 WSPOLCZYNNIJ( SKALI -X ~ 0,70E OO • T= O.SOE OO

, •• "''0!'1•1111~

''• ... 'lłlłłl'llld

" • , Il 11 ! I

I I

I I 11 ~;

"•,'''"'1!"11!':

I} • • ~ t: t f łl I' li~ I t

"• /'lfl•1111_10' f'

t .'. (ł

{111!111 \I

~:

1 • .„5 IJl/lł!lll~

i'• 41)1JUIJUłl„

~: ~~:::::'.:;::;-

,.. .1 ,'t'łll"''''' f'l,1·'>11•1!'"''' r. /11on110••'

" .

(~'''·''-

•1111

h: -~~~~::::. ;; ^::~~

"•''f\01,l("!'l•'·

l',1''illl.'f'l1••'·

r.1r.•11111-1:11r

(" t , fi~)

0 !

jl

I•

¹

f·. i 1 ou1·1111r f•.'!1">

¹¹

••f

¹

llf

~: ~ ~~:::,;:;:;

('·, '~łll!l•tt11F (\ 1:11"\"l'lillf

~:~~:~~::::;:;f

f.,

^~11111•

^llttf;

".,"i ~11111111[

(l, 1/.łJłJltllOf ,, • , 1,~111•1111!"

n., ^7P

•11 · ~11·,

r., 7'.' .,,,,",,::

r., g111•1111••F.

r., ~! ~ •11· 1111'

n. ir>o1111011r

r.. 1C..'~'ll'll(ą;

(1, /nllłl!JIHT[

r, /l"~Plllłflf f', 11 l.'l'l llOJ:

r .11 '.'>~.·;~11ur

n. l2tlO~JUll[

(), t?"Jtt1

¹

01.

¹

f:

r • ,'~OłJtlllOE ,.. ,,!>t't't.1111:

,...: 1411t'(:l11Jf r, 14:>1111t1ll[

(', l,ł"1 1Jlll.'lltJr.

f" ,l5)'...łllllO[

(~, /ldlO(łllOl I"!, l6:>1Jr.111t1E (l, /.7UOfl!1U(

0, l.7~0

¹

1Pl1t:

"· 1.g111.11~t•11r r ,1.r-.:>11t11"•r.

r, /~łJllflllH[

~. ;~;,!;:::;:~~

'': .so-.o•.•o•_if.

11 ,.S1 i1t•lltJ11E

r, ³¹ ~Ultlll!f:

Il • .Słl!UIJllO[

g:~~e:w::::r

~~:.~L~~~:~::::f

r •. ~4!>op11or

o .S5ll011(11)f

„: .S5!>011f!Uf:

I' ,.S60Ul;011r

!" •. ~

¹

,) 'Jl.'llUF.

~I.~ 700(10\IE r· • .s1:>01111or

~: ~~~~!~::~~F

('I, .s9011011ur·

(', ,SO)tJ(!fl\!f:

r, 1.oot1t111or

~ =~~~:;:_:g:;~

r , 41) 1.•l1\.

¹

•1ł~

~: ~~g:~:::;::ę

<'• '•3t•t11·110E n ,4~)t•f•Pqf

<'•'•4000llOE

I\•

'''·~''11lłłJ~

r,, 45llUPllłlf.

r ,1,5)UOllOL

<'• '•"llOt'l1t1r n 1 4f,)Vł)•111r

(\, l1 "flll1(1

!11.'~

r.

^1,

^7)t1t

¹

^{1111 [}

f'I

1 /, -.~O

li!' 1.

¹

Ut:

~:: ~;~~~!:~::~~r.

r·, '~n'

¹

.1,'l!I, I.

r, ~/:11(14• 1111~

;.co~ oc• ~,lor. P1 ~.1H Q;! 2,7.H'07 ~,2!F. o? z.3SE Ol. ~.4~~ ó? ~,49Ę Ul z.56E o.:! 2,f,JE Ol +

· ~ 1 ,,,. r•o 11 1 ,. nr .. 11 0111• 1q1 ,„. ,, ··•>O'l'IO" ''"''"''fl•l'lPflOf'OOll n1_,o 11noouonnno 11nuo11·101)1100011ooou1Juoonuounnooono;JOO '):)o 'ł0•

¹

0ł~OVJ

pq i•11 ;• J 1 I I I I I I I I

::~ :: ·' , .1F :l ^{l l : : l l l l}

::·~ :; ł ¹ : 3 F 2 I I I I ^ł l

:'. ~ :~ I l , , ³ ! ~ l ⁷ I z l I l I l

::i :: l : 1 I l ³ ! ² l ² ! ! I !

::~ ::---;---;---;-' ---; ---;---;---5--;-:'3·---;---;---;---1

::ł :: ^ł :, , ^ł I I zF ³ 1 I I l

:: ^~ l ^{1. , :} I I I l I ^:i 1 I l

n ^{I : I I ~ I l I ~ f 3 l}

1•1 11 I~ 1 I I 2 I I I I I 3 I

~'.~ :: _______ , _1 _________ ! _________ .! __ ---2- -~ !---____ r _________ 1 _________ 1 _________ r _________ ^{1 _____ ~ .s-- ~}

:::: (• 1, I ^{~ l} I I I I ! I ³ .il

;:~ , : I l I I -1 I I il

C'I. 1 1 12 I I I I I I I :li

:~~ 1 1 I ;>: l l : l l l l ^{-~ J} l

::~ :;, 1 I ^{;> 2 :} I I I l I I l :. ³ l

::~ ~:~ --- ---~z---;---;---;---- ---;--- ----~---;--- ---r---;-3-J---l

::~ ~· ¹ /I : l I I I I I ³ 1 l

Ili. 1 ~ . I I I I I I I I I I

::~ ł / I : I l l l ! ! l l

P,'. 1 l I I I I I I I 3 I I I

::~ 1 „ ² I I ^ł I I I I .s ³ I I I

;; ^{~ ~} -2 ^~ ---;---;---;---; ---;---;--- --:;-i---;---;---:

::~ 1 ^~ l : l l l : l I l l

I'<'. 12 I I I I I _\ I I I I

•·<'. ~ ~ I : I : : ^{3 : : : : I}

::~ , ;: I ^{I I I I I I I I I}

::5. ~ I l l l .s I~ I l l I I

:: ^~ r---;---;---;---;--- ^~ -- ^~- -;---;---;---;---;---:

llł 7. I I I I J I I I I I I

::~ ~ I I i I I l ł I I

"'- 2 ł : I ~ ^j 1 I I l I l I

:;~ ~ I I I -~ - ! I I I I I I

::~ ~ l : - P I l l I I l l

:: ~ ~ ---;---;---33- ~- ---;---;---;---;---;--- --i---:

Il<'. 2 I : 3 l I I l l l I l

::~ ~ l ^{I ~:'I I I I I I I I I}

:;~ ~ l ~ ³ „ l l I l l l l I

::~ ~ l ^{33 :} l l I l I l l l

::~ ~---i-- - ^{-; 3} ~--7---f---f--- --f---;---;---f---- ---t---l

::~ ;?. I 3 I I I I t I 1- 05 mu I I I

::~ ~ łl / \ ¹ l \ \ l ^\2 ~ WmU l \ !

11<'. 2 3 I I I I 13 2 ^{l I '}

tit'. 7. J3 ^I l ^{I I I I ~ .OmU I I I}

::~ ~ ~ ^{3 : 1} l l : : l : j

:~~ ~---11 _________ ! _________ ! _________ ! ___ . --· -__ 1 ____ -____ ! _________ l_ ---___ r __ --_____ 1 _________ J

t•t'. 2 j I I I I I I I I I I

l'ł 2 3 l ^{I I I I I I I I} I

11.: z ^~ t I l \ l \ l I l ^l

::~ ~ J I I I I I I I I l

::~ ~ ~ l l l l : : l l ^I

::~ ~ ::~ ~---~---;--- 1 I ----7---·; ---___ : l l i ___ ---;--- l : ^~ ---;---;---;--- l I I --1 I

::~ ~ ~ I : l l \ : l l l i

::~ ~ f. ^{I :} l ^{I : : : :} l l

U ~ 3 ł 1 I I I I I I I

t•.: l 3 t l I I I I I I I

:·~ ~ ~ l : ; l l : l l l J

f• L ~ :; ~; ;; ;.~ ;: ;; ;;; ;; ; ; : ; ; ; ; ; ;;_;;; ~;;;;.;;;;:;;;;;;;;;:; ~;; ~;;;;:;;;;;;;;; ~;; ;; ~;; ;~ ;; ; ~ ;; ; ; ; ; ! :~ =~:; ~; ~· =:

~.l)('f 00 ~. 7nF 01 ~.HE 02 ~.21Ę O~ o.281: 02. ~.35E Ol Q,4ZE 02 ~.~9E oz ~.55f ^{02 0.63E} oz Rys. 3 Charakterystyki impulsowe obliczone na podstawie wyznaczonych parametrów modeli

(b) stałą czasową T 1 odpowiedzialną za zanikanie hormonu ACTH we krwi ^należy

przyjąć z nowszych opracowań poświęconych problemowi.

(c) Gdy przyjmie się stałą czasową T 1 zgodnie z punktem (b), wówczas czynnik

(9) K

I ^ł

równania (2) wyraża zmianę w czasie ^stężenia kortikosteronu nadnercza. Jest on wynikiem ^działania dwu czynników: produkcji hormonu przez ^gruczoł i jego wy- dzielania [2].

Prace cytowane

[1] I. E. W. Be nek a n, I. A. B l o m, Trend prediction as a basis for optima/ therapy, TH-report 18-E-86 ISBN 90-6144-089-6, Edinhoven, August 1978.

[2] I. C. G i 11 e, I. G. Pa q u et, Modelowanie ^układów regulacji hormonalnej na maszynie analogowej-funkcja ^przejść kory nadnercza, Prace Instytutu Automatyki PAN, z. 13, War- szawa 1964.

(3] F. S. Gr od i n s, Control theory and biologica! systems, Columbia University Press, New York 1963.

[4] T. Kaczorek, Teoria ^układów regulacji automatycznej, WNT, Warszawa 1974.

(5) K. M a ^ń cz a k, Technika planowania eksperymentu, WNT, Warszawa 1976.

[6] Z. Mitek, Modele matematyczne układów "biologicznych, praca magisterska, Gliwice 1979.

[7] T. H. Ney 1 or, Modelowanie cyfrowe układów ekonomicznych, PWN, Warszawa 1975.

(8) C. D. Turner, J. T. Bagna r a, Endokrynologia ogólna, PWRiL, Warszawa 1978.

(9) F. E. Y at e s, J. Ur q u h ar t, Control of plasma concentrations of adrenocortical hor- mones, Pysio!. Rev. (1962), 359-443.

[10] III Krajowe Sympozjum TIP- „Analiza i modelowanie procesów fizjologicznych", ^Wrocław

13.IX-14.IX.1974.