�ݺ�ߣ

MATLAB: soluzione di equazioni
differenziali e fitting non lineare
Focus principale sull’utilizzo delle funzioni
ode45 e lsqcurvefit
Michele Scipioni
Ph.D. Student
Dip. Ingegneria dell’informazione
Università di Pisa
Corso di Immagini Biomediche, 02 Dicembre 2016

2
Outline
 PARTE 1 - Soluzione di un sistema di equazioni differenziali in MATLAB
 PARTE 2 - Fitting di modelli non lineari in MATLAB
 PARTE 3 - Esercitazione: implementazione e fitting di modelli bi-compartimentali
A.A. 2016/2017 MATLAB: soluzione di equazioni differenziali e fitting non lineare

3
PARTE 1
Soluzione di un sistema di equazioni differenziali in MATLAB

4
Obiettivo
MODELLO BI-COMPARTIMENTALE
DESCRITTO DA UN SISTEMA DI 2
EQUAZIONI DIFFERENZIALI
1. Prendere confidenza con le funzioni integrate in MATLAB dedicate alla soluzione delle ODE.
2. Impostare le opzioni e i parametri richiesti dal motore di soluzione ODE.
3. Passare il nostro modello ed eventuali parametri aggiuntivi al risolutore.

5A.A. 2016/2017 MATLAB: soluzione di equazioni differenziali e fitting non lineare
Risoluzione numerica di equazioni differenziali
MATLAB fornisce un gran numero di strumenti per risolvere numericamente delle equazioni differenziali.
Noi ci concentraremo sui due metodi principali già integrati in tutte le versioni, ossia le funzioni ode23 e
ode45, entrambe basate sul metodo di risoluzione Runge-Kutta.
Entrambe queste funzioni si presentano con la medesima interfaccia e si aspettano gli stessi parametri
in ingresso.
• y: vettore di soluzione, in cui ogni colonna è associata ad una variabile (più di una se lavoriamo con
un sistema di ODE), ed ogni riga rappresenta un istante temporale riferito al vettore t.
• odefun: funzione esterna (m-file) che restituisce un vettore colonna con I’uscita del sistema ODE
• tspan: specifica il vettore dei tempi di integrazione per risolvere il Sistema
 se tspan è un vettore di due elementi, questi sono trattati come tempi di inizio e fine, e gli
step di integrazione intermedi sono scelti arbitrariamente dal risolutore
 se tspan ha più di 2 valori, i risultati sono calcolati soltanti per gli specifici istanti temporali
richiesti
• y0: vettore delle condizioni iniziali per (t=0)
[t, y] = ode45(odefun, tspan, y0)

Esempio: equazione differenziale del primo ordine
𝑑𝑦
𝑑𝑥
= 𝑥𝑦2
+ 𝑦; 𝑦 0 = 1 𝑥 ∈ [0; 0.5]
1) Definire la funzione 𝐨𝐝𝐞𝐟𝐮𝐧: 𝒇 𝒙, 𝒚 = 𝒚′

𝑑𝑦
𝑑𝑥
= 𝑥𝑦2
+ 𝑦; 𝑦 0 = 1 𝑥 ∈ [0; 0.5]
2) Passare la funzione appena creata in input al risolutore ode45()

𝑑𝑦
𝑑𝑥
= 𝑥𝑦2
+ 𝑦; 𝑦 0 = 1 𝑥 ∈ [0; 0.5]
3) Specificare un tspan personalizzato
Quello che è interessante puntualizzare riguardo il
tspan è che MATLAB in background continua ad
usare più o meno lo stesso timestep predefinito
per risolvere l’equazione differenziale, variando
soltanto il modo in cui restituisce in output il
risultato finale. Lavorare con un tspan
personalizzato quindi non altera in modo eccessivo
l’accuratezza della soluzione numerica.

𝑑𝑦
𝑑𝑥
= 𝑥𝑦2
+ 𝑦; 𝑦 0 = 1 𝑥 ∈ [0; 0.5]
4) Opzioni di personalizzazione
Ad ogni iterazione del risolutore viene calcolato un errore. Se chiamiamo 𝑦 𝑘 l’approssimazione di
𝑦 𝑥 𝑘 allo step k, ed 𝑒 𝑘 l’errore corrispondente, MATLAB adatta il suo passo di integrazione in
modo che risulti:
𝒆 𝒌 ≤ 𝒎𝒂𝒙(𝑹𝒆𝒍𝑻𝒐𝒍 ∗ |𝒚 𝒌|, 𝑨𝒃𝒔𝑻𝒐𝒍),
Dove i valori di default sono 𝑅𝑒𝑙𝑇𝑜𝑙 = 10−3
= .001 and 𝐴𝑏𝑠𝑇𝑜𝑙 = 10−6
= .000001.
N.B.
Con questa convenzione, se la soluzione |𝑦 𝑘| diventa molto grande, di conseguenza lo è anche l’errore, e per
evitare che il risolutore si arresti prima di convergere è necessario ridurre RelTol. Viceversa, se la soluzione
tende ad avere valori molto piccoli (< 10−6
), allora è AbsTol a dover essere ridotto.

𝑑𝑦
𝑑𝑥
= 𝑥𝑦2
+ 𝑦; 𝑦 0 = 1 𝑥 ∈ [0; 0.5]
4) Opzioni di personalizzazione
𝑑𝑦
𝑑𝑥
= 𝑥𝑦2
+ 𝑦; 𝑦 𝑥 =
1
1 − 𝑥

Esempio: sistemi di ODE
𝛼 = 1,2 𝑦1 0 = 2
𝛽 = −0,6 𝑦2 0 = 1
𝜎 = −0,8
𝜌 = 0,3
Modello Predator-Prey
Vettore colonna
212
2
211
1
yyy
dt
dy
yyy
dt
dy





Esempio: passaggio di parametri aggiuntivi
Immaginiamo di voler continuare a lavorare all’esempio del modello predatore-preda, ma variando il
valore delle costanti del modello 𝛼, 𝛽, σ e 𝜌, ma ovviamente senza modificare manualmente il file
MATLAB.

13
PARTE 2
Fitting di modelli non lineari in MATLAB

Obiettivo
STIMA AI MINIMI QUADRATI
1. Prendere confidenza con le funzioni integrate in MATLAB dedicate al fitting non lineare.
2. Impostare le opzioni e i parametri richiesti dalla funzione lsqcurvefit.
3. Passare il nostro modello ed eventuali parametri aggiuntivi alla funzione di fitting.

Least Squares Estimator
L’obiettivo di uno stimatore ai minimi quadrati è trovare un vettore x che sia un minimizzatore locale della
funzione data dalla somma degli scarti quadratici tra curva misurata e uscita del modello teorico,
possibilmente sottoposta ad alcuni vincoli:
𝑚𝑖𝑛 𝑥 𝐹(𝑥) 2
2
= 𝑚𝑖𝑛 𝑥
𝑖
𝐹𝑖
2
(𝑥)
tale che 𝐴 · 𝑥 ≤ 𝑏, 𝐴𝑒𝑞 · 𝑥 = 𝑏𝑒𝑞, 𝑙𝑏 ≤ 𝑥 ≤ 𝑢𝑏.

lsqcurvefit
[x,resnorm,residual,exitflag,output] =
lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)INPUT
• 𝒇𝒖𝒏(𝒙, 𝒙𝒅𝒂𝒕𝒂): funzione non lineare che fornisce in output la stima della curva da fittare
 𝒙: vettore dei parametri incogniti
 𝒙𝒅𝒂𝒕𝒂: variabile indipendete (generalmente il vettore temporale a cui è campionato ydata)
• 𝒙𝟎: stima inziale dei coefficienti del modello da fittare
• 𝒚𝒅𝒂𝒕𝒂: curva misurata (deve essere della stessa misura dell’output di fun)
• 𝒍𝒃, 𝒖𝒃: lower e upper bounds tali che 𝑙𝑏 ≤ 𝑥 ≤ 𝑢𝑏 (devono avere stessa dimensione di x0, o essere matrici vuote se
non si è interessati a vincolare la stima parametrica)
• 𝒐𝒑𝒕𝒊𝒐𝒏𝒔: opzioni di ottimizzazione
OUTPUT
• 𝒙: parametri del modello stimati (stessa dimensione di x0 e del vettore di parametri preso in input da fun(x,xdata))
• 𝒓𝒆𝒔𝒏𝒐𝒓𝒎: somma degli scarti quadratici (residui) corrispondenti all’output del fitting
• 𝒓𝒆𝒔𝒊𝒅𝒖𝒂𝒍: vettore dei residui [fun(x,xdata)-ydata] alla soluzione x
• 𝒆𝒙𝒊𝒕𝒇𝒍𝒂𝒈: valore che descrive la condizione di uscita
• 𝒐𝒖𝒕𝒑𝒖𝒕: struttura contenente informazioni riassuntive sul processo di ottimizzazione

Esempio: fitting di una funzione esponenziale
𝑦 = 𝑎1 exp −𝑏1 𝑡 + 𝑎2 exp −𝑏2 𝑡 𝑡 ∈ [0; 2]
1) Definire la funzione 𝐦𝐨𝐝𝐞𝐥𝐥𝐨: 𝒚 = 𝒇𝒖𝒏(𝒙, 𝒙𝒅𝒂𝒕𝒂)

𝑦 = 𝑎1 exp −𝑏1 𝑡 + 𝑎2 exp −𝑏2 𝑡 𝑡 ∈ [0; 2]
2) Usiamo il modello implementato per simulare una misura rumorosa ydata

𝑦 = 𝑎1 exp −𝑏1 𝑡 + 𝑎2 exp −𝑏2 𝑡 𝑡 ∈ [0; 2]
3) Passare la funzione modello e il vettore misurato con input ad lsqcurvefit

𝑦 = 𝑎1 exp −𝑏1 𝑡 + 𝑎2 exp −𝑏2 𝑡 𝑡 ∈ [0; 2]
4) Passare parametri aggiuntivi alla funzione modello

𝑦 = 𝑎1 exp −𝑏1 𝑡 + 𝑎2 exp −𝑏2 𝑡 𝑡 ∈ [0; 2]
5) Principali opzioni per personalizzare il risultato del fitting
x0 = [1 1 1 0]
options = optimset('Display', 'none');
options.TolFun = 1e-6;
options.TolX = 1e-6;
options.MaxFunEvals = 1e6;
options.MaxIter = 1e6;
lb = [0. 5. 5. 0.];
ub = [5. 15. 10. 3.];

22
PARTE 3
Esercitazione: implementazione e fitting di modelli
bi-compartimentali ad un dataset PET
Link per scaricare i dati sperimentali per l’esercitazione

Punti principali
1. Implementare modello bicompartimentale mediante ode45
2. Implementare modello bicompartimentale mediante la soluzione analitica bi-esponenziale
3. Verificare che i due modelli si comportano in modo analogo simulando delle curve usando i medesimi
parametri cinetici
4. Estrarre una input function ed una curva tissutale dal dataset clinico fornito, e stimare i parametri
cinetici usando entrambe le implementazioni del modello compartimentale
5. Usando l’implementazione più veloce delle due, applicare il fitting a tutti i voxel dell’immagine,
generando delle mappe parametriche
___________________________________
FACOLTATIVO: sfruttando l’implementazione dell’algoritmo ml-em sviluppato in una precedente esercitazione,
cimentarsi nell’implementazione del metodo diretto di stima da sinogrammi, descritto nell’ultima parte della
precedente lezione.

1. Implementare modello bicompartimentale
mediante ode45
𝐶1 𝑡 = 0 = 𝐾1
𝐶2 𝑡 = 0 = 0
Consiglio: implementare due funzioni distinte nello stesso m-file
function Ct = modello_ode(k, t)  funzione principale che calcola l’uscita della ode45 e da in
output Ct=C(:,1)+C(:,2);
function Cp = sistema_ode(t,C,param)  funzione interna in cui vengono esplicitate le due equazioni
differenziali del sistema

2. Implementare modello bicompartimentale
mediante la soluzione analitica bi-esponenziale
𝑪 𝑻 = 𝜶 𝟏 𝒆−𝜷 𝟏 𝒕 + 𝜶 𝟐 𝒆−𝜷 𝟐 𝒕
𝛼1 = 𝐾1
𝑘3 + 𝑘4 − 𝛽1
𝛽2 − 𝛽1
𝛼2 = 𝐾1
𝛽2 − 𝑘3 − 𝑘4
𝛽2 − 𝛽1
𝛽1 =
𝑘2 + 𝑘3 + 𝑘4 − 𝑘2 + 𝑘3 + 𝑘4
2 − 4𝑘2 𝑘4
2
𝛽2 =
𝑘2 + 𝑘3 + 𝑘4 + 𝑘2 + 𝑘3 + 𝑘4
2 − 4𝑘2 𝑘4
2
function Ct = modello_exp(k, time)
Questa funzione deve prendere gli stessi
ingressi della precedente (le 4 costanti del
modello compartimentale e il vettore dei tempi)
e svolgere al suo interno la conversione nel set
di variabili ausiliarie prima di calcolare 𝑪 𝑻(t).

3. Verificare che i due modelli si comportano in
modo analogo
Simulare le risposte impulsive (IRF) dei
due tessuti con i seguenti valori di
ingresso:
𝑡 = 0: 35; % 𝑖𝑛 𝑚𝑖𝑛𝑢𝑡𝑖
𝑘 = [0.3 0.5 0.6 0.1];

4. Stimare i parametri cinetici dai dati sperimentali (1/2)
1. Importare il dataset allegato all’esercitazione
2. Estrarre 2 curve tempo-attività dalla matrice 4D “Volume”
• ROI Input function: slice 7, time frame 4
• ROI tessuto: slice 23, time frame 24
3. Importare dalla struttura “PETInfo” il vettore dei tempi:
• time = mean(PETInfo.time,2)./60;
4. Implementare le formule riportate nella diapositiva 33 della lezione precedente
utilizzando la funzione riportata nella prossima slide. Usare quest’ultima come
ingresso all’algoritmo di ottimizzazione lsqcurvefit
5. Impostare le seguenti opzioni di ottimizzazione:
• lb = [ 0. 0. 0. 0. 0. ]; % K1 k2 k3 k4 vB
• ub= [ 1. 1. 1. 1. 1. ]; % K1 k2 k3 k4 vB
• k0= [ 0.1 0.1 0.01 0.01 0.1]; % K1 k2 k3 k4 vB

function Ct = modello_conv_IF(k, time, Cp)
vB = k(5);
dt = 0.1;
t = 0:dt:time(end);
IF = max(0,interp1(time,Cp,t,'linear',0) );
sol = conv(modello_exp(k, t), IF)*dt;
tsol = (0:dt:2*max(t));
% Taking into account natural radioactive decay
dk = log(2)/109.8; % radioactive decay constant for F-18
sol = sol .* exp(-dk * tsol);
% Taking into account blood fraction in tissue
sol = max(0,interp1(tsol,sol,time,'linear',0) );
Ct = (1-vB)*sol + vB* Cp;
Ct(Ct<=0) = eps; K1 k2 k3 k4 vB
0.2820 0.5087 0.6218 0.0804 0.0656
4. Stimare i parametri cinetici dai dati sperimentali (2/2)

5. Generazione di mappe parametriche
2. Estrarre 1 curva tempo-attività dalla matrice 4D “Volume” relativa alla input function (ROI da slice 7, time
frame 4)
4. Scegliere una fetta del volume PET importato (ad esempio la 24 da cui si è estratta la ROI precedente o
qualsiasi altra a piacere) e fittare voxel per voxel il modello, come fatto al punto precedente.
• img = squeeze(Volume(:,:,24,:));

5. Generazione di mappe parametriche
Consiglio
Impostare un “filtro” che riduca il numero di pixel da fittare, per
accelerare l’esecuzione del codice ed anche per ottenere delle mappe
più pulite (evitando di fittare anche pixel di sfondo). Le mappe qui
riportate sono ottenute:
• Moltiplicando il volume importato x1000 (img = img.*1000)
• Fittando solo le curve estratte da img tali che (max(curva)>=15)

6. Stima diretta dei parametri da sinogramma
(FACOLTATIVO)
2. Estrarre 1 curve tempo-attività dalla matrice 4D “Volume” relativa all input function (slice 7, time frame 4)
4. Scegliere una fetta del volume PET importato (ad esempio la 24 da cui si è estratta la ROI precedente o qualsiasi altra
a piacere) e fittare voxel per voxel il modello, come fatto al punto precedente.
5. Richiamare l’esercitazione sull’algoritmo di ricostruzione MLEM (28 Ottobre 2016) e:
• utilizzare la matrice di sistema (A) per convertire il volume PET in un sinogramma dinamico (come fatto con il fantoccio l’altra volta)
• “copiare” il blocco di codice che esegue 1 iterazione di MLEM, facendo l’ipotesi che il sinogramma non abbia artefatti da correggere
6. Implementare il metodo diretto di stima parametrica, alternando una iterazione MLEM ad un fitting dell’intero
volume ricostruito, come svolto al punto 5.

6. Stima diretta dei parametri da sinogramma
(FACOLTATIVO)

�ݺ�ߣ

Kinetic_Modeling_02_12_2016

Recommended

More Related Content

Similar to Kinetic_Modeling_02_12_2016 (20)

Kinetic_Modeling_02_12_2016