20/11/2014

1

SYSTEMS BIOLOGY

BIOTECNOLOGIE MOLECOLARI E BIOINFORMATICA

Metodi di stima dei parametri

Marco S. Nobile 21 novembre 2014

Overview • Ruolo e rilevanza dei parametri di un modello (reprise)

• Derivata e gradiente (reprise)

• Metodi di ricerca locale – Metodo del gradiente, metodo di Levenberg-Marquardt

• Metodi di ricerca globale – Simulated annealing (e tabu search)

– Particle Swarm Optimization

– Algoritmi genetici

– Evolution Strategy

20/11/2014

2

Parametrizzazione del sistema • Per modellare un sistema biochimico (con approccio bottom-up)

non basta la formalizzazione delle reazioni

• Necessaria una opportuna parametrizzazione – Stato iniziale del sistema 𝐗0

– Costanti cinetiche 𝐤 (o stocastiche 𝐜, nel caso di modellazione stocastica)

• Con la scelta corretta dei parametri possiamo produrre simulazioni predittive – Ad esempio, utilizzando Runge-Kutta o SSA

Parametri di un sistema • Determinano il comportamento dinamico del sistema

• Esempio: il sistema biochimico Brusselator – Caratterizzato da dinamica oscillante

– Quattro specie: X, Y, A, B

– Quattro reazioni: r1 : A → X

r2 : B + X → Y

r3 : 2X + Y → 3X

r4 : X → λ

– Parametri: stato iniziale delle specie, costanti cinetiche delle reazioni

A

B

X

Y

20/11/2014

3

Importanza dei parametri

c1 = 1 c2 = 5 x 10−3

c3 = 2.5 x 10−5

c4 = 1.5

c1 = 1 c2 = 5 x 10−3

c3 = 2.5 x 10−4

c4 = 1.5

c1 = 1 c2 = 5 x 10−3

c3 = 2.5 x 10−5

c4 = 1.5 x 10−1

Quantità iniziali: A = 200 molecole (costante), X = 200 molecole, B = 600 molecole (costante), Y = 300 molecole

Costanti di reazione e fenomeni emergenti • Secondo esempio: modello di Schlögl

– Bistabilità determinata da scelta dei parametri stocastici

Con queste costanti, entrambi gli stati hanno la stessa probabilità

Nobile et al., PLoS ONE, 2014

20/11/2014

4

Come determinare i parametri? • La generazione di dati sperimentali adatti per la parametrizzazione,

calibrazione e validazione del modello è – Costosa in termini di tempo e soldi

– Non sempre possibile con le tecnologie attualmente disponibili • Problemi di precisione, sistematicità, ampiezza

• Necessità di un numero esauriente di dati sperimentali quantitativi come punto di riferimento per la modellazione, la simulazione e l’analisi

– Controllo della qualità delle misure

– Automazione

– Standardizzazione

– Condivisione

La misurazione dei dati • Requisiti ottimali per la misurazione dei dati:

– Precisione

– Ampiezza

– Sistematicità

20/11/2014

5

La misurazione dei dati - requisiti [1/2] • Ampiezza:

– A livello dei target numero di geni e proteine nel sistema di interesse, oltre a geni e proteine correlati

– Su scala temporale misure dei target a intervalli di tempo costanti per catturare la dinamica del sistema

• non solo prima/dopo evento, non solo punti campione

– Su scala funzionale/spaziale livelli di trascrizione, interazione fra proteine, modificazioni post-traduzionali, localizzazione

La misurazione dei dati - requisiti [2/2] • Sistematicità:

– Perturbazione sistematica

– Misure simultanee di diversi fattori per lo stesso campione biologico, nelle stesse condizioni sperimentali

– Integrazione di dati consistenti fra di loro e misurati a livelli diversi (mRNA, proteine, ecc.)

20/11/2014

6

La misurazione dei dati – un esempio [1/2] • Costruzione di una rete di regolazione genica:

– Precisione:

• Misura del sistema in vivo

• Scelta dei parametri (binding constants, transcription/translation rates, degradation rate, diffusion rate, ecc.) forte dipendenza dalle caratteristiche del sistema, che può mostrare risposte diverse a diverse perturbazioni (concetto di robustezza)

– Sistematicità:

• Misura a livello di espressione genica e a livello di proteine sullo stesso identico campione (es. interazione proteina-proteina con yeast 2-hybrid method)

La misurazione dei dati – un esempio [2/2] • Ampiezza a livello di target:

– Dati di espressione wild-type + delezione/sovraespressione di tutti i geni (costoso…)

– Misurazione massiccia di tutti i geni, o mirata per i geni più rilevanti (misure “casuali” introducono ambiguità nella ricostruzione della rete di regolazione)

• Ampiezza su scala temporale:

– Misura della variazione (sincronizzata) dei livelli di espressione genica a intervalli di tempo ridotti (fine-grain time-series) conoscenza più precisa del sistema, creazione di ipotesi più corrette, eliminazione di interazioni non consistenti

– Validazione dei modelli e dei parametri tramite confronto con i dati sperimentali

20/11/2014

7

Sottocampionamento Experimental data

Simulation

t

?

?

conce

ntr

ati

on

Scarsità di informazioni

Exp. data species 1

Simulation

t

Exp. data species 2

?

20/11/2014

8

Dati sperimentali: serie temporali

Dati sperimentali: serie temporali Serie temporali delle concentrazioni (o quantità) relative a diverse ripetizioni di uno stesso esperimento perturbando gli stati iniziali

Domanda: quali parametrizzazioni consentono di riprodurre queste dinamiche?

Ogni esperimento viene ripetuto più volte in ogni condizione iniziale

20/11/2014

9

Stimare i parametri (cinetici) • Per ogni reazione 𝑅𝜇 (𝜇 = 1,… ,𝑀)

– Modellazione deterministica: cerchiamo la costante cinetica 𝑘𝜇 associata

– Modellazione stocastica: cerchiamo la costante stocastica 𝑐𝜇 associata

• In entrambi i casi, si tratta di una «collezione» di 𝑀 numeri reali – Dunque, la parametrizzazione di un sistema di 𝑀 reazioni corrisponde a un

punto in un semi-spazio 𝑀-dimensionale

– Matematicamente, è rappresentato da un vettore a 𝑀 dimensioni (ad es., 𝐤 = 𝑘1, … , 𝑘𝑀 ∈ ℝ

𝑀)

Come stimare i parametri? • Una proposta di algoritmo:

1. Scegliamo un vettore iniziale di parametri 𝐤

2. Inseriamo 𝐤 nel modello ed effettuiamo (ad es., con Runge-Kutta, SSA) una simulazione del sistema

3. Confrontiamo i dati sperimentali con la simulazione

4. Se le due dinamiche non sono uguali, «miglioriamo» la parametrizzazione 𝐤 e ripartiamo dal passo 2; altrimenti 𝐤 è la parametrizzazione ottima

• Come effettuare i punti 3 e 4?

20/11/2014

10

Confronto delle dinamiche [1/3] • Least squares fitting

– Denotiamo con 𝐘 = (𝐘1, … , 𝐘𝑄) la serie temporale di 𝑄 campioni osservata

sperimentalmente

• I vari campioni sono presi negli istanti 𝑡1, … , 𝑡𝑄

– Denotiamo con 𝐗 = 𝑓(𝐱0, 𝐤) la serie temporale simulata a partire dallo stato 𝐱0 utilizzando il vettore di parametri 𝐤

• 𝐗𝑞 denota lo stato del sistema all’istante 𝑡𝑞

– La parametrizzazione è corretta se le due dinamiche si sovrappongono, ovvero se

𝐘𝑞 − 𝐗𝑞2

𝑄

𝑞=1

= 0

Confronto delle dinamiche [2/3]

|𝐘𝑞 − 𝐗𝑞|

Simulazione Dati sperimentali

20/11/2014

11

Confronto delle dinamiche [3/3]

𝑓𝑖𝑡𝑛𝑒𝑠𝑠 𝐜 = (𝐘1−𝐗1)2 +⋯+ (𝐘10−𝐗10)

2 1

2

3

4 5

6 7

8 9 10

Q=10

La funzione di fitness • La funzione appena descritta può essere utilizzata per quantificare

quanto la dinamica simulata corrisponda alle osservazioni sperimentali – 𝑓: 𝐤 → ℝ+ ovvero associa a un vettore 𝐤 la sua «qualità»

– 𝑓 prende il nome di funzione di fitness

• La stima di parametri può essere ri-formulata come il problema di minimizzare la funzione di fitness – La soluzione ottima è dove 𝑓(𝐤) = 0

20/11/2014

12

Problemi di ottimizzazione • In generale, in un problema di ottimizzazione si cerca la soluzione

migliore (ottima) 𝐨 in uno spazio di possibili soluzioni 𝒮 – In questo specifico caso, si tratta di un problema di minimizzazione

pertanto cerchiamo la soluzione 𝐨 ∈ 𝒮 tale che

𝑓 𝐨 ≤ 𝑓 𝐱 , ∀𝐱 ∈ 𝒮

– Lo spazio delle possibili soluzioni è detto spazio di ricerca

Spazio di ricerca • Immaginiamo di dover stimare una sola costante

– Una soluzione corrisponde a un punto sulla retta dei numeri reali positivi: lo spazio di ricerca è tutta la semi-retta ℝ+

• Come detto all’inizio, nel caso della stima di 𝑀 parametri, una soluzione corrisponde a un punto nell’iperspazio 𝑀-dimensionale e lo spazio di ricerca è ℝ𝑀

• La funzione di fitness associa un vettore di parametri a un valore reale: definisce una superficie chiamata fitness landscape

20/11/2014

13

Esempio di fitness landscape

Soluzione ottima (ad hoc) c0 = 0.0025 c2 = 5.0

METODI DI RICERCA LOCALE Migliorare progressivamente una stima iniziale

20/11/2014

14

Metodi di ricerca locale • In prossimità del minimo di una funzione la derivata si

annulla (teorema di Fermat)

• Se lo spazio di ricerca è un intervallo chiuso e limitato e la funzione è continua allora ha almeno un punto di massimo e uno di minimo (teorema di Weierstraß)

• Risolvendo 𝑓′ 𝑥 = 0 possiamo identificare l’insieme (finito) dei minimi e massimi della funzione

Derivata

Derivata positiva

Derivata negativa

Derivata uguale a 0

20/11/2014

15

Gradiente [1/2] • Nelle funzioni a più variabili, la derivata è generalizzata dalla

derivata direzionale 𝐷𝐮𝑓 𝐱 = limℎ→0 𝑓 𝐱+ℎ𝐮 −𝑓(𝐱)

ℎ calcolata lungo un

vettore 𝐮 unitario

• Se 𝑓 è differenziabile in 𝐱 allora 𝐷𝐮𝑓 𝐱 = 𝛻𝑓 𝐱 ⋅ 𝐮 dove l’operatore 𝛻 è detto nabla e rappresenta il gradiente

Gradiente [2/2] • Il gradiente è un vettore formato dalle derivate parziali della

funzione 𝑓 – Come nel caso della derivata 𝑓 cresce in

direzione del gradiente

– Individua un campo vettoriale: «seguendo» le frecce giungiamo ai massimi

– Seguendole nel verso opposto, raggiungiamo i minimi

20/11/2014

16

Gradient descent • Se scegliamo una soluzione di partenza arbitraria 𝐤0

possiamo calcolare il gradiente e determinare una nuova soluzione 𝐤1 = 𝐤0 −𝛾𝛻𝑓(𝐤0) con fitness migliore (in questo caso, minore)

• Iterando il processo si converge a un minimo dove il gradiente è zero

• Nella variante di Levenberg-Marquardt il fattore di smorzamento 𝛾 viene corretto ad ogni iterazione per adattarsi alla velocità di convergenza

Marquardt, SIAM J Appl Math, 1963

Esempio di esecuzione di Levenberg-Marquardt

Fitness migliore Fitness peggiore Soluzione iniziale

20/11/2014

17

Il problema dei minimi locali

Fitness migliore Fitness peggiore Soluzione iniziale

• Nell’esempio, mediante la discesa del gradiente siamo giunti a

un minimo locale

• Si evince come la qualità della soluzione identificata dipenda dalla

scelta della soluzione iniziale

• Se l’algoritmo non parte in un intorno dell’ottimo globale non possiamo convergere all’ottimo

Limitazioni dei metodi di ricerca locale • I sistemi biochimici sono spesso multi-modali

e fortemente non lineari – Presenza di molti minimi locali (anche infiniti) – Rumore biologico = altri minimi locali

• Metodi basati sulla discesa del gradiente si intrappolano nei minimi locali – Possibile soluzione: ripetere l’algoritmo più volte,

scegliendo sempre una soluzione iniziale differente

• L’algoritmo si intrappola perché la scelta iterativa delle nuove soluzioni è effettuata con metodi deterministici – Cambierebbe qualcosa se, anziché prendere SEMPRE la direzione del gradiente, «ogni

tanto» andassimo altrove?

20/11/2014

18

METODI DI RICERCA GLOBALE Esplorare lo spazio di ricerca

Simulated annealing (SA) • Euristica ispirata alla «tempratura»

dei metalli – Simile al gradient descent, ma l’algoritmo

consente di prendere direzioni che «peggiorano» la funzione di fitness con una certa probabilità

– SA possiede un teorema di convergenza: dato un tempo infinito, l’algoritmo converge alla soluzione ottimale con probabilità prossima a 1

Kirkpatrick et al., Science, 1983

20/11/2014

19

Algoritmo di SA • Il Simulated Annealing funziona iterativamente

1. Si parte da una posizione casuale 𝐤0 e si setta 𝐤𝑏𝑒𝑠𝑡 = 𝐤0

2. Viene generata una nuova posizione 𝐤𝑛𝑒𝑤 a partire da 𝐤0 (ad es., mediante una perturbazione dei valori)

3. Se 𝑓 𝐤𝑛𝑒𝑤 < 𝑓 𝐤0 allora 𝐤𝑛𝑒𝑤 prende il posto di 𝐤0

4. Se 𝑓 𝐤𝑛𝑒𝑤 ≥ 𝑓 𝐤0 allora 𝐤𝑛𝑒𝑤 prende il posto di 𝐤0 con probabilità uguale a

exp −𝑓 𝐤𝑏𝑒𝑠𝑡 −𝑓 𝐤𝑛𝑒𝑤

𝑇

5. 𝑇 è il parametro «temperatura» e decresce nel corso delle iterazioni (probabilità di accettazione alta all’inizio e tende a zero verso la fine)

6. L’algoritmo itera da 2; si ferma quando la temperatura raggiunge una soglia o dopo un certo numero di iterazioni, restituendo 𝐤𝑏𝑒𝑠𝑡 come soluzione ottimale

Performance di SA • La forza di SA è di evitare i minimi locali «saltando» stocasticamente in

regioni alternative dello spazio di ricerca

• Le capacità di ottimizzazione di SA tuttavia sono molto legate alle impostazioni utilizzate – Metodo di scelta della soluzione iniziale (per motivi simili al gradient descent)

– Temperatura di partenza e velocità di raffreddamento (cooling factor)

– Metodo di modifica della soluzione (qual è il modo migliore di generare il vicinato di una parametrizzazione candidata?)

• È possibile mantenere un insieme di soluzioni scartate (tabu search) per evitare di «ricadere» negli stessi minimi locali

20/11/2014

20

Metodi di ricerca globale basati su popolazioni • La natura ci offre molti esempi di risoluzione

di problemi complessi – Generalmente, però, questi metodi richiedono il

contributo di più individui (ad es., evoluzione, super-organismi, ecc.)

• Anziché modificare iterativamente una singola soluzione, potremmo utilizzarne un insieme – Facciamo «collaborare» le soluzioni tra loro, in

modo che esplorino al meglio lo spazio di ricerca

Swarm Intelligence (SI) • I metodi di SI sfruttano il comportamento emergente di gruppi di

individui auto-organizzanti – Generalmente ispirati a sistemi biologici

(ad es., formicai, alveari)

• Il metodo di SI più noto si chiama Particle Swarm Optimization (PSO) – Meta-euristica ispirata a banchi di pesci e stormi di uccelli

– Uno «sciame» di 𝑁 «particelle» (ovvero possibili parametrizzazioni) si muove nello spazio di ricerca, condividendo informazioni e convergendo a una soluzione ottimale (rispetto alla funzione di fitness)

Kennedy et al.,Proc. IEEE Int Conf Neural Netw, 1995

20/11/2014

21

Algoritmo PSO [1/2] • Algoritmo iterativo

1. Si generano le particelle in maniera casuale nello spazio di ricerca

2. Si aggiorna la velocità delle particelle in base a due attrazioni

• Attrazione sociale: le particelle sono attratte dal migliore individuo 𝐠 ∈ ℝ𝑀

• Attrazione cognitiva: ogni particella è attratta dalla migliore posizione in cui si è trovata 𝐛𝑖 ∈ ℝ

𝑀 (dove 𝑖 = 1,… , 𝑁 è l’indice della particella)

3. Si aggiornano le posizioni delle particelle

4. Se il criterio di terminazione non è stato raggiunto si itera da 2

5. Altrimenti si restituisce 𝐠 come soluzione ottimale

Algoritmo PSO [2/2] • L’intensità delle due attrazioni può essere bilanciata agendo su due

fattori, detti fattore sociale e cognitivo

• La velocità cambia con una certa inerzia – L’inerzia può essere mantenuta costante, decrementata linearmente,

determinata sulla base di regole

– La velocità può essere limitata per evitare comportamenti caotici

• Lo spazio di ricerca è limitato – Se una particella «sconfina», viene riposizionata con una strategia decisa

dall’utente (ad es., rimessa «sul bordo», rimbalzo «elastico», ecc.)

20/11/2014

22

Aggiornamento della velocità PSO [1/3]

Migliore soluzione identificata dallo sciame

Migliore soluzione identificata dalla particella

Velocità precedente

Nuova velocità

Posizione corrente della particella Inerzia

Aggiornamento della velocità PSO [2/3] • Consideriamo la stima di 𝑀 parametri (costanti cinetiche), utilizzando uno

sciame di 𝑁 particelle – 𝐱𝑖 ∈ ℝ

𝑀 , 𝑖 = 1,… , 𝑁 posizione della particella 𝑖-esima

– 𝐯𝑖 ∈ ℝ𝑀 velocità della particella 𝑖-esima

– 𝐠 ∈ ℝ𝑀 posizione migliore trovata dallo sciame

– 𝐛𝑖 ∈ ℝ𝑀 posizione migliore trovata dalla particella 𝑖-esima

– 𝑤 ∈ ℝ+ è il valore dell’inerzia

– 𝑐𝑜𝑔 ∈ ℝ+ è il fattore cognitivo: modula l’attrazione della posizione migliore della particella

– 𝑠𝑜𝑐 ∈ ℝ+ è il fattore sociale: modula l’attrazione della posizione migliore identificata dallo sciame

20/11/2014

23

Aggiornamento della velocità PSO [3/3] • Formula di aggiornamento della velocità della 𝑖-esima particella 𝐯𝑖 𝑡 = 𝑤 ⋅ 𝐯𝑖 𝑡 − 1 + 𝐫1 ∘ 𝑠𝑜𝑐 ∘ 𝐠 − 𝐱𝑖 + 𝐫2 ∘ 𝑐𝑜𝑔 ∘ 𝐛𝑖 − 𝐱𝑖

• 𝐫1, 𝐫2 ∈ ℝ𝑛 sono due vettori di numeri pseudo-casuali campionati

uniformemente in [0,1) – Perturbazioni casuali utilizzate per evitare ottimi locali e una migliore

esplorazione dello spazio di ricerca

• Formula di aggiornamento posizione

𝐱𝑖(𝑡) = 𝐱𝑖(𝑡 − 1) + 𝐯𝑖(𝑡)

Algoritmo PSO

Calcola fitness delle

particelle

Aggiorna vettori g e b

Calcola vettori velocità v

Aggiorna posizioni x

Raggiunta condizione di terminazione

?

Posiziona le particelle in maniera casuale

Il vettore g è la soluzione ottimale

VERO

FALSO

20/11/2014

24

Esempi PSO

https://vimeo.com/17407010

PSO • Potente euristica, computazionalmente efficiente, efficace anche

con problemi di grandi dimensioni

• Non esiste (ancora?) un teorema di convergenza – A differenza di SA, non c’è nessuna garanzia di convergere ad un ottimo

globale, nemmeno in termini asintotici (=aspettando un tempo infinito)

• PSO è molto sensibile alle impostazioni – Dimensione dello sciame, inerzia, fattori sociali e cognitivi, dimensione della

regione ammissibile, limite di velocità, metodo di generazione delle particelle

20/11/2014

25

PSO multi-sciame

https://vimeo.com/19133989

PSO per la stima di parametri [1/3] • PSO può essere utilizzata per la stima delle costanti cinetiche delle

reazione

• Molto spesso, gli esperimenti sono ripetuti più volte perturbando il sistema e vedendo come reagisce – Ad es.: diverse concentrazioni iniziali, livelli di irradiazione UV, ecc.

– Possiamo utilizzare una PSO multi-sciame per gestire le diverse condizioni dell’esperimento: la fisica soggiacente il sistema è la medesima, ci aspettiamo che le costanti di reazione stimate siano valide per tutte le condizioni

– Gli esperimenti in ogni condizione iniziale sono generalmente ripetuti più volte per poterne valutare la variabilità statistica: possiamo estendere la funzione di fitness per considerare le ripetizioni

20/11/2014

26

PSO per la stima di parametri [2/3] • Come visto in precedenza, nel caso di sistemi i cui reagenti hanno

concentrazioni elevate il rumore intrinseco è negligibile: si può utilizzare la modellazione deterministica e le equazioni differenziali – In questo caso le particelle corrispondono alle costanti cinetiche

• Nel caso di sistemi caratterizzati da rumore è conveniente utilizzare la modellazione stocastica e relativi metodi di simulazione – In questo caso le particelle corrispondono alle costanti stocastiche

– Al fine di ridurre le fluttuazioni dovute alla stocasticità possiamo ripetere più volte le simulazioni e calcolare una dinamica media

PSO per la stima di parametri [3/3] • Consideriamo il caso stocastico. Una possibile funzione di fitness

per il singolo sciame/esperimento 𝑑 può essere la seguente

𝑓𝑑 𝐜 =1

𝑆 ⋅ 𝐸 ⋅ 𝑄⋅ |𝐘𝑠

𝑑,𝑒(𝑡𝑞) − 𝐗𝑠𝐜(𝑡𝑞) |

𝑄

𝑞=1

𝐸

𝑒=1

𝑆

𝑠=1

– 𝑆: numero di specie misurate, 𝐸 numero di ripetizioni del 𝑑-esimo esperimento, Q numero di campioni nella serie temporale

– 𝐘𝑠𝑑,𝑒(𝑡𝑞) denota la quantità della specie 𝑠-esima, nella condizione

sperimentale 𝑑, durante la 𝑒-esima ripetizione, al tempo 𝑡𝑞

– Per mitigare la stocasticità confrontiamo i target con una dinamica media di più simulazioni denotata da 𝐗𝑠

𝐜(𝑡𝑞)

20/11/2014

27

Serie temporali a tempi discreti delle specie osservate (Y)

Marco S. Nobile - Systems Biology (AA 2013-14) - LM BMB – Università degli Studi di Milano

Simuliamo utilizzando la parametrizzazione c e otteniamo la dinamica simulata X

Marco S. Nobile - Systems Biology (AA 2013-14) - LM BMB – Università degli Studi di Milano

20/11/2014

28

In ogni punto campionato determiniamo la distanza dalla simulazione media 𝐗𝑠𝐜(𝑡𝑞)

nello stesso istante temporale 𝑡𝑞

𝐘𝑠𝑑,𝑒(𝑡𝑞)

𝐗𝑠𝐜(𝑡𝑞)

𝑓𝑑 𝐜

=1

𝑆 ⋅ 𝐸 ⋅ 𝑄 |𝐘𝑠

𝑑,𝑒(𝑡𝑞) − 𝐗𝑠𝐜(𝑡𝑞) |

𝑄

𝑞=1

𝐸

𝑒=1

𝑆

𝑠=1

Marco S. Nobile - Systems Biology (AA 2013-14) - LM BMB – Università degli Studi di Milano

Ripetiamo la procedura per tutti i Q istanti temporali

𝑓𝑑 𝐜 =1

𝑆 ⋅ 𝐸 ⋅ 𝑄⋅ |𝐘𝑠

𝑑,𝑒(𝑡𝑞) − 𝐗𝑠𝐜(𝑡𝑞) |

𝑄

𝑞=1

𝐸

𝑒=1

𝑆

𝑠=1

Marco S. Nobile - Systems Biology (AA 2013-14) - LM BMB – Università degli Studi di Milano

20/11/2014

29

Ripetiamo la procedura per tutte le E ripetizioni dell’esperimento

𝐘𝑠𝑑,𝑒(𝑡𝑞)

𝑓𝑑 𝐜

=1

𝑆 ⋅ 𝐸 ⋅ 𝑄 |𝐘𝑠

𝑑,𝑒(𝑡𝑞) − 𝐗𝑠𝐜(𝑡𝑞) |

𝑄

𝑞=1

𝐸

𝑒=1

𝑆

𝑠=1

𝐘𝑠𝑑,𝑒(𝑡𝑞)

Marco S. Nobile - Systems Biology (AA 2013-14) - LM BMB – Università degli Studi di Milano

Ripetiamo la procedura per tutte le S specie

𝐘𝑠𝑑,𝑒(𝑡𝑞)

𝑓𝑑 𝐜

=1

𝑆 ⋅ 𝐸 ⋅ 𝑄 |𝐘𝑠

𝑑,𝑒(𝑡𝑞) − 𝐗𝑠𝐜(𝑡𝑞) |

𝑄

𝑞=1

𝐸

𝑒=1

𝑆

𝑠=1

Marco S. Nobile - Systems Biology (AA 2013-14) - LM BMB – Università degli Studi di Milano

20/11/2014

30

Risultati PSO per la stima di parametri • La validità del metodo appena descritto è stata mostrata

empiricamente con alcuni modelli di test: – Michaelis-Menten (sx)

– Modello di espressione genica nei procarioti (dx)

Nobile et al., Proc. EvoBIO, 2013

Algoritmi evolutivi • Un’altra categoria di algoritmi basati su popolazioni di soluzioni

candidate sono gli algoritmi evolutivi

• Basati sui principi della selezione naturale – Riproduzione

– Adattamento

– Ereditarietà

– Variazione

– Competizione

20/11/2014

31

Algoritmi genetici [1/2] • Tecnica più nota di computazione evolutiva

– Si sfruttano i meccanismi di selezione naturale per evolvere una popolazione di 𝑁 soluzioni candidate, facendola adattare al problema in esame (=fitness)

Popolazione 𝒫 Popolazione 𝒫′

Popolazione 𝒫′′

Processo di selezione 2) selezione 3) adattamento

Holland, 1975

Algoritmi genetici [2/2] • Il genoma delle soluzioni, nel caso della stima di parametri, è il

vettore ignoto di costanti cinetiche

• I processi di selezione e variazione sono probabilistici e hanno il compito di simulare la pressione evolutiva – Questi processi avvengono mediante applicazione di operatori genetici che

riproducono i meccanismi tipici della riproduzione sessuata: selezione, crossover, mutazione, riproduzione

– Gli operatori genetici modificano il genoma delle soluzioni, generazione dopo generazione

20/11/2014

32

Operatore di selezione • Ad ogni generazione dell’algoritmo genetico vengono scelti alcuni

individui per formare la nuova popolazione – Negli algoritmi genetici standard, la dimensione della popolazione rimane

costante: ripetiamo la selezione fino a quando 𝒫 = |𝒫′|

– La selezione deve privilegiare gli individui migliori

– Ma non deve ricreare una popolazione identica a quella di partenza

– Non si vuole nemmeno selezionare solo gli individui con fitness migliore, pena la perdita di diversità

– Utilizziamo dunque una strategia probabilistica

Esempi di meccanismi di selezione Roulette wheel Ranking Tournament

La probabilità di selezionare un individuo è proporzionale alla fitness

Viene effettuato un ordinamento delle soluzioni dalla peggiore alla migliore La probabilità di selezionare un individuo è proporzionale alla posizione in elenco: per la posizione i-esima P 𝑖 = 𝑖

𝑗𝑁𝑗=1

Vengono estratti K individui dalla popolazione di partenza, in maniera casuale con ripetizione Si sceglie deterministicamente il migliore di questi K individui (torneo) e viene immesso nella popolazione P’

20/11/2014

33

Operatori di riproduzione e crossover • Riproduzione e crossover sono mutualmente

esclusivi – Operatori applicati agli individui della popolazione

selezionata 𝒫′

– Si applica il crossover con probabilità 𝑃𝜒, mentre

la riproduzione è effettuata con probabilità 1 − 𝑃𝜒

– Nel primo caso, il genoma di due individui «genitori» promettenti viene ricombinato per produrre due nuove soluzioni «figlie»

– Nel secondo caso, due individui passano direttamente alla fase di mutazione

Operatore di mutazione • Simula la mutazione casuale del genoma

– Si altera il genoma di una soluzione «figlia» con modifiche casuali

– Utilizzato per introdurre nuovo materiale genetico nella popolazione

– Consente di esplorare nuove regioni dello spazio di ricerca, evitando minimi locali

– Effettuato con probabilità 𝑃𝑚

– Nel caso della stima di parametri, può essere implementato come la perturbazione casuale dei valori (ad es., sommando un numero casuale)

20/11/2014

34

Performance degli algoritmi genetici • Esiste un teorema di convergenza

• Algoritmo computazionalmente più oneroso di PSO e non pensato per gestire soluzioni a valori reali – Necessaria definizione di operatori di mutazione e crossover specifici

• Come in PSO, prestazioni legate alla scelta delle impostazioni – dimensione della popolazione, tipo di selezione, tipo di rappresentazione,

scelta degli operatori e delle probabilità (𝑃𝜒, 𝑃𝑚), criterio di arresto, numero

di iterazioni, ecc.

PSO vs algoritmi genetici • Messi a confronto per il problema della stima di parametri

– Mostrato empiricamente che PSO converge «meglio» e più velocemente

Bes

ozz

i et

al.,

Pro

c. E

voB

IO, 2

012

Drä

ger

et a

l., B

MC

Sys

t B

iol,

2009

20/11/2014

35

Evolution strategy [1/4] • Algoritmo evolutivo basato su una popolazione di individui (come

gli algoritmi genetici) – Si basa anch’esso su evoluzione e adattamento alla fitness

– Differentemente dagli algoritmi genetici, di solito sfruttano solo meccanismi di mutazione gaussiana (no crossover)

– Differentemente dagli algoritmi genetici, la tecnica nasce per ottimizzare vettori a valori reali (ad es., le costanti di reazione)

Rechenberg, 1973

Evolution strategy [2/4] • Meccanismo di mutazione

– Alla soluzione 𝐜𝑖 ∈ ℝ𝑀 viene sommato un vettore di 𝑀 numeri casuali scelti

con distribuzione normale e una deviazione standard fissata a priori (mutation strength)

– Nelle versioni più avanzate, la mutation strength si adatta dinamicamente, in modo da esplorare meglio lo spazio di ricerca

– Nella versione considerata stato dell’arte (covariance matrix adaptation evolution strategy, CMA-ES) la mutation strength è specifica per ogni parametro e tiene conto di eventuali correlazioni tra i parametri

20/11/2014

36

Evolution strategy [3/4]

Evolution strategy [4/4] • Meccanismo di selezione

– La selezione è deterministica e basata sul ranking delle soluzioni

– I parametri possono essere limitati in un sotto-spazio, ma le normali «corrono» su tutto l’asse reale: è possibile escludere le soluzioni non valide o considerarle ugualmente con una certa probabilità

• Evolution Strategy: nessun drawback? – Facile dal punto di vista utente: poche impostazioni e ottime performance

– Computazionalmente molto oneroso e più complesso da implementare (specialmente CMA-ES), ma assicura convergenza rapida

20/11/2014

37

OLTRE LA STIMA DI PARAMETRI

Oltre la stima di parametri • In alcuni casi non solo le costanti di reazione sono ignote: anche la

struttura del modello è sconosciuta – Oppure una parte del modello è ignota o non si è sicuri degli esatti

meccanismi molecolari

• Necessaria una metodologia automatica per il reverse engineering del sistema biochimico – Naturalmente, la ricostruzione del modello richiede anche la sua

parametrizzazione! Il problema del reverse engineering incorpora il problema della stima di parametri

20/11/2014

38

Programmazione genetica • Gli algoritmi genetici consentono l’ottimizzazione di individui

– Rappresentano soluzioni candidate di uno specifico problema

– È possibile identificare soluzioni di intere classi di problemi, introducendo delle variabili nella struttura della soluzione

• Programmazione genetica: invece di evolvere una popolazione di soluzioni candidate, si evolve una popolazione di programmi – Codificati con opportune strutture dati algebriche

(ad es., alberi, grafi)

Reverse engineering di sistemi biochimici • Dato un insieme di osservazioni del sistema, quale modello

consente di riprodurre (mediante simulazione) i fenomeni? – Possibile affrontare questo problema con gli algoritmi evolutivi

– Invece di avere una popolazione di parametrizzazioni, abbiamo una popolazione di modelli

• E i parametri? – Possono co-evolvere con il modello

– Possono essere stimati per ogni modello (richiede nidificazione di algoritmi evolutivi)

Koza et al., Comp Discov, 2007

20/11/2014

39

Esempio di reverse engineering • Nidificazione di Genetic Programming e PSO

Target network Reverse engineered network

Comparison of target data and

simulated dynamics

Nobile et al, Proc. CEC 2013

Indistinguishability • Si può mostrare che per ogni network biochimico possono esistere

molteplici sistemi matematicamente equivalenti – Equivalenti nel senso che la dinamica è identica

– Network equivalenti sono matematicamente (e logicamente) indistinguibili

• Come risolvere? – Attualmente, sembra impossibile senza aggiungere conoscenza di dominio

– Possibile farlo in maniera automatica?

– Idee?

Craciun and Pantea, J Math Chem, 2008

20/11/2014

40

Bibliografia • Particle Swarm Optimization: Kennedy et al. (1995). “Particle Swarm Optimization”. Proc.

of IEEE Int. Conf. on Neural Networks. IV. pp. 1942–1948

• Multi-swarm PSO per stima di parametri di modelli stocastici: Nobile et al. (2012). “A GPU-Based Multi-Swarm PSO Method for Parameter Estimation in Stochastic Biological Systems Exploiting Discrete-Time Target Series”. Proc. EvoBIO 2012, LNCS. Vol. 7264, pp. 1421–1422

• Parameter estimation w/ global optimization: Moles et al. (2003), “Parameter estimation in biochemical pathways: a comparison of global optimization methods”. Genome Research 13: 2467–2474

• Parameter estimation also w/ gradient-based methods: Chou et al. (2009), “Recent Developments in Parameter Estimation and Structure Identification of Biochemical and Genomic Systems”. Mathematical Biosciences. Vol. 19, pp. 57–83

• Parameter estimation w/ various methods: Dräger et al. (2009), “Modeling metabolic networks in C. glutaminicum: a comparison of rate laws in combination with various parameter optimization strategies”. BMC Systems Biology, 3:5

• Indinguishability: Craciun and Pantea (2008), “Identifiability of chemical reaction networks”. Journal of Mathematical Chemistry, 44, pp. 244–259

• Reverse engineering: Nobile et al. (2013), “Reverse engineering of kinetic reaction networks by means of cartesian genetic programming and particle swarm optimization”, Proc. IEEE Congress on Evolutionary Computation (CEC2013), Vol. 1, pp. 1594–1601, 2013

Altre letture • Gradient Descent: Levenberg K. (1944). “A Method for the Solution of Certain Non-Linear

Problems in Least Squares”. Quarterly of Applied Mathematics 2: 164–168

• Simulated Annealing: Kirkpatrick et al. (1983). “Optimization by Simulated Annealing”. Science 220 (4598): 671–680

• Genetic Algorithms: Holland J. (1975), Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor

• Confronto tra GA e PSO per stima di parametri di modelli stocastici: Besozzi et al. (2009), “A Comparison of Genetic Algorithms and Particle Swarm Optimization for Parameter Estimation in Stochastic Biochemical Systems”, Proc. EvoBIO 2009, LNCS, Vol. 5483, pp. 116–127

• Estimation of Distribution Algorithms: Larrañaga et al. (2002), Estimation of distribution algorithms: A new tool for evolutionary computation. Kluwer Acad. Pub., Boston

• Covariance Matrix Adaptation Evolution Strategy: Hansen, N. (2006), “The CMA evolution strategy: a comparing review”, Towards a new evolutionary computation. Advances on estimation of distribution algorithms, Springer, pp. 1769–1776