Modelli Matematici a Compartimenti

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15 MODELLI MATEMATICI A COMPARTIMENTI Definizione di compartimento compartimento

Il concetto di compartimento è essenzialmente di tipo intuitivo ed operativo ed il suo sviluppo costituisce un tentativo di semplificazione e schematizzazione di sistemi biologici anche molto complessi. È difficile dare una definizione rigorosa di compartimento; qui ci si limita a darne una abbastanza generica da comprendere tutti i casi possibili: il compartimento è una quantità di sostanza che si comporta in un modo caratteristico ed omogeneo dal punto di vista dei fenomeni di trasporto e/o delle trasformazioni chimiche; per omogeneità si intende la impossibilità di distinguere campioni diversi prelevati nello stesso istante all'interno di uno stesso compartimento. Non necessariamente il compartimento coincide con una regione di spazio fisiologicamente

identificabile, contornato da membrane o altre barriere; infatti il compartimento può anche essere una determinata specie chimica, sia considerata a se stante, sia considerata in una certa regione di spazio che essa occupa. Appare chiaro come la definizione di compartimento sia legata ad un approccio macroscopico; macroscopico; in essa non si tiene conto di variazioni locali della quantità di sostanza e/o di specie chimica, funzioni delle coordinate spaziali. I modelli a compartimenti sono quindi a parametri concentrati. concentrati.

Secondo la Teoria dei Sistemi si potrebbe dire che un compartimento dà origine ad una sola variabile di stato, che non dà alcuna informazione sulla distribuzione della sostanza stessa nell'interno del volume considerato. Modelli a compartimenti

Un modello a compartimenti è costituito da un numero finito di compartimenti interagenti attraverso lo scambio/trasformazione di materiale di cui si tiene conto con relazioni intercompartimentali. Di solito i compartimenti vengono rappresentati rappresentati con blocchi collegati da frecce. È da sottolineare che la schematizzazione del sistema allo studio con un modello a compartimenti corrisponde ad una impostazione che congloba tutti i fenomeni fisico-chimici di ciascun trasferimento e/o di ciascuna trasformazione in una relazione intercompartimentale. Si tratta di un approccio molto intuitivo, ma semplificato; quindi la utilizzazione dei risultati forniti dal modello a compartimenti, adottato per la descrizione del sistema biologico, va fatta con cautela. Nella trattazione che segue, oltre a quelle insite nella definizione di compartimento, vengono fatte le seguenti ipotesi semplificative: - costanza del volume di ciascun compartimento, in cui è distribuita la sostanza in esame; - linearità: la quantità di sostanza che, nell'unità di tempo, esce da un compartimento è proporzionale alla concentrazione della sostanza nel medesimo compartimento; - tempo invarianza (o stazionarietà): oltre ai volumi, anche i coefficienti di proporzionalità nelle relazioni intercompartimentali sono costanti nel tempo. Esempio di modello con un solo compartimento

Si consideri una sostanza che, una volta immessa nel plasma, viene fissata o semplicemente eliminata dal circolo ad un ritmo proporzionale alla sua concentrazione plasmatica. Per studiare come varia nel tempo la concentrazione plasmatica (cioè la cosiddetta cinetica della sostanza) si può fare riferimento al seguente semplice modello: u(t)

x(t)

k

V plasma

Modelli a compartimenti

Pagina 1

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15

in cui si ipotizza che la sostanza si trovi in quantità x(t) nel plasma, distribuita uniformemente nel volume V di plasma; k (detta anche coefficiente di eliminazione) è una costante di proporzionalità positiva fra portata istantanea (quantità nell’unità di tempo) di sostanza eliminata e concentrazione plasmatica della sostanza, u(t) è la portata istantanea di sostanza entrante nel plasma (per es. per iniezione). Nell’ipotesi di distribuzione uniforme della sostanza nel volume V di plasma, la concentrazione y(t) avrà lo stesso valore in ogni punto del compartimento plasmatico e sarà pari a y(t) = x(t)/V .

Il principio generale della conservazione della massa è rispettato se è soddisfatta l’equazione del bilancio di massa (massa entrata – massa uscita = massa accumulata ), che va scritta con riferimento al sistema di interesse e ad un determinato intervallo di tempo. Nel caso in esame, con riferimento al compartimento e ad un intervallo di tempo infinitesimo dt (da t a t +dt ), il bilancio di massa può essere scritto nella forma seguente x (t ) u (t ) dt − k dt = dx ( = x& (t ) dt ) . V

Dividendo tutti i termini per dt e riordinando, si ottiene l’equazione differenziale: x ( t ) x& ( t ) = − k + u ( t ) ,

(1)

V

x (t )

dove k

V

è la portata di sostanza eliminata all’istante t e x& (t ) è la quantità di sostanza

accumulata nel compartimento nell’unità di tempo, all’istante t . L’equazione (1) rappresenta la cosiddetta equazione di stato del modello, in forma normale (o canonica). Risolvendo questa equazione differenziale, si può calcolare x(t) (e quindi anche y(t)) per prefissati valori di k , V , x(0) e per ogni andamento temporale di u(t). k

Equazione caratteristica

λ+

Autovalore

λ=−

Costante di tempo

τ=

V

=0

k

V V k Infusione continua di durata infinita

• u(t) = U 1(t) = gradino di ampiezza U

u(t) U 0 0

x(0) ≠ 0 , cioè

t

presenza di sostanza nel plasma prima dell'infusione.

Risulta: x(t) U τ

t t  −  − Quantità di sostanza x(t ) = U τ 1 − e τ  + x(0)e τ  

τ

x(0) 0 0

t

x ( ∞ ) = U τ


Pagina 2


Concentrazione

t t  −  −   U τ 1 − e τ + x(0)e τ     y (t ) =

V

Nel caso particolare in cui x(0) = 0 , cioè assenza di sostanza nel plasma prima dell'infusione, risulta: x(t)

•

τ

U τ t  −  Quantità di sostanza x(t ) = U τ 1 − e τ   

0 0

t

x ( ∞ ) = U τ



Concentrazione

y (t ) =

U τ  V

1 − e

−

t 

τ 

 



Infusione continua di durata finita u(t) U • u(t) = impulso rettangolare di ampiezza U e durata T

0 0

T

t

x(0) = 0 , cioè assenza di sostanza nel plasma prima dell'infusione.

Risulta: Quantità di sostanza x(t)

τ

U τ

0 < t < T

t > T

t  −    x(t ) = U τ 1 − e τ    

x(t ) = x (T )e

−

0 0

T

t

t −T

τ

T  −    x(T ) = U τ 1 − e τ    

Concentrazione y(t ) =


x (t ) V

Pagina 3

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15 Infusione impulsiva (di un bolo di massa M) • u(t) = M δ(t) = impulso di Dirac con area sottesa pari a M

x(0) = 0, cioè assenza di sostanza nel

plasma prima dell'infusione.

L'impulso di Dirac di area M può essere ottenuto dall'impulso rettangolare considerato nel caso precedente ponendo U = M/T al limite per T → 0. Risulta: lim



T 





− M  τ 1 − e τ  = M x (T ) =  . T 

T → 0

Dunque: x(t) M

Quantità di sostanza

x (t ) = Me

−

k V

t

0 0

τ

t

x ( ∞ ) = 0

Concentrazione

y (t ) =

M V

−

e

k V

t

Le unità di misura delle grandezze e dei parametri utilizzati nel modello possono essere le seguenti. - x espresso in unità di massa le portate d'infusione e di eliminazione sono espresse in massa/tempo e M è la massa totale iniettata - x espresso in unità di quantità di materia (moli) le portate d'ingresso e di eliminazione sono espresse in moli/tempo e M è la quantità totale di materia iniettata (in moli) in entrambi i casi, la costante k è espressa in volume (di soluzione)/tempo Osservazione :

il parametro k ha il significato di clearance, cioè volume di soluzione che viene depurata completamente dalla sostanza nell'unità di tempo . Modelli con due compartimenti Esempio n.1 - Reazione chimica reversibile

La cinetica della semplice reazione chimica reversibile fra due sostanze S 1 e S 2 in soluzione S 1

→ ←

S 2

può essere interpretata con un modello a due compartimenti: uno relativo alla sostanza S 1, l'altro alla sostanza S 2. Poichè le due sostanze si trovano nella stessa soluzione di volume V , si dovranno considerare due compartimenti con lo stesso volume di distribuzione V .


Pagina 4


x1 (t)

k 1

x2 (t)

V

k 2

V

soluzione

soluzione

Indicando con x1 e x2 la massa dei soluti S 1 e S 2 e con k 1 e k 2 i coefficienti di scambio fra i due compartimenti, il bilancio di massa per i due compartimenti porta a scrivere le equazioni: x&1 = − k 1 x&2 = k 1

x1

V x1 V

+ k 2

− k 2

x2

V x2 V

(2)

,

(3)

che rappresentano le equazioni di stato per il modello, in forma normale. Poichè non è presente alcuna variabile di ingresso, si dice che il sistema è autonomo. In queste condizioni, x1(t) e x2(t ) possono essere non nulle soltanto se almeno una delle due ha valore iniziale diverso da zero. Si noti che, se le quantità di soluto sono espresse in moli e il volume V è espresso in litri, i due coefficienti k 1 e k 2 assumono il significato di costanti di velocità della reazione chimica reversibile. Risolvendo, si possono calcolare x1(t) e x2(t) per ogni t > 0, per fissati valori di k 1, k 2, V , x1(0) e x2(0). Equazione caratteristica

 1 0 1  − k 1 det λ   − V  k 0 1   1  

Autovalori

λ1 = 0

Costanti di tempo

τ1 = ∞

e e

λ2 = −

τ2 =

k 2  

  = 0 ovvero − k 2  

2

λ +

1 V

( k1 + k 2 )λ = 0

k1 + k 2 V

V k1 + k 2

Risulta: x1 (t) B C A1++C 2

Quantità di sostanza

x1 (t ) = C 1 + C 2 e

−

> 0 C B2>0

t

τ 2

C A1 0 τ

0

x2 (t ) = C 3 + C 4e

−

τ 2

L’andamento è analogo a quello di x (t): parte da C +C e tende asintoticamente a C 1

3

a regime

x1 (∞ ) = C 1 ,

t

2

t

4

3

x2 (∞) = C 3

Le costanti C 1, C 2, C 3 e C 4 si determinano imponendo le condizioni iniziali. Devono infatti soddisfare le condizioni (4) C 1 + C 2 = x1(0) (5) C 3 + C 4 = x2(0).


Pagina 5


Per calcolare le quattro costanti occorrono tuttavia altre due condizioni, che si possono scrivere imponendo anche i valori di x&1 (0) e x&2 (0) . In particolare, si ha λ 2 C 2 = x&1 (0) (6) λ 2 C 4 = x& 2 (0) . (7) Si conclude quindi che dall’equazione (7) si può ricavare C 4, dalla (6) C2, dalla (5) C 3 e dalla (4) C 1. Si può verificare facilmente che all'equilibrio (cioè per t → ∞) risulta x 2 ( ∞ ) k 1 = , x1 ( ∞ ) k 2 che costituisce la costante di equilibrio della reazione chimica esaminata. La costante di tempo τ2 della reazione chimica risulta tanto più piccola quanto più grandi sono i coefficienti di scambio. Questo corrisponde al fatto generale per cui le reazioni chimiche raggiungono tanto più rapidamente la condizione finale di equilibrio, quanto più grandi sono le costanti di velocità. Caso particolare. Se la reazione è unidirezionale, cioè per es. k 2 = 0, le equazioni (2) e (3) si semplificano. In particolare x1(t) diviene indipendente da x2(t), essendo determinata dalla equazione differenziale Osservazione.

x&1 = −k 1

in cui non compare x2. Risolvendo si ottiene:

x1 V

,

x1 (t ) = x1 (0)e

− k 1

x 2 (t ) = x 2 (0) +

t V

k 1 V

t

∫0 x1 (ξ )d ξ . ∞

A regime risulta

x1 ( ∞ ) = 0

e

k 1   − ξ k 1  e V  x 2 ( ∞ ) = x 2 ( 0) +  x1 (0) k 1  V   − V  0 

= x 2 ( 0) + x1 (0) ,

cioè tutta la sostanza S 1 inizialmente presente si trasforma nella sostanza S 2. Esempio n. 2 - Test della Bromosulftaleina (BSF) La BSF è un pigmento usato come indicatore in

prove di funzionalità epatica. Queste prove consistono in una rapida iniezione endovenosa di una piccola quantità di colorante, di solito proporzionale al peso del paziente e nella successiva misura, ad istanti prefissati, della concentrazione plasmatica del colorante; in questo modo si ottiene la cosiddetta curva di scomparsa o di decadimento plasmatico della sostanza. Per l'analisi di questa curva si assume di solito un modello a due compartimenti (Fig. 1) nel quale si tiene conto soltanto degli scambi di colorante fra plasma sanguigno e tessuto epatico e della eliminazione dal circolo mediante fissazione nel fegato. Viene trascurato l'assorbimento del colorante da parte dell'intestino, dei reni e dei tessuti diversi da quello epatico.


Pagina 6

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15 1

x1 (t) u(t)

Fig. 1

k 1

V1

k 2 plasma

2

x2 (t)

k

V 2

3

fegato

k 1, k 2 sono coefficienti di scambio ( ≥ 0) e k 3 è il coefficiente di eliminazione o di fissazione (≥ 0) x1(t) è la quantità di colorante nel compartimento plasmatico e x2(t) è la quantità di colorante nel

compartimento epatico, V 1 è il volume di plasma e V 2 è il volume di distribuzione del compartimento epatico, u(t) è la portata di colorante iniettata nel plasma. Scrivendo il bilancio di massa del colorante per ciascun compartimento, si ottengono le equazioni di stato del modello: x&1 = − k 1 x& 2 = k 1

x1 V 1

x1 V 1

+ k 2

x2 V 2

+u

x − ( k 2 + k 3 ) 2

,

V 2

che di solito, per semplicità di notazione, sono scritte come segue: x&1 = − m1 x1 + m2 x2 + u x&2 = m1 x1 − ( m2 + m3 ) x2

dove

m1 =

k 1 V 1

, m2 =

k 2 V 2

, m3 =

k 3 V 2

(8) (9)

,

.

A queste equazioni si deve aggiungere la cosiddetta equazione di uscita x y = 1

(10)

V 1

per tenere conto del fatto che si misura la concentrazione del colorante nel plasma. Soluzione delle equazioni differenziali (8) e (9) nel caso di iniezione rapida di un bolo di BSF. Se l'iniezione è sufficientemente rapida, u(t) può essere approssimata con un impulso di Dirac M δ(t) di area pari alla massa M di colorante iniettato. Supponendo nulle le condizioni iniziali (cioè assenza di colorante prima dell'iniezione), trasformando secondo Laplace le equazioni (8) e (9) si ottiene il seguente sistema di due equazioni algebriche lineari nelle incognite X 1(s) = [ x1(t)] e X 2(s) = [ x2(t)]: ( s + m 1 ) X1 ( s ) − m 2 X 2 ( s ) = M (11) (12) − m1 X1 ( s ) + ( s + m2 + m3 ) X 2 (s ) = 0 , da cui, risolvendo X1 ( s ) = M

s + m 2 + m3 ∆( s )

,

X 2 ( s ) = M

m1 ∆(s)

,

(13)

dove ∆ ( s ) = ( s + m1 )( s + m2 + m3 ) − m1 m2 = s + ( m1 + m2 + m3 ) s + m1 m3 2

.

(14)

Risolvendo l'equazione caratteristica ∆ ( s ) = 0 si ottengono gli autovalori del modello λ 1, 2 =

dove Modelli a compartimenti

−( m1 + m2 + m3 ) ±

Discr

2 Discr = ( m1 + m2 + m3 ) 2 − 4 m1 m3 = (m1 − m2 − m3 ) 2 + 4 m1 m2

(15) .

(16) Pagina 7


Se i coefficienti di scambio e fissazione sono > 0, dalla (16) appare evidente che Discr > 0. Dunque, gli autovalori risultano sempre reali e distinti . Inoltre, poiché i coefficienti del polinomio caratteristico (14) hanno tutti lo stesso segno, gli autovalori sono sempre negativi. Le equazioni differenziali (8) e (9), con u(t) = M δ(t), ammettono la soluzione seguente per t > 0+: t t x1 (t ) = C 1e λ 1 + C 2eλ 2

(17) λ t λ t (18) x2 (t ) = C 3e 1 + C 4 e 2 dove le costanti C 1 , C 2 , C 3 e C 4 possono essere calcolate antitrasformando le (13), oppure imponendo le condizioni iniziali. Volendo seguire il secondo procedimento, occorre considerare che, subito dopo l'applicazione dell'ingresso impulsivo (istante 0+), risulta x1(0+) = x1(0-) + M = M e x2(0+) = x2(0-) = 0. Tenendo conto delle equazioni (8), (9), (17) e (18), si ottiene: x1 (0+ ) = C 1 + C 2 = M x&1 (0+) = C 1λ 1 + C 2λ 2 = − m1 M ,

da cui

C 1 = M

λ 2 + m1 , λ 2 − λ 1

C 2 = − M

λ 1 + m1 . λ 2 − λ 1

x2 (0+ ) = C 3 + C 4 = 0 x& 2 (0+) = C 3λ 1 + C 4 λ 2 = m1 M

da cui

C 3 = M

− m1

λ 2 − λ 1

,

C 4 = −C 3 .

Per valutare i segni delle costanti C 1, C 2, C 3 e C 4, basta osservare che: λ 2 − λ 1 = − Discr < 0

(m1 − m2 − m3 ) + (m1 − m2 − m3 ) 2 + 4m1 m2 λ 1 + m1 = 2 (m1 − m2 − m3 ) − (m1 − m2 − m3 ) 2 + 4m1 m2 λ 2 + m1 = 2

>

0

<

0 ,

per cui risulta C 1 > 0, C 2 > 0, C 3 > 0, C 4 < 0 . Dunque, x1(t) e x2(t) hanno l'andamento nel tempo del tipo seguente: x1 (t)

x2 (t)

M A

0

0 0

t

È ovvio che y(t) ha un andamento analogo a quello di x1(t).

0

t

Osservazione.

Il polinomio caratteristico (14) può essere determinato, come nell’esempio precedente, senza utilizzare le trasformate di Laplace. Risulta infatti:  1 0 − m1 det λ  − m 0 1    1  Modelli a compartimenti

   = λ 2 + (m1 + m2 + m3 )λ + m1m3 .  − ( m2 + m3 )  m2

Pagina 8

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15 Esempio n. 3 - Test della Bilirubina

La bilirubina è una sostanza presente normalmente nell'organismo: essa deriva dalla distruzione dei globuli rossi come prodotto del catabolismo dell'emoglobina e viene eliminata dal fegato con la bile. Essa viene utilizzata per valutare la funzionalità epatica, con le stesse modalità della BSF. La curva di scomparsa della bilirubina può essere interpretata come risposta di un modello a due compartimenti uguale a quello usato per la BSF. Valgono quindi tutti i risultati ottenuti in precedenza, con la sola accortezza di considerare x1(t) e x2(t) come le variazioni rispetto ai corrisponmdenti valori di equilibrio stazionario prima dell'iniezione. Esempio n.4 - Cinetica del glucosio

Il glucosio attraversa con difficoltà la membrana cellulare per diffusione semplice, poichè è idrosolubile e molto poco liposolubile. In presenza di insulina, ormone prodotto dal pancreas in alcune formazioni cellulari specializzate dette isole di Langherans , il passaggio del glucosio attraverso la membrana cellulare aumenta decisamente. Una descrizione a compartimenti deve quindi prevedere un compartimento del glucosio e uno dell'insulina, distinto dal primo, anche se entrambe le sostanze si trovano nel plasma. Sia il glucosio che l'insulina possono essere eliminati dal plasma, per es. per via renale, e quindi occorrerà prevedere delle vie di eliminazione, così come potranno essere iniettate con portate q1(t) e q2(t), rispettivamente. q1(t) terrà conto anche della introduzione di glucosio con l'alimentazione. Infine, occorre tenere conto dei meccanismi che regolano la produzione di insulina da parte del pancreas e il trasporto del glucosio attraverso la membrana cellulare. Per far questo, si introducono dei termini proporzionali, rispettivamente, alla variazione di glucosio nel sangue e alla variazione di insulina nel sangue. È da sottolineare che questi sono termini anomali, che non rappresentano scambi di sostanza fra i due compartimenti, bensì azioni di controllo. Se indichiamo con X 1(t) la quantità di glucosio nel plasma, X 2(t) la quantità di insulina nel plasma, X 10 il valore di X 1(t) a riposo, X 20 il valore di X 2(t) a riposo, x1(t) = X 1(t) - X 10 , x2(t) = X 2(t) - X 20 , il modello può essere rappresentato come segue 1

- m x 21 2

2

x

x

1

m x 12 1

2 q 2

q1 (t)

plasma

plasma

m 10

m 20

Le equazioni relative al bilancio per i due compartimenti sono x&1 = −m10 x1 − m21 x2 + q1

(35)

,

(36)

x& 2 = m12 x1 − m20 x2 + q2

dove m10 =

k 10 V

, m20 =

k 20 V

, m12 =

k 12 V

, m21 =

k 21 V

sono costanti > 0.

A queste si deve aggiungere l'equazione di uscita, relativa alla variabile misurata. Per esempio, se si misura la concentrazione plasmatica di glucosio, si ha: y (t ) =


x1 (t ) V

.

Pagina 9


Analizzando l'equazione caratteristica, si vede che gli autovalori possono essere complessi coniugati, anche se sempre a parte reale negativa. Infatti: Equazione caratteristica  1 0 − m10 − m21    = 0 ovvero λ2 + ( m10 + m20 ) λ + m10 m2 0 + m12m21 = 0 − det λ      0 1  m12 − m 20   Autovalori λ 1, 2 =

− ( m10 + m20 ) ±

Discr

dove

2

Discr = ( m10 + m20 )2 − 4( m10 m20 + m12 m21 ) .

(37)

Il discriminante può essere sia > 0 che < 0. Discr > 0 Le risposte ad un ingresso impulsivo sono non oscillatorie. Ad es., per q2 = 0 e q1 = Q δ(t), cioè per una iniezione impulsiva di glucosio, le risposte a partire da condizioni di riposo sono le seguenti x 1

x2

Q

0

0 0

0

t

t

Per q1 = 0 e q2 = Q δ(t), cioè per una iniezione impulsiva di insulina, le risposte a partire da condizioni di riposo diventano x2

x1 0

Q

0 0

t

0

t

Discr < 0

Le risposte ad un ingresso impulsivo sono analoghe alle precedenti ma oscillatorie. Per q1 = Q δ(t) e q2 = 0, x1

x2

Q

0

0 t

0

0

t

Per q1 = 0 e q2 = Q δ(t), x

2

x1 0

0 0


Q

t

0

t

Pagina 10

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15 IDENTIFICAZIONE DEI PARAMETRI

I modelli a compartimenti vengono impiegati in campo clinico soprattutto allo scopo di stimare i valori dei coefficienti di scambio e di eliminazione, nonchè dei volumi di distribuzione. Infatti, questi parametri possono essere molto utili nella formulazione della diagnosi, ma, di solito, è difficile (a volte impossibile) misurarli in vivo. Per determinare (identificare) i suddetti parametri, generalmente si effettua un esperimento ingrssouscita consistente nella infusione di un indicatore (cioè una sostanza con caratteristiche fisiche o chimiche che rendono facile la determinazione della sua concentrazione) e nella misura della concentrazione dell'indicatore in qualche compartimento dell’organismo in un certo numero di istanti successivi all'infusione. Identificazione dei parametri del modello con un solo compartimento

Si consideri il modello monocompartimentale visto in precedenza, nel caso di ingresso impulsivo u(t) = M δ(t). L'uscita da misurare sia la concentrazione nel compartimento. Si è ricavato che la risposta all'impulso del modello è espressa come y (t ) =

M V

−

e

k t V

.

Il problema da risolvere è: determinare i parametri del modello, k e V , note la massa M infusa rapidamente all’istante 0 e la concentrazione y (t i) misurata in un numero finito N di istanti t i > 0, con i = 1, 2, .... N . Un possibile procedimento è illustrato di seguito. 1. Si riportano i valori misurati y (t i) su di un piano in cui t è la variabile in ascissa, ln[ y ] è la variabile in ordinata. In altre parole, si riportano i valori misurati in scala semilogaritmica. In questo modo, poichè la concentrazione y(t) predetta dal modello, nella stessa scala logaritmica, è S

S

S

M

ln[ y(t )] = ln

V

−

k V

(19)

t

i punti misurati potranno essere approssimati da una retta con pendenza (0, ln

M V

−

k V

e intercette sugli assi

V M

) e ( ln k

V

, 0).

Dai valori delle due intercette, o di una intercetta e della pendenza, si possono ricavare k e V . Se N = 2 , si potranno facilmente determinare k e V in modo che la retta (19) passi esattamente per i due punti misurati. Tuttavia, i valori di k e V così determinati saranno comunque affetti da errori a causa degli inevitabili errori che si commettono nella misura dei due valori di ys. In pratica, si preferisce effettuare un numero di misure maggiore del numero dei parametri incogniti, cioè, in questo caso, N > 2. Di conseguenza è molto poco probabile che tutti i punti misurati appartengano alla stessa retta, sempre a causa degli inevitabili errori di misura, ma anche di modello (cioè della non perfetta corrispondenza fra modello e sistema biologico). Un possibile approccio è allora quello di determinare la retta che meglio approssima i punti misurati nel senso dei minimi quadrati . L'equazione (19) può essere scritta per gli N istanti di misura, ottenendo N equazioni algebriche nelle due incognite k e V : 2.

M

ln

V

−

k V

t i ≅ ln[ y S (t i )] , con i = 1, 2, .... N .

Le suddette equazioni non sono lineari rispetto a k e V , ma sono facilmente trasformabili in M k e w2 = − . Con questo accorgimento è equazioni lineari rispetto alle nuove incognite w1 = ln V

possibile applicare il metodo dei minimi quadrati. Risulta: Modelli a compartimenti

V

Pagina 11


1 t 1  .  ; A=  .    1 t N 

 w1  ; w  2

x= 

ln[ y S (t 1 )]  .  ; b=  .    ln[ y S (t N )]

x = (AT A)-1AT b.

Una volta determinate w1 e w2 , si possono calcolare V =

M e

w1

(20)

e k = − w2V .

Nel caso di soli due punti misurati ( N = 2) le (20) si riducono alle seguenti w1 =

t 2 ln[ y S ( t 1 )] − t 1 ln[ y S (t 2 )]

w2 =

t 2 − t 1

ln[ y S (t 2 )] − ln[ y S (t 1 )] t 2 − t 1

,

che corrispondono ad una risposta all'impulso del modello che passa esattamente per i punti (t 1 , y S (t 1 ) ) e (t 2 , y S (t 2 ) ) . In modo analogo, ci si può porre il problema di determinare i parametri k e V nel caso di un diverso andamento temporale della variabile d'ingresso u(t). Per esempio, si può considerare la risposta ad un ingresso a gradino di ampiezza U , cioè alla infusione continua di una sostanza. Come si è visto, partendo da condizione iniziale nulla, la concentrazione y(t) è data dall’espressione 

y (t ) =

U τ 

1− e V  

−

t 

τ 

 , dove τ = V/k è la costante di tempo. 

In questo caso, anche passando ai logaritmi, non è possibile ricondursi all’applicazione del metodo dei minimi quadrati. Occorre perciò utilizzare un algoritmo diverso per trovare i valori di V e τ che rendono minima la somma degli scarti al quadrato fra i valori misurati di concentrazione e i valori calcolati come uscita del modello. Osservazione.

Se si considera soltanto la risposta a regime, si può imporre che y S (∞) = y ( ∞) =

U τ V

.

È evidente che, noto U , dalla misura di y (∞) si può ottenere soltanto il rapporto τ / V = 1/k . In altre parole, non è possibile determinare separatamente τ e V. È possibile determinare solo il coefficiente di eliminazione k , mentre non è possibile determinare il volume V . In conclusione, non è possibile identificare due parametri incogniti utilizzando una sola equazione: sarà dunque necessario misurare il valore di y almeno in un istante finito t oltre che a regime. S

S

i

Identificazione dei parametri del modello per la BSF/Bilirubina

Analoghe considerazioni possono essere ripetute nel caso di cinetica di una sostanza, come la BSF/Bilirubina, descrivibile con un modello a due compartimenti. L'obiettivo è di determinare i cinque parametri del modello, che sono k 1, k 2, k 3, V 1 e V 2 , a partire da un esperimento ingresso-uscita. L'esperimento generalmente adottato è la prova ad impulso. Anche in questo caso, l'andamento nel tempo della concentrazione nel plasma effettivamente misurata non segue esattamente l'espressione analitica ricavata per la risposta all'impulso del modello a due compartimenti y (t ) =

1 V 1

(C 1eλ

1t

+ C 2e

λ 2t

) = C 1 Aeλ

1t

+ C 2e

λ 2t

.

(21)

Infatti, il modello matematico è solo un'approssimazione della cinetica del colorante ed inoltre i


Pagina 12


valori misurati sono inevitabilmente affetti da errori. Si tratta allora di approssimare nel modo migliore l'insieme dei valori misurati (curva sperimentale) con la relazione (21). Per comodità di trattazione, possiamo pensare di pervenire all'identificazione dei cinque parametri del modello in due passi: 1) si determinano i valori dei coefficienti C 1 , λ 1 , C 2 , λ 2 per cui l'uscita del modello si adatta meglio alla curva sperimentale ( best-fitting); 2) si calcolano i coefficienti di scambio e i volumi di distribuzione che è possibile determinare utilizzando le espressioni che legano C 1 , λ 1 , C 2 , λ 2 a k 1, k 2, k 3, V 1 e V 2 . Passo 1)

Questo problema può essere affrontato usando un metodo approssimato detto del " peeling", oppure uno dei tanti algoritmi numerici di ottimizzazione. Metodo del peeling

Questo metodo si può applicare se i due autovalori sono molto diversi, cioè se ha: per t molto piccolo

y (t ) ≅ C 1 + C 2e λ 2t

per t molto grande

y (t ) ≅ C 1eλ 1

t

,

,

λ 2 << λ 1 < 0 .

Si

ovvero

ln[ y (t ) − C 1 ] ≅ ln C 2 + λ 2t

(22)

ovvero

ln[ y(t )] ≅ ln C 1 + λ 1t

(23)

Le equazioni (22) e (23) suggeriscono la seguente procedura per la determinazione approssimata dei coefficienti ottimi: 1. si traccia, in scala semilogaritmica, la curva di scomparsa del colorante, quindi si individua l'asintoto, per t → ∞ della curva; 2. il coefficiente angolare dell'asintoto fornisce il valore di λ1 e la sua intercetta sull'asse delle ordinate il valore di C 1 ; 3. si sottrae alla curva di scomparsa il valore C 1 calcolato al punto 2 e si traccia, in scala semilogaritmica, la curva risultante; si individua l'asintoto per t → 0 ; 4. il coefficiente angolare di questo asintoto fornisce il valore di λ2 e la sua intercetta sull'asse delle ordinate il valore di C 2 . Esempio

Nella tabella sono forniti i valori misurati su un paziente per il test della BSF. t (min) y (mg/L)

4 96,0

8 75,0

12 69,8

16 61,0

45 41,8

75 37,5

120 31,3

La procedura descritta porta a scrivere: λ1 ≅

ln y7 − ln y6 t7 − t 6

C 1 ≅ y7 e λ2 ≅


− λ 1t 7

≅

≅ −0, 004

min −1

51 mg/L

ln( y2 − A ) − ln( y1 − A ) t 2 − t 1

≅ −0, 155

min −1

Pagina 13

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15 − t C 2 ≅ ( y1 − A )e λ 2 1 ≅ 84

mg/L

Metodo numerico

Lo stesso problema può essere affrontato in modo più generale riconducendolo al problema di trovare il minimo di una funzione criterio (o costo) che misura la "distanza" fra curva sperimentale e curva fornita dal modello. Di solito, analogamente al caso dei minimi quadrati, si adotta come funzione criterio la somma dei quadrati degli scarti fra la risposta del modello y (t ) = C 1eλ 1t + C 2eλ 2t (funzione dei parametri) e la risposta sperimentale yS (t) in un numero finito N di istanti t i, ovvero E (C 1 , λ 1 , C 2 , λ 2 ) =

N

∑1 [ y

S

(t i ) − y (t i )]2 .

(24)

i=

Dunque, si tratta di determinare i quattro parametri C 1 , λ 1 , C 2 , λ 2 in modo da rendere minima E . Essi dovranno soddisfare le condizioni necessarie per avere un estremo di E ∂ E ∂ E ∂ E = 0; = 0; ∂ C 1 ∂λ 1 ∂ C 2

=

0;

∂ E = 0. ∂λ 2

(25)

Poichè y(t) è una funzione non lineare dei parametri, è arduo risolvere il sistema (25) procedendo per via analitica. Per questo, in pratica si ricorre ad algoritmi di aggiustamento dei parametri (o di minimizzazione, o di o ttimizzazione), che sono generalmente disponibili nelle più diffuse librerie di programmi matematici, incluso MATLAB. La determinazione dei parametri avviene mediante un processo iterativo secondo lo schema illustrato in Fig. 2. Sistema

y (t ) S i

biologico Calcolo della E funzione criterio

u(t)

E

è minima ?

SI

STOP

E

Modello

NO

y(t ) i

Parametri modificati

Algoritmo di aggiustamento dei parametri

Fig. 2 Passo 2)

Noti C 1 , λ 1 , C 2 , λ 2 relativi al minimo di E , per trovare i parametri del modello occorre considerare le espressioni ricavate in precedenza, che costituiscono il seguente sistema di equazioni algebriche non lineari, nelle incognite V 1, m1, m2, m3 . C 1 =

λ 2 + m1 M λ 2 − λ 1 V 1


(26)

Pagina 14


C 2 = −

λ 1 + m1 M λ 2 − λ 1 V 1

(27)

 λ 1 + λ 2 = −(m1 + m2 o  λ 2 = λ 2 (m1 , m2 , m3 ) λ 1λ 2 = m1 m3 λ 1

= λ 1 ( m1 , m2 , m3 )

+ m3 )

(28)

Dividendo la prima per la seconda equazione si elimina V 1 e si ricava m1. λ + m1 =− 2 λ 1 + m1 C 2 C 1

; m1 = −

C 1λ 1 + C 2λ 2

.

C 1 + C 2

Sommando la prima e la seconda si elimina m1 e si ricava V 1 . C 1 + C 2 =

M V 1

;

V 1 =

M C 1 + C 2

.

Dalla quarta equazione si ottiene

m3 =

Dalla terza equazione si ottiene m2

λ 1λ 2

m1

.

= −(λ 1 + λ 2 + m1 + m3 ) .

Dunque, si può concludere che con l'esperimento ingresso-uscita adottato, cioè misurando la curva di scomparsa del colorante nel plasma conseguente ad una infusione impulsiva, è possibile determinare V 1 e i rapporti k 1 / V1 , k 2 / V 2, k 3 / V2 , ma non è possibile determinare V 2. Si supponga di effettuare un esperimento diverso, consistente nel misurare la curva di concentrazione del colorante nel tessuto epatico in risposta ad una infusione impulsiva nel plasma. In questo caso l'uscita del modello a due compartimenti risulta così espressa: y (t ) =

C 3 V 2

(eλ

1t

−e

λ 2t

) = C 3 (eλ

1t

−e

λ 2t

) .

Procedendo in due passi come in precedenza, ci si riduce a dover risolvere il sistema di tre equazioni algebriche non lineari seguente: C 3 =

− m1 M

λ 2 − λ 1 V 2

(29)

 λ 1 + λ 2 = −(m1 + m2 o  λ 2 = λ 2 (m1 , m2 , m3 ) λ 1λ 2 = m1 m3 λ 1

= λ 1 ( m1 , m2 , m3 )

+ m3 )

.

(30)

È evidente che questa volta è possibile calcolare soltanto tre parametri in funzione degli altri. Ad esempio, è possibile ricavare V 2, m2 ed m3 in funzione di m1 : dalla (29)

V 2 =

dalle (30)

m3 =

− m1 M

λ 2

− λ 1 C

λ 1λ 2

m1

, e

m2 = −(λ 1 + λ 2 + m1 + m3 ) .

Si può dunque concludere che con questo esperimento non è possibile determinare univocamente né il volume V 1, nè il volume V 2, né i tre coefficienti m1, m2 ed m3. Per riuscire a determinare tutti i parametri del modello è necessario ricorrere ad eventuali altre tecniche in grado di fornire direttamente il valore di alcuni dei parametri incogniti, oppure ad altri esperimenti ingresso-uscita. Per esempio, nel caso della BSF/Bilirubina, si potrebbero misurare Modelli a compartimenti

Pagina 15


contemporaneamente la concentrazione plasmatica ed epatica del colorante a seguito di una infusione impulsiva. Così facendo, infatti, prima si potrebbero determinare V 1, m1(e quindi k 1), m2 ed m3 dalla curva di scomparsa plasmatica, poi V 2, k 2 e k 3, dalla curva di concentrazione epatica. Occorre tuttavia sottolineare che la misura di concentrazione epatica non viene effettuata in pratica sia per il disagio procurato al paziente, sia per la scarsa affidabilità delle misure. Identificabilità a priori Definizione:

per identificabilità a priori di un modello matematico (anche non a compartimenti) si intende la possibilità teorica di determinare univocamente tutti i parametri incogniti del modello per mezzo di un esperimento ingresso-uscita. Alla luce degli esempi presentati, è evidente che in generale sarà opportuno cercare di verificare l'identificabilità del modello già nella fase di progetto degli esperimenti, prima di effettuare concretamente gli esperimenti stessi. Nel seguito mostreremo che è possibile procedere alla verifica di identificabilità di un modello con un prefissato esperimento ingresso-uscita, senza necessariamente risolvere le equazioni differenziali. Questo è particolarmente utile nel caso di ordine elevato del modello o di esperimenti che richiedano forme d'onda d'ingresso diverse da quelle canoniche. Nel caso di sistemi lineari e tempo-invarianti si può considerare la funzione di trasferimento, che, come è noto, descrive in modo completo il legame ingresso-uscita. Per meglio chiarire il procedimento da seguire, si faccia riferimento al modello della BSF di Fig. 1 e relative equazioni. La funzione di trasferimento risulta espressa da: G(s ) =

X1 ( s ) MV1

=

1

s + m2 + m 3

V 1 s 2 + ( m1 + m2 + m3 )s + m1m3

=a

s+b

.

s 2 + cs + d

(31)

La (31) rappresenta il legame ingresso-uscita per qualunque andamento del segnale di ingresso ed è univocamente determinata dai valori dei coefficienti a, b, c, d. Come in precedenza, si può procedere in due passi per determinare i parametri del modello. Nel primo passo si determinano i coefficienti a, b, c, d che corrispondono al segnale di ingresso applicato e alla risposta misurata. Nel secondo passo si calcolano i parametri del modello tenendo conto delle relazioni che li legano ad a, b, c, d , ovvero: 1 V 1

=a

m2 + m3 = b

(32) m1 + m2 + m3 = c m1 m3 = d

Il sistema (32) è costituito da 4 equazioni e ammette una soluzione univoca per le 4 incognite V 1, m1, m2 e m3 . Dunque, si può dire che, rispetto a questi parametri, il modello è identificabile con l'esperimento ingresso-uscita considerato. Ovviamente, anche il sistema (32) conferma il fatto che non è possibile identificare il volume V 2. Considerando la concentrazione nel fegato come uscita si ottiene: X ( s ) 1 m1 1 G(s) = 2 = =a 2 . 2 MV2 V 2 s + ( m1 + m2 + m3 ) s + m1 m3 s + bs + c Modelli a compartimenti

(33) Pagina 16


La (33) è univocamente determinata dai valori dei coefficienti a, b, c, che sono legati ai parametri del modello come segue. m1 V 2

=a

m1 + m2 + m3 = b

(34)

m1 m3 = c

In questo caso si hanno solo tre equazioni, che non consentono di determinare univocamente nemmeno le 4 incognite V 2, m1, m2 e m3. Con questo esperimento ingresso-uscita il modello non è identificabile. In generale, quando un modello non risulta identificabile, si possono tentare tre strade: - scegliere un esperimento ingresso-uscita diverso; - modificare il modello in modo che risulti identificabile con lo stesso esperimento ingressouscita; - cercare relazioni, non dipendenti dal modello adottato, fra alcuni dei parametri da determinare.


Pagina 17

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15 Esempio

Verificare l'identificabilità a priori del modello a compartimenti illustrato in figura, considerando l'ingresso u(t) indicato e come uscita la concentrazione nel compartimento n. 1. 1

x1 u(t

V 1

k 1

k y=x /V 1

2

x2 2

2

1

Le equazioni differenziali del modello sono: x1

x&1 = −k 2 x&2 = − k 1

che, ponendo m1 =

k 2 V 1

, m2

V 1 x2 V 2

k 1

=

V 2

+ k 1

x2 V 2

+u

,

, possono essere scritte come segue:

x&1 = − m1 x1 + m2 x2 + u x&2 = −m2 x2 .

Trasformando secondo Laplace, con condizioni iniziali nulle, si ottiene ( s + m1 ) X1 (s ) − m2 X2 (s ) = U (s ) ( s + m 2 )X 2 ( s ) = 0 , da cui, risolvendo X1 ( s ) = U ( s )

1 s + m1

,

X2 ( s ) = 0

.

La funzione di trasferimento relativa all'esperimento ingresso-uscita considerato risulta quindi G ( s ) =

1

1

V1 s + m1

.

Dunque, si può concludere che è possibile determinare univocamente soltanto V 1 e m1. Gli altri parametri V 2 e m2 non sono identificabili a priori.


Pagina 18

G. Gnudi - Dispense di BIOINGEGNERIA – sede di Cesena – 2014/15 Esempio

Si consideri un farmaco, la cui cinetica è descrivibile con il modello a due compartimenti illustrato in figura, q(t)

m 3

x1 (t)

m 1

x2 (t)

V 1

m 2

V 2

plasma

tessuti

dove x1 e x2 indicano la massa di farmaco contenuta, rispettivamente, nel plasma e nei tessuti, V 1 e V 2 sono i volumi del plasma e dei tessuti, q(t) è la portata in massa di infusione nel plasma del farmaco stesso, i coefficienti m sono definiti al solito modo. In particolare, m3 tiene conto della escrezione renale. Le equazioni differenziali che descrivono la cinetica del farmaco sono: x&1 = −(m1 + m3 ) x1 + m2 x2 + q;

x& 2 = m1 x1 − m2 x2

Si supponga di misurare la concentrazione del farmaco nel plasma: y(t) = x1 /V 1. La funzione di trasferimento risulta: G(s ) =

s + m2

1

V 1 s 2 + ( m1 + m2 + m3 )s + m2m3

.

Essa è caratterizzata da una costante moltiplicativa e tre coefficienti; dunque, è possibile determinare univocamente V1 , m1, m2 e m3 . Si supponga di misurare la concentrazione del farmaco nei tessuti: y(t) = x2 /V 2. La funzione di trasferimento risulta: G( s ) =

m1

1

V2 s + ( m1 + m2 + m3 ) s + m2m3 2

,

dunque, poichè la funzione di trasferimento è caratterizzata da una costante moltiplicativa e da due coefficienti, non è possibile determinare tutti e quattro i parametri V2 , m1, m2 e m3 ; è possibile soltanto calcolare tre di loro in funzione del quarto.


Pagina 19

Modelli Matematici a Compartimenti

Recommend Documents