Analisi QSAR

di

Chiara Cappelletti, Davide Cittaro, Morena Spreafico, Marco Stefani

INTRODUZIONE Dato un insieme di molecole congeneri vogliamo quantificare le differenze dell’attività biologica in termini di differenze strutturali. A tale scopo abbiamo utilizzato il metodo di regressione lineare dei Minimi Quadrati Ordinari formulando dei modelli QSAR. Le differenze strutturali vengono rappresentate da un insieme di descrittori molecolari appositamente scelti in modo da rendere il più completa possibile la rappresentazione delle molecole. Con un numero di descrittori grande, probabilmente riusciremo ad ottenere una buona descrizione della realtà in esame, ma il modello quasi sicuramente incorrerà nel fenomeno dell’over-fitting, riducendo anche drasticamente la capacità predittiva per le nuove molecole. Un buon modello mantiene generalmente il rapporto variabili-oggetti inferiore a 5, perciò diventa indispensabile selezionare adeguatamente i descrittori in modo da ottenere un modello con una buona capacità sia descrittiva che predittiva. Un ulteriore vincolo è dato dall’uso del metodo dei minimi quadrati, che impone l’eliminazione delle variabili troppo correlate. L’esplorazione esaustiva di tutte le possibili combinazioni di un insieme realistico (varie decine o centinaia) di descrittori richiederebbe un tempo di calcolo troppo elevato; di fatto, ciò impedisce di utilizzare degli algoritmi deterministici in grado di garantire la soluzione ottimale, e quindi s’impone l’utilizzo di metodi approssimati. La nostra scelta ricade sugli algoritmi genetici che permettono di esplorare in tempi ragionevoli lo spazio delle soluzioni. Ottenuto il modello, un passo molto importante consiste nella validazione del modello stesso, cioè la verifica della sua capacità predittiva nei confronti di molecole che non sono state utilizzate per la costruzione del modello. Se il modello è robusto, il passo successivo consiste nell’interpretazione del ruolo dei descrittori scelti. ALGORITMI GENETICI PER LA SELEZIONE DELLE VARIABILI Gli algoritmi genetici, proposti da Holland nel 1975, sono una strategia di ottimizzazione che si ispira alla genetica delle popolazioni. Il processo di ricerca di una soluzione avviene tramite un particolare bilanciamento tra l’esigenza di esplorare nuove regioni nello spazio delle soluzioni, e quella di sfruttare con una ricerca locale le informazioni a disposizione. La descrizione dei valori che ogni variabile può assumere avviene mediante un codice binario, ove ciascun termine binario costituisce un bit, cioè un numero uguale a 0 o 1. Il linguaggio degli algoritmi genetici ha un lessico appropriato in analogia con la genetica delle popolazioni:

o Gene: ciascuna codifica binaria di un numero (una variabile, un parametro numerico) o Individuo o Cromosoma: l’insieme di più geni, cioè la concatenazione di più bits, ed è

associato ad una particolare soluzione o Popolazione: l’insieme di individui o cromosomi o Riproduzione, Cross-over, Mutazione: rappresentano gli operatori genetici che regolano

il processo di ricerca della soluzione ottima Gli algoritmi genetici sono organizzati in tre fasi:

1. Fase iniziale: viene effettuata una codifica delle variabili mediante il codice binario. Ciascuna condizione viene rappresentata da una stringa di 0 o di 1. Nella selezione delle variabili che ottimizzano un modello di regressione, lo 0 indica l’assenza della variabile nel modello mentre l’1 indica la presenza. Viene estratto un certo numero di cromosomi a caso e per ciascuno di essi viene calcolata la risposta corrispondente; i migliori entrano a far parte della popolazione iniziale.

2. Fase evolutiva: è costituita da una fase di generazione di cromosomi che derivano dall’accoppiamento dei cromosomi presenti nella popolazione e calcolo delle risposte corrispondenti, e da una fase di mutazione dove ciascun

cromosoma può casualmente subire delle mutazioni. In questo caso, viene calcolata la corrispondente risposta. Se le risposte calcolate sono migliori della risposta corrispondente all’ultimo cromosoma della popolazione (il peggiore della popolazione esistente), i cromosomi corrispondenti entrano a far parte della popolazione, nella posizione che loro compete in base alla qualità della risposta. In caso contrario, i cromosomi vengono scartati. Esaminando più nel dettaglio, durante la riproduzione sono estratti dei cromosomi destinati all’accoppiamento, con probabilità proporzionale alla qualità del cromosoma. Durante il cross-over vengono generati i cromosomi derivanti dall’accoppiamento, scambiando i geni dei cromosomi genitori con una probabilità prefissata (cross-over probability). Durante la mutazione alcuni geni cambiano il loro valore da 0 a 1 e viceversa, secondo una determinata probabilità (mutation probability). La fase evolutiva procede fino a che non si raggiunge una condizione che pone termine alla procedura.

3. Fase finale: la procedura viene fermata in base a criteri prestabiliti, quali, ad

esempio, il raggiungimento di un numero massimo di iterazioni o quando la popolazione non si rinnova per un certo numero di iterazioni.

Al termine della procedura si ottiene la selezione delle variabili migliori tra quelle che l’algoritmo ha considerato; se tutto lo spazio delle soluzioni possibili è stato indagato, quella proposta è la soluzione ottima.

DESCRIZIONE DEL PROBLEMA Le molecole utilizzate per lo studio sono 14 isomeri delle PCDD (Poli-Cloro-dibenzo-diossina); è importante sottolineare che le molecole analizzate sono strettamente congeneri e differiscono solamente nel numero e nella posizione degli atomi di Cloro, pertanto rappresentano un set molto omogeneo di composti per il quale ci si può aspettare una buona correlazione tra le strutture e la attività biologica. La risposta biologica per la quale viene costruito un modello di regressione è l’affinità di binding. Al fine di calcolare i descrittori molecolari le 14 diossine sono state disegnate con il programma MOE (Molecular Operating Environment) e successivamente sono state minimizzate attraverso i metodi della meccanica molecolare. In particolare è stato utilizzato il campo di forza MMFF94 adatto a descrivere le proprietà di piccole molecole organiche. I DESCRITTORI MOLECOLARI La conformazione spaziale delle 14 molecole minimizzate è il punto di partenza per il calcolo dell’insieme dei descrittori molecolari, effettuato con DRAGON. Il linguaggio con cui descriviamo la realtà non coincide con la realtà stessa: per questo motivo le parole che utilizziamo per esprimerci svolgono un ruolo importantissimo nella rappresentazione della realtà. Nella costruzione di modelli QSAR i descrittori molecolari svolgono lo stesso compito delle parole in un discorso: la qualità del modello finale, che ha lo scopo di rappresentare una parte di realtà, nella fattispecie l’attività biologica dei composti analizzati, dipende strettamente dai descrittori utilizzati. Essi sono il modo attraverso il quale l’informazione chimica è trasformata e codificata. Ognuno di essi è un numero che descrive un aspetto strutturale della molecola: più un descrittore sarà in grado di spiegare bene quell’aspetto, ovvero più il linguaggio è accurato e chiaro, più il modello sarà buono.Per questo motivo la scelta dei descrittori riveste un ruolo fondamentale nella formulazione di un modello. Le classi di descrittori utilizzate in questo lavoro sono:

COSTITUZIONALI Sono descrittori molto semplici, indipendenti dalla connettività e dalla conformazione. Indicano il numero di alcuni atomi, o di legami, o la somma di alcune proprietà (volumi, elettronegatività, polarizzabilità, etc.). Sono facili da ricavare ma la descrizione che danno della molecola è piuttosto generica. Tuttavia, ciò non vuol dire che apportino un’informazione poco importante, sono infatti descrittori molto semplici che si ritrovano spesso nei modelli calcolati.

TOPOLOGICI Largamente utilizzati negli studi QSAR, sono descrittori 2D che derivano dal grafo molecolare (solitamente tolti gli atomi di idrogeno) e dalle matrici di connettività. Sono prevalentemente fondati sui concetti di entropia e di contenuto di informazione ed hanno la proprietà di essere invarianti al grafo molecolare, cioè il loro valore è indipendente da come vengono numerati i vertici del grafo. Sono semplici da calcolare ma i loro limiti più evidenti sono legati alla totale mancanza di informazione riguardo alla geometria 3D e agli aspetti conformazionali, e alla non sempre chiara interpretabilità chimica. A dispetto di ciò, si dimostrano spesso estremamente efficaci per il loro alto potere modellante.

WHIM Il calcolo di questi descrittori prevede la centratura della molecola nel centro delle coordinate e, mediante l’analisi delle componenti principali, la determinazione degli assi principali. La rappresentazione della molecola avviene mediante una proiezione degli atomi sui tre assi principali e sul calcolo di alcuni parametri statistici (gli indici WHIM) per ciascuno di questi assi, senza applicare schemi di pesatura oppure pesando attraverso le masse atomiche, i volumi di Van der Waals, l’elettronegatività atomica di Mulliken, la polarizzabilità atomica o gli indici elettrotopologici di Kier. Si ottengono così dei descrittori WHIM direzionali, che, per ogni schema di pesatura, possono poi dar luogo anche ad indici globali, cioè non direzionali.

GETAWAY Sono dei descrittori basati su una matrice d’influenza, chiamata Molecular Influence Matrix, simile alla matrice dei Leverages solitamente usata nella regressione. Un primo set di descrittori, chiamato H-GETAWAY, viene ottenuto utilizzando solo le informazioni provenienti dalla Molecular Influence Matrix, mentre un secondo set, chiamato R-GETAWAY, combina questa informazione con le distanze interatomiche nella molecola.

Le tre tipologie di descrittori sono state trattate separatamente considerando per ciascun insieme i descrittori con una correlazione pairwise inferiore al 95%. Sono state ottenute rispettivamente 39, 11 e 22 variabili. I PARAMETRI DI PERFORMANCE Per ciascuno dei tre insiemi di descrittori abbiamo utilizzato gli algoritmi genetici per selezionare la migliore combinazione di descrittori. Ogni modello è identificato da un cromosoma di lunghezza pari al numero di descrittori per ogni insieme; il numero di descrittori selezionati in ciascun cromosoma può variare tra 3 e 5. La popolazione è costituita da 100 individui e sono state effettuate 10000 iterazioni per ogni insieme di descrittori. La bontà dei modelli è valutata mediante un parametro di fitting, il coefficiente di determinazione R2, e due parametri di validazione, R2

cv e R2

LOO. Il coefficiente di determinazione R² viene definito come:

2

22

)()ˆ(

11ii

ii

yyyy

TSSRSSR

∑∑

−

−−=−=

Questo valore, moltiplicato per 100, rappresenta la varianza percentuale spiegata dal modello in fitting. La radice quadrata del parametro precedente è una quantità nota come coefficiente di correlazione multipla, R, e rappresenta la correlazione tra la risposta sperimentale e i predittori. Oltre a questi parametri, che misurano il fitting, ve ne sono altri che misurano la capacità predittiva del modello. Il parametro definito PRESS (Predictive Error Sum of Squares) è formalmente analogo a RSS, con la differenza che al posto delle risposte calcolate vengono utilizzate le risposte predette.

PRESS = 2/ )ˆ( iii i yy −∑

Analogamente al parametro R², sostituendo RSS con PRESS, viene definita la varianza spiegata dal modello in predizione, ovvero:

Q2 = R2cv = 1-

TSSPRESS

Diversamente da R2, questo parametro presenta un massimo per la complessità ottimale del modello e ridiscende ogni volta che aggiungiamo al modello variabili non predittive. Inoltre Q2 viene esplicitamente valutato rispetto al potere predittivo e non al livello di fitting del modello. E’ perciò un parametro molto importante, soprattutto nel caso in cui il modello presuppone un suo utilizzo per applicazioni successive su dati la cui risposta non è ancora nota. La validazione è stata fatta secondo due diverse tecniche: attraverso il metodo Leave One Out, sintetizzata dal parametro R2

LOO, e attraverso il metodo Leave Many Out, descritta da R2cv.

Entrambe sono tecniche di cross validation; con questo termine si indica una validazione incrociata, in cui un oggetto o un gruppo di oggetti viene utilizzato per la predizione una sola volta. METODO LEAVE ONE OUT Secondo questo metodo, dati n oggetti, si calcolano n modelli in ciascuno dei quali viene escluso un oggetto alla volta: ciascun modello viene calcolato con gli n-1 oggetti restanti e viene utilizzato per predire la risposta, che può essere di tipo qualitativo in un modello di classificazione o quantitativo in un modello di regressione. Lo scarto tra la risposta predetta e quella sperimentale viene accumulato per tutti gli oggetti che a turno vengono esclusi dal modello. Solitamente il modello finale viene costruito considerando comunque tutti gli oggetti, in modo da sfruttare tutta l’informazione disponibile. Questo metodo è l’unico che consente un confronto univoco tra modelli ottenuti con metodi diversi, ed è quello che comporta la minima perturbazione in quanto viene escluso un solo oggetto per volta, e per questo motivo talvolta fornisce risultati troppo ottimistici. METODO LEAVE MANY OUT E’ una generalizzazione del metodo precedente. Infatti i dati vengono suddivisi in G gruppi di cancellazione in modo tale che k oggetti pari a n/G vengono assegnati all’insieme di valutazione e quindi sono esclusi dal training set: il modello, stimato da un training set costituito da n-k oggetti, viene utilizzato per predire la risposta dei k esclusi. I k oggetti vengono esclusi una volta soltanto; una volta esclusi, vengono successivamente reinseriti nel modello e ne vengono esclusi altri k. Rispetto al metodo precedente, questo è più perturbativo, e valuta in modo più severo la predittività, ma i risultati ottenuti dipendono dalla ripartizione degli oggetti all’interno dei gruppi di cancellazione.

La funzione obiettivo che guida la selezione delle variabili negli Algoritmi Genetici è il parametro R2

LOO e utilizza il metodo di regressione lineare dei Minimi Quadrati Ordinari.

RISULTATI Riportiamo il miglior modello a 3, 4 e 5 variabili per ciascun insieme di descrittori. I modelli scelti sono quelli che si comportano meglio per i tre parametri.

R2 R2LOO R2

cv N Costituzionali e Topologici 0.994 0.938 0.984 5 0.943 0.894 0.88 4 Getaway 0.978 0.93 0.891 5 0.921 0.858 0.845 4 0.952 0.912 0.372 3 WHIM 0.924 0.125 -0.921 5 0.638 0.173 0.122 4 0.635 0.335 0.288 3

Tabella 1 – Performance dei migliori modelli calcolati

Il fatto che nel gruppo Costituzionali e Topologici non vengano selezionati modelli a tre variabili ci fa supporre che tre sole variabili non siano sufficienti a raggiungere la qualità dei modelli a 4 e 5 variabili. I descrittori WHIM hanno una performance decisamente peggiore rispetto agli altri due insiemi. Il valore negativo di R2

cv per il modello a 5 variabili significa che i residui del modello sono tanto superiori ai residui calcolati sul valor medio da far preferire il modello Zero, che seglie sempre il valore medio. Il suo alto valore di fitting non è sufficiente a renderlo un buon modello per la pessima capacità predittiva. I modelli con basso numero di variabili sono preferibili perchè più semplici e più facilmente interpretabili, e per la consuetudine di mantenere il rapporto oggetti/variabili superiore a 5, i modelli a 3 variabili sarebbero stati preferibili. Nel nostro caso però il modello a 3 variabili ha un R2

cv così basso da escluderne la scelta. I descrittori Costituzionali e Topologici si comportano complessivamente meglio e quindi saranno analizzati più nel dettaglio.

ID R2 R2LOO R2

cv N Variabili 1 0.994 0.938 0.984 5 1 3 20 23 37 2 0.983 0.844 0.95 5 1 3 20 21 37 3 0.986 0.905 0.933 5 1 3 15 20 37 4 0.943 0.893 0.888 5 3 11 20 25 36 5 0.943 0.894 0.88 4 3 11 25 36 6 0.952 0.823 0.821 4 4 17 25 29 7 0.928 0.861 0.794 4 3 11 20 25 8 0.92 0.812 0.089 4 1 20 25 29

Tabella 2 – Modelli per i descrittori Costituzionali e Topologici

Osserviamo che negli otto modelli selezionati solo 1/3 dei descrittori totali è rappresentato. Il fatto che pochi descrittori ricorrano frequentemente è indice di robustezza dei modelli. Secondo noi il miglior modello è il numero 5, perché ha 4 variabili e perché tra i modelli a 4 variabili ha la migliore performance, sia in fitting che in predizione. L’equazione del modello è la seguente:

y = 36.324 - 0.009 * MSD -328.954 * TIE -10.725 * Xindex + 2.351 * PCR Riportiamo anche l’equazione in forma standardizzata:

y = -2.692 * MSD -1.053 * TIE -1.377 * Xindex + 0.422 * PCR

Risposta predetta - sperimentale

3

4

5

6

7

8

9

3 4 5 6 7 8 9

Valore sperimentale

Valo

re p

rede

tto

Figura 1 – Risposta predetta – sperimentale per il modello scelto

Nel grafico si osserva il confronto tra risposta sperimentale e risposta predetta dal modello. I punti più vicini alla bisettrice sono quelli meglio predetti. Si osserva un buon accordo del modello con i dati sperimentali per tutte le molecole:

Molecola Valore sperimentale Valore predetto 1 8.00 7.73 2 7.15 7.39 3 7.10 7.05 4 6.79 6.96 5 6.92 6.79 6 6.55 5.94 7 6.10 6.03 8 5.96 6.13 9 5.88 5.91 10 5.49 5.27 11 5.19 5.69 12 5.00 5.10 13 4.89 4.86 14 4.00 3.83

Tabella 3 – Valori predetti dal modello

I descrittori coinvolti sono:

• MSD: mean square distance index (Balaban) • TIE: E-state topological parameter • Xindex: Balaban X index • PCR: ratio of multiple path counts to path counts

I descrittori selezionati appartengono tutti alla categoria dei Topologici e con le nostre conoscenze non è possibile interpretare il loro significato riguardo all’attività biologica. L’omogeneità dei descrittori all’interno di ciascuno dei tre gruppi studiati ci offre una visione parziale della realtà analizzata, limitando la comprensione del fenomeno biologico. Un approccio differente potrebbe essere quello di creare dei gruppi “misti”, che tengano conto di informazioni di diverso tipo, per avere una migliore visione d’insieme e un modello predittivo più completo. Nonostante l’utilizzo di descrittori relativamente omogenei, i risultati ottenuti sono decisamente buoni. Questo comportamento potrebbe essere spiegato da un lato con la congenericità dei composti analizzati, dall’altro con l’efficacia degli Algoritmi Genetici nella selezione delle variabili.