Analisi QSAR
di
Chiara Cappelletti, Davide Cittaro, Morena Spreafico, Marco Stefani
L’analisi QSAR descritta in questa relazione è stata effettuata utilizzando un insieme di 71
molecole derivate dal 3-amidinofenilalanina selezionate in base alle costanti di inibizione per la
tripsina, la trombina e il fattore Xa.
Un QSAR (Qantitative Structure-Activity Relationship) è un modello matematico che mette in
relazione in termini quantitativi la struttura delle molecole con l’attività biologica, attraverso
l’utilizzo di un insieme di descrittori molecolari, grandezze in grado di rappresentare diverse
proprietà e aspetti delle stesse molecole.
L’obiettivo della nostra analisi è trovare un modello predittivo in grado di descrivere l’attività
biologica dell’insieme dato di leganti nei confronti della tripsina, e di ricavare da questo modello,
se possibile, delle informazioni che permettano di chiarire meglio il meccanismo d’azione.
I descrittori sono stati scelti prendendo in considerazione tutte le proprietà delle molecole in grado
in qualche misura di spiegare l’effetto biologico che vogliamo studiare, o che perlomeno possono
essere coinvolte nel generare l’attività. Le proprietà che abbiamo selezionato descrivono le
caratteristiche steriche, elettroniche, idrofobiche e di flessibiltà delle molecole (Tabella 1).
ID
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Descrittore
Dimensionalità
b_rotN
Weight
a_nN
a_nO
a_nS
Q_VSA_HYD
Q_VSA_POL
KierFlex
apol
dipole
pmi
a_acc
ASA_H
ASA_P
VSA
vol
logP(o/w)
2D
2D
2D
2D
2D
2D
2D
2D
2D
i3D
i3D
2D
i3D
i3D
i3D
i3D
2D
Significato
Numero di legami rotabili
Peso molecolare
Numero di atomi di azoto
Numero di atomi di ossigeno
Numero di atomi di zolfo
Area della superficie idrofobica totale di VdW
Area della superficie polare totale di VdW
Flessibilità molecolare
Somma delle polarizzabilità atomiche
Momento di dipolo
Momento principale di inerzia
Numero di atomi accettori di legami idrogeno
Area totale della superficie idrofobica
Area totale della superficie polare
Area della superficie di VdW
Volume di VdW
Logaritmo del coefficiente di ripartizione ottanolo/acqua
Proprietà descritte
Flessibilità
Steriche
Steriche
Steriche
Steriche
Idrofobicità
Steriche
Flessibilità
Elettroniche
Elettroniche
Steriche
Elettroniche
Idrofobicità
Elettroniche
Steriche
Steriche
Idrofobicità
Tabella 1 - I descrittori
Per valutare l’effettiva capacità dei descrittori selezionati di spiegare il comportamento delle
molecole, e quantificarne l’informazione portata, abbiamo eseguito l’Analisi delle Componenti
Principali (PCA). In Tabella 2 sono riportati i valori di varianza cumulata per le componenti
principali. Si può osservare come le prime due componenti portino più del 75% dell’informazione
totale; è pertanto ragionevole utilizzare solo queste due componenti per le considerazioni
successive.
PCA
1
2
3
4
5
6
7
8
9
10
11
Varianza
52.6667
75.0736
83.3191
90.6100
94.2003
97.0389
98.2187
98.7768
99.2158
99.5344
99.7769
12
13
14
15
16
17
99.8815
99.9449
99.9713
99.9834
99.9938
100.0000
Tabella 2 - Varianza cumulata
•
•
•
•
Sterici
Idrofobici
Flessibili
Elettronici
Figura 1 - Loading PC1-PC2
La Figura 1 mostra i loading dei descrittori, che risultano ben raggruppati secondo le rispettive
proprietà, steriche, elettroniche, idrofobiche e di flessiblità. La prima componente è caratterizzata
soprattutto dal contributo di alcuni dei descrittori sterici (2, 11 e 16); nella seconda componente
sono i descrittori di idrofobicità (in particolare 6, 13 e 17) che appaiono più importanti e decisivi.
I descrittori di flessibilità (1 e 8) sono relativamente vicini nel grafico perché portano lo stesso tipo
di informazione. Per tutti gli altri la separazione non è altrettanto netta, ma è bene sottolineare che
la distinzione sterico/elettronico per i descrittori analizzati non è sempre così marcata:
Q_VSA_POL, per esempio, può essere considerato un descrittore sia sterico che elettronico.
Figura 2 - Score PC1-PC2
L’analisi esplorativa dei dati permette di osservare come si distribuiscono i composti rispetto alla
loro struttura molecolare: con l’Analisi delle Componenti Principali possiamo evidenziare la
presenza di raggruppamenti di molecole con caratteristiche strutturali simili e di individuare
molecole con un comportamento da outlier, ovvero che risultano molto dissimili rispetto alle altre
molecole considerate.
Analizzando la distribuzione degli score possiamo osservare che la maggior parte delle molecole
forma un gruppo compatto al centro del grafico (Figura 2). Alcune molecole si discostano dal
gruppo: le molecole 25, 54 e 57, per esempio, sono in un’area del grafico caratterizzata da una
elevata importanza dei descrittori di tipo idrofobico; le molecole 1, 2, 5, 6, 11 (la più attiva), 16 e
22 si collocano in una regione in cui pesano di più i descrittori che contano il numero degli atomi di
ossigeno, azoto e zolfo.
La molecola 45, interessante perché cocristallizzata con la tripsina nel file PDB del recettore, si
situa in una regione affollata in cui non prevale alcun descrittore particolare.
Dopo l’analisi esplorativa dei dati sono stati calcolati dei modelli utilizzando due differenti metodi
di regressione: il metodo PCR (Principal Component Regression) e il metodo PLS (Partial Least
Square regression), i cui risultati sono riportati di seguito.
PCR
Activity Field
: pki-trypsin
Condition Limit
: 1e+006
Component Limit
: 0
Observations
: 71
Descriptors
: 17
Components Used
: 7
Condition Number : 488359.36
ROOT MEAN SQUARE ERROR (RMSE):0.64878
CORRELATION COEFFICIENT (R2) :0.35025
CROSS-VALIDATED RMSE
:0.74028
CROSS-VALIDATED R2
:0.18674
PLS
Activity Field
: pki-trypsin
Condition Limit
: 1e+006
Component Limit
: 0
Observations
: 71
Descriptors
: 17
Components Used
: 7
Condition Number : 488359.36
ROOT MEAN SQUARE ERROR (RMSE):0.57875
CORRELATION COEFFICIENT (R2) :0.48296
CROSS-VALIDATED RMSE
:0.74521
CROSS-VALIDATED R2
:0.21547
Con il metodo PLS si ottengono valori di correlazione R2 e RCV2 leggermente migliori
ma in entrambi i casi non sono pienamente soddisfacenti. Tutte le molecole
considerate condividono uno scheletro comune a cui vengono aggiunti combinazioni
da due insiemi R1 e R2 di sostituenti: R1 è composto da 14 molecole prevalentemente
policicliche aromatiche con un alto grado di similarità; R2 è un insieme di 43
molecole costruite su poche architetture comuni. Vista l’elevata omogeneità delle
strutture molecolari, quindi, ci aspetteremmo dei valori di correlazione più alti di
quelli ottenuti.
Di seguito sono riportati i grafici che mostrano il confronto tra la risposta
sperimentale e la risposta predetta per i due differenti modelli, applicati al training set
e al test set.
Figura 3 - PCR e PLS sui dati di training
In grafici di questo tipo le molecole che più si avvicinano alla bisettrice sono quelle
meglio predette; le molecole che si trovano al di sopra hanno una attività biologica
sovrastimata, mentre quelle che si trovano al di sotto hanno un’attività biologica
sottostimata dal modello considerato.
Accanto ad alcune molecole del training set è riportato il numero identificativo,
conforme alla numerazione all’interno del database di origine e sul precedente grafico
degli scores. Le molecole numero 25, 54 e 57 (segnate in rosso) sono quelle per cui
rivestono particolare importanza i descrittori di lipofilia. Nel grafico degli scores si
distinguevano abbastanza bene dal resto del gruppo; la loro diversità strutturale
comunque, come si osserva in Figura 3, è spiegata abbastanza bene dal modello; esse
sono caratterizzate da un grado di attività intermedio.
Le molecole 2, 3, 28, 67, 71, 72 (segnate in nero) sono quelle peggio predette dal
modello. E’ stata paragonata la loro posizione rispetto alle altre molecole nel grafico
degli scores, per valutare se la cattiva predizione potesse essere associata ad una
sostanziale differenza strutturale dalle altre molecole. La molecola 2 effettivamente è
collocata fuori dal gruppo centrale, mentre le molecole numero 3 e 67 sono isolate ma
nella regione compresa tra i due addensamenti principali di molecole.
Infine sono state evidenziate in blu le molecole numero 1, 5, 16 e 66, che sono nelle
regioni del grafico degli scores meno affollate. Si osserva una buona predizione per
tutte loro ad eccezione della molecola 66, la cui notevole differenza strutturale
rispetto alle altre non è stata sufficientemente spiegata dal modello.
I valori di R2 per i modelli PCR e PLS riferiti alle molecole del test set, riportate in
Figura 3, sono rispettivamente 0,59 e 0,49.
Test
8
7.5
pKi Predetto
7
6.5
PCR
6
PLS
5.5
5
4.5
4
4
4.5
5
5.5
6
6.5
7
7.5
8
pKi Sperimentale
Figura 3 - PCR e PLS sui dati di test
Abbiamo deciso di riportare i valori ed i grafici dei due modelli per mostrare come,
nel nostro caso, la PLS sia in grado di descrivere meglio il training set, mentre la PCR
ha una maggiore capacità di predizione sui dati del test set. In ogni caso, consideriamo
più importante la capacità predittiva di un modello, quindi la scelta definitiva cade sul
modello PCR.
Il modello scelto viene definito dalla seguente equazione
pki-trypsin =
6.23757
+0.00026
-0.00430
+0.00029
-0.00023
-0.00005
-0.00717
+0.00752
+0.00008
+0.00043
-0.00010
+0.00019
-0.00023
+0.00141
-0.00912
+0.00270
+0.00033
-0.00061
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
b_rotN
Weight
a_nN
a_nO
a_nS
Q_VSA_HYD
Q_VSA_POL
KierFlex
apol
dipole
pmi
a_acc
ASA_H
ASA_P
VSA
vol
logP(o/w)
L’importanza relativa dei descrittori è mostrata in Tabella 3. Risulta immediatamente
evidente come i primi sei descrittori abbiano un coefficiente relativo di un ordine di
grandezza superiore rispetto agli altri e come, escludendo la flessibilità, essi facciano
parte di tutte le tipologie di descrittori da noi usate.
1.000000
0.676809
0.545686
0.449416
0.360005
0.303888
0.193123
0.036336
0.007932
0.001444
0.001357
0.000339
0.000335
0.000285
0.000151
0.000110
0.000019
pmi
Q_VSA_HYD
ASA_P
Weight
Q_VSA_POL
VSA
ASA_H
vol
apol
logP(o/w)
b_rotN
a_acc
a_nO
a_nN
KierFlex
dipole
a_nS
Tabella 3 - PCR: importanza relativa dei descrittori
Dati i risultati così poco soddisfacienti, abbiamo pensato ad un metodo alternativo per
la scelta dei descrittori, basata su considerazioni puramente statistiche. Abbiamo
scelto 75 dei descrittori definiti da MOE in modo da ottenere una ripartizione
proporzionata delle caratteristiche elettroniche, steriche e lipofile delle molecole.
Abbiamo calcolato i valori dei descrittori per le molecole del training set, e per ogni
descrittore abbiamo ottenuto media, deviazione standard e il loro rapporto, scegliendo
i descrittori con un coefficiente di variazione relativo (rapporto deviazione
standard/media) superiore al 30%, partendo dall’ipotesi (arbitraria), che questi
descrittori ci possano dare maggiori informazioni. L'insieme dei descrittori ottenuti è
così sceso a soli 16 elementi (Tabella 4)
ID
Descrittore
Significato
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
PEOE_VSA+0
PEOE_VSA+2
PEOE_VSA+4
PEOE_VSA-0
PEOE_VSA-1
PEOE_VSA-2
PEOE_VSA-4
E
dipoleX
dipoleY
dipoleZ
a_base
a_don
ASADASA
FCASA-
Somma di V degli atomi con carica compresa nell’intervallo [0.00,0.05)
Somma di V degli atomi con carica compresa nell’intervallo [0.10,0.15).
Somma di V degli atomi con carica compresa nell’intervallo [0.20,0.25).
Somma di V degli atomi con carica compresa nell’intervallo [-0.05,0.00).
Somma di V degli atomi con carica compresa nell’intervallo [-0.10,-0.05).
Somma di V degli atomi con carica compresa nell’intervallo [-0.15,-0.10).
Somma di V degli atomi con carica compresa nell’intervallo [-0.25,-0.20).
Energia potenziale
Coordinata X del momento di dipolo (coordinate esterne)
Coordinata Y del momento di dipolo (coordinate esterne)
Coordinata Z del momento di dipolo (coordinate esterne)
Numero di atomi basici
Numero di atomi donatori nei legami idrogeno
Area della superficie accessibile all’acqua di tutti gli atomi con carica parziale negativa
Valore assoluto della differenza tra ASA+ e ASACASA-/ASA dove CASA- = ASA- * max { qi < 0 }
Tabella 4 - I descrittori rivisti
Con questi descrittori abbiamo prodotto vari modelli utilizzando i metodi PCR e PLS,
variando il numero limite di componenti utilizzate. I valori di R2 e Rcv2 sono
direttamente influenzati dal limite imposto al numero di componenti utilizzate.
PCR
0.70
0.60
valore
0.50
0.40
R2
0.30
R2cv
0.20
0.10
0.00
1
2
3
4
5
6
7
8
9
10 11 12 13 14 15 16
componente limite
Figura 4 - Performance dei modelli PCR
PLS
0.70
0.60
valore
0.50
R2
0.40
0.30
R2cv
0.20
0.10
0.00
1
2
3
4
5
6
7
8
9
10 11 12 13 14 15 16
componente limite
Figura 5 - Performance dei modelli PLS
Dall’analisi dei grafici osserviamo che in entrambi i modelli i valori di R2 crescono
all’aumentare del numero di componenti ritenute, mentre i valori di R2cv hanno un punto di
massimo attorno ad un preciso numero di componenti considerate. Questo comportamento è
dovuto al fatto che aumentando la complessità del modello la capacità descrittiva aumenta
sempre mentre lo stesso non vale per la capacità predittiva che, a causa dell’overfitting, cala
oltre un certo limite di componenti considerate. In particolare il modello PCR raggiunge un
massimo in predizione con 14 componenti, mentre il modello PLS raggiunge un massimo
con 13 componenti.
I migliori modelli, calcolati rispettivamente con i metodi PCR e PLS, sono i seguenti:
PCR
PLS
Activity Field
: pki-trypsin
Condition Limit
: 1e+006
Component Limit
: 14
Observations
: 71
Descriptors
: 16
Components Used
: 14
Condition Number
: 35504.335
ROOT MEAN SQUARE ERROR (RMSE):0.49493
CORRELATION COEFFICIENT (R2) :0.62187
CROSS-VALIDATED RMSE
:0.63081
CROSS-VALIDATED R2
:0.40975
Activity Field : pki-trypsin
Condition Limit
: 1e+006
Component Limit
: 13
Observations
: 71
Descriptors
: 16
Components Used
: 13
Condition Number
: 897548.35
ROOT MEAN SQUARE ERROR (RMSE):0.49434
CORRELATION COEFFICIENT (R2) :0.62278
CROSS-VALIDATED RMSE
:0.63773
CROSS-VALIDATED R2
:0.39805
Non c’è una differenza sostanziale tra i due modelli sia sul training set che sul test set.
Rispetto ai modelli precedenti notiamo un miglioramento della performance che passa da
un R2cv di 0,19 a 0,41 per la PCR e da 0,21 a 0,40 per la PLS. Anche il valore di R2
aumenta notevolmente, passando da 0,35 a 0,62 per la PCR e da 0,48 a 0,62 per la PLS.
Training
8
7.5
pKi predetto
7
6.5
PLS
6
PCR
5.5
5
4.5
4
4
5
6
7
8
pKi sperimentale
Figura 6 - PCR e PLS sui dati di training
Test set
8
7.5
pKi predetto
7
6.5
PCR
PLS
6
5.5
5
4.5
4
4
4.5
5
5.5
6
6.5
7
pKi sperimentale
Figura 7 - PCR e PLS sui dati di test
7.5
8
La capacità di predizione dei due modelli sul test set è valutabile in base ai
coefficienti R2 di seguito riportati
PCR: 0,7192
PLS: 0,7189
Tali valori indicano dei modelli più che soddisfacenti rispetto ai risultati dei modelli
precedenti. Analizziamo nel dettaglio il modello PCR che ha una performance
leggermente superiore al PLS, e che ci permette di fare un confronto diretto con il
miglior modello precedente.
L’equazione del modello è:
pki-trypsin
2.87512
+0.00861
+0.00677
+0.04605
-0.00160
+0.00167
+0.01000
+0.02511
-0.00067
+0.02468
+0.03160
-0.03931
+0.03858
-0.64576
+0.00832
-0.00199
-0.08915
=
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
PEOE_VSA+0
PEOE_VSA+2
PEOE_VSA+4
PEOE_VSA-0
PEOE_VSA-1
PEOE_VSA-2
PEOE_VSA-4
E
dipoleX
dipoleY
dipoleZ
a_base
a_don
ASADASA
FCASA-
L’importanza realtiva dei descrittori è:
1.000000
0.957393
0.867060
0.623861
0.230970
0.142853
0.127369
0.120230
0.114637
0.104425
0.100137
0.096977
0.089519
0.084710
0.072327
0.056106
a_don
ASAPEOE_VSA+4
PEOE_VSA+0
PEOE_VSA-2
dipoleY
PEOE_VSA-4
FCASAPEOE_VSA-1
E
PEOE_VSA+2
dipoleX
dipoleZ
DASA
PEOE_VSA-0
a_base
Tabella 5 - PCR: importanza relativa dei descrittori
Osserviamo che i primi quattro descrittori sono di gran lunga i più importanti, ad essi
competono proprietà elettroniche della molecola. Questo modello dà una certa
rilevanza alla presenza di atomi a carica parziale negativa e quindi alla tendenza a
formare legami idrogeno (a_don e ASA-) con la tripsina; in ordine di importanza
troviamo successivamente proprietà
(PEOE_VSA+4 e PEOE_VSA+0).
di
atomi
a
parziale
carica
positiva
Per capire meglio l’attività biologica descritta dai modelli abbiamo analizzato la tasca
in cui si inserisce l’inibitore. La funzione di “site finder” di MOE conferma che il sito
recettoriale coincide con la regione in cui è cocristallizzato l’inibitore. La tasca è
cosituita dai seguenti aminoacidi:
HIS57 ASP189 SER190 CYS192 GLY193 ASP194 SER195 VAL213 SER214
TRP215 GLY216 SER217 GLY218 CYS219 GLY226 VAL227 TYR228
Lo scaffold del legante si inserisce nella tasca. I due atomi di azoto creano legami
idrogeno con ASP189, SER190, GLY218. GLY218 crea un legame idrogeno anche
con uno degli ossigeni del gruppo solfato del sostituente R2.
Sulla superficie del recettore, HIS57 caratterizza una “protuberanza” idrofilica verso
cui è orientato l'azoto del sostituente R2. Questo sostituente è costituito da un anello
idrofobico comunque orientato verso l'istidina. E' ragionevole supporre che un
sostituente di minori dimensioni comporti una più facile interazione con il recettore.
Questa osservazione potrebbe spiegare l’importanza del descrittore Weight (il quarto
per importanza) nel primo modello discusso e la sua correlazione inversa con l’attività
biologica. Sempre sulla superficie, TRP215 è responsabile di una piccola area
idrofobica verso cui si orienta il sostituente R1. Nella stessa zona altri aminoacidi
estendono la regione idrofobica permettendo di supporre che un sostituente idrofobico
di maggiori dimensioni potrebbe comunque accomodarsi ed anzi aumentare la
capacità di interagire con il recettore. Alla definizione di tale superficie idrofobica
partecipa l'anello aromatico della TYR172 e la catena alifatica di LEU99.
Sup. idrofobica
Sup. idrofilica
Figura 8 - Sito recettoriale con legante
L’immagine mostra come si posiziona il legante nel recettore: in verde lo scaffold
inserito nella tasca sopra descritta, in blu il sostituente R1 orientato verso la regione
idofobica e in giallo il sostituente R2 orientato verso la regione idrofilica.
L’analisi del recettore conferma l’importanza dei descrittori del tipo donatori o
accettori di legami idrogeno presenti in entrambi i modelli visto l’alto numero di
legami idrogeno che si formano tra il recettore e il legante, sia ad opera dello scaffold
che ad opera del sostituente R2.
La presenza nel modello di descrittori di tipo idrofobico/idrofilico può essere
ricondotta alle caratteristiche della superficie del sito recettoriale, costituito sia da una
parte idrofobica che da una idrofilica. Infine una correlazione inversa con il peso
molecolare potrebbe indicare che molecole troppo grosse non si collocano
adeguatamente nel sito recettoriale.
L’analisi del sito ci permette di confermare quindi l’importanza di alcuni descrittori
usati nei modelli.
I due approcci utilizzati in questo lavoro sono molto diversi tra loro: se da un lato
sono stati scelti descrittori piuttosto semplici e interpretabili dal punto di vista
dell’attività biologica, nel secondo sono stati premiati dei descrittori in relazione alle
loro caratteristiche statistiche indipendentemente dal significato biologico. I modelli
ottenuti sono altrettanto diversi: con il primo metodo otteniamo dei modelli
maggiormente interpretabili ma poco predittivi, mentre con il secondo metodo
abbiamo ottenuto dei modelli con performance maggiore ma che non ci danno
informazioni riguardo al meccanismo d’azione, in quanto difficilmente interpretabili.
Il passo successivo dovrebbe essere quello di costruire dei modelli più raffinati sulla
base di queste informazioni, per esempio selezionando maggiormente descrittori
relazionabili alle caratteristiche del sito recettoriale, oppure cercando di scandagliare
altre regioni nello spazio dei descrittori, per scoprire nuove proprietà importanti per la
descrizione dell’attività biologica.