UNIVERSITA’ DEGLI STUDI DI MILANO
FACOLTA’ DI FARMACIA
Corso di laurea in Chimica e Tecnologia Farmaceutiche
NUOVO METODO PER IL DOCKING AUTOMATICO
DI LIGANDI CON MACROMOLECOLE
A STRUTTURA 3D NOTA
Relatore: Chiar.mo Prof. Luigi Villa
Tesi di laurea di: Alessandro Pedretti
matricola n° 398885
anno accademico 1994/1995
Stato dell’arte
1. Introduzione
La ricerca tramite banca dati di strutture tridimensionali [1] è ormai uno degli strumenti più validi per fornire informazioni per l'identificazione di composti biologicamente attivi. Tale tipo di
ricerca può essere di due tipi: la ricerca geometrica e il docking molecolare. La ricerca geometrica evidenzia quelle molecole che soddisfano una serie di relazioni geometriche fra specifici insiemi di atomi o di gruppi funzionali (gruppo farmacoforo, tossicoforo o, più in generale,
bioforo). Il docking, invece, ha la funzione di individuare quelle molecole che possono occupare un specifico sito di legame in una macromolecola e di studiarne la distribuzione spaziale.
Quest'ultima tecnica richiede la conoscenza della struttura del sito di legame, che può essere
acquisita tramite tecniche cristallografiche a raggi X, NMR e/o modellazione per omologia con
strutture note. I metodi di docking posizionano un ligando nel sito di legame di una proteina
mediante ottimizzazione della complementarità sterica, idrofobica ed elettrostatica e, in alcuni
casi, mediante la stima dell'energia libera d'interazione (scoring).
L'incremento del numero disponibile di strutture enzimatiche, proteiche e di nucleotidi di potenziale interesse terapeutico determinate sperimentalmente tramite diffrazione di raggi X e
spettroscopia NMR [2] è fonte d'interesse per la ricerca farmaceutica ed ha reso possiblie la
progettazione di nuovi ligandi (“de novo design”). I metodi di “de novo design” hanno lo scopo di generare nuove strutture chimiche per estrapolazione diretta delle informazioni contenute
nella struttura del sito recettoriale [3]. La costruzione di nuove molecole è un processo iterativo che attualmente si basa largamente sull’intuizione chimica del ricercatore e sul docking interattivo tramite l'impiego della computer grafica. Lo sviluppo interattivo di strutture ad attività
1
farmacologica è solitamente guidato dalla visualizzazione delle proprietà con mappe codificate
opportunamente a colori sulle superfici molecolari. Fra tali metodi di visualizzazione bisogna ricordare le mappe di contorni energetici ottenute tramite sonde molecolari [4,5], il potenziale
elettrostatioco molecolare, le mappe di idrofobicità - idrofilia locale, l’utilizzo di griglie di punti
per la valutazione rapida dell'energia [6] e la minimizzazione energetica di semplici strutture [7].
Il calcolo dell’energia può essere condotto mediante metodiche quantomeccaniche, più rigorose ma riservate per ora a sistemi di dimensioni molto limitate, a causa dell’impegno di calcolo e della quantità di memoria fisica richiesti. In alternativa, sono state sviluppate molte metodiche basate su algoritmi di meccanica classica (campi di forze, vedi oltre) che, se da un lato
consentono di riprodurre le geometrie molecolari con accettabile grado di accuratezza,
dall’altro presentano limitate esigenze di potenza di calcolo e di memoria. Pertanto, con tale
approccio, possono essere trattati in modo qualitativamente accettabile anche i sistemi macromolecola - ligando o addirittura proteina - proteina. La minimizzazione dell’energia stessa
procede poi in diverse forme e con metodiche più o meno raffinate (algoritmi “steepest desent”, VAØ4A, gradienti coniugati, ecc) a seconda dello scopo da ottenere. Essa può riguardare l’ottimizzazione della geometria di valenza dei due partners dell’interazione, oppure la
loro geometria torsionale ( in questo senso si apre il problema conformazionale che è a
tutt’oggi il più complesso da ridurre,) oppure, ed è questo il caso più frequente nel design di ligandi bioattivi, minimizzare la sola componente intermolecolare dell’energia del sistema, mantenendo costanti le componenti intramolecolari dei due partners di interazione. In questo caso
lo scopo è ovviamente quello di ottimizzare l’orientazione reciproca, cioè la complementarità
sterica, lipofilica ed elettronica.
Un cenno a parte merita anche il problema dell’evasione da eventuali minimi locali
sull’ipersuperficie di energia potenziale. I tradizionali metodi di minimizzazione a gradiente sono
2
infatti in grado di individuare solo i minimi più prossimi al punto di partenza. Per ovviare a questo inconveniente sono state sviluppate diverse metodiche, alcune delle quali sono già state discusse nella sezione 2 della parte introduttiva.
Attualmente, il “de novo drug design” rappresenta un settore di ricerca “di frontiera” [8,9,10]
in forte sviluppo anche a causa dell’elevatissiomo numero di strutture 3D disponibili nelle banche dati di svariate industrie appartenenti al settore chimico, farmaceutico ed agricolo. Queste
banche dati rappresentano quindi uno strumento insostituibile per ricerche di docking; ad esse
si aggiungono le banche dati commerciali, pubbliche o create da enti accademici come nel caso del Cambridge Structural Database (CSD) [11], presso il quale sono reperibili le strutture
cristallografiche di oltre 50000 piccole molecole organiche.
Molti sistemi di docking automatico sono in grado di selezionare composti teoricamente interessanti [12] a partire da quelli presenti nelle banche dati, assegnando a ciascuno di essi un
punteggio che ne sintetizzi la potenziale abilità d'interazione col sito di legame. Le molecole
dotate di un punteggio più elevato, quindi potenzialmente candidate ad essere nuovi farmaci,
possono essere sperimentate farmacologicamente e dai risultati ottenuti, si può trarre spunto
per nuove modifiche alla struttura.
In molte tecniche di docking automatico, il ligando è posizionato in una macromolecola a conformazione bloccata, poiché la difficoltà di effettuare l'analisi di una macrostruttura conformazionalmente flessibile rende impossibile l'ottenimento di risultati in tempi ragionevoli. Il docking
in strutture rigide permette di individuare solo le conformazioni a bassa energia d'interazione e
non consente di evidenziare quelle molecole che, per sfortunata scelta del conformero, sono
impossibilitate ad interagire col sito recettoriale benché siano potenzialmente attive [12]. Ad
esempio, il metotressato è disponibile in due conformazioni sperimentali presso la banca dati
3
CSD, ma nessuna è simile a quella che si trova nel complesso diidrofolato reduttasi - metotressato cristallizzato [13].
4
2. Metodi proposti
La selezione di situazioni spaziali potenzialmente interessanti dipende dall'attribuzione di un
punteggio (score) capace di evidenziare la bontà dell’interazione ligando - proteina. Infatti,
soltanto le situazioni con il punteggio più elevato vengono prese in considerazione. Per questo
motivo, l'attribuzione del punteggio risulta essere una delle fasi più critiche del processo di
docking.
I fattori che influenzano l'interazione ligando - recettore comprendono i fenomeni elettrostatici
(inclusi i legami d’idrogeno), le forze di Van der Waals e l'interazione idrofobica [14,15]. In
particolare quest'ultima costituisce la “forza motrice” in grado di promuovere la formazione del
legame, mentre i legami d'idrogeno e le interazioni elettrostatiche sono responsabili della specificità [16].
La formazione del legame intermolecolare richiede la desolvatazione del ligando e della macromolecola con contemporanea rottura dei legami d'idrogeno acqua - recettore e acqua - ligando. Poiché molti di questi legami che si spezzano vengono poi riformati fra il ligando e il recettore, la differenza netta di energia libera è spesso prossima a zero [17,18].
I due principali problemi del docking derivano dalla valutazione della struttura ligando - proteina e dal calcolo dell'energia libera di legame. Il computo dell'energia libera di legame è
estremamente complicato e pertanto solo pochi sistemi di docking sono in grado di realizzare
questo. Le funzioni di attribuzione del punteggio delle singole situazioni spaziali generate dall'interazione ligando - recettore, possono essere di vario tipo: si va da calcoli teorici, molto
onerosi dal punto di vista computazionale, a calcoli empirici, molto più semplici e rapidi da
eseguire. Gli algoritmi di docking devono essere molto veloci per poter avere un utilizzo prati5
co nella ricerca di nuovi farmaci in quanto spesso ci si trova di fronte a numerosi potenziali ligandi che devono essere selezionati dagli archivi delle banche dati. Per questo motivo, il calcolo dell'energia non può che essere approssimativo e non può essere utilizzato in maniera affidabile per confrontare ligandi molto diversi fra loro. In ogni caso, i metodi empirici di calcolo
dell'energia sono quelli maggiormente utilizzati per il docking per i problemi di velocità di calcolo precedentemente menzionati, in quanto veloci e in grado di fornire dati sufficientemente
affidabili di energie relative.
La valutazione della componente di solvatazione nel computo dell'energia libera di legame è
molto difficile. I calcoli di meccanica molecolare tramite campo di forze sovrastimano l'importanza dei legami d'idrogeno anche utilizzando una costante dielettrica distanza - dipendente
[19]. Infatti, a causa di tale sovrastima, è stato ossevato che molecole altamente polari come
gli zuccheri, ottengono spesso elevati punteggi per l’interazione con molti siti di legame, perchè
in grado di essere ottimi accettori - donatori di legami d’idrogeno. Una soluzione a questo
problema può essere quella di introdurre nel campo di forze dei fattori di correzione per riprodurre con maggiore accuratezza il legame d’idrogeno [20]. Tale procedimento è però alquanto
laborioso, in quanto la modifica va adeguata ad ogni diverso tipo di legame d’idrogeno, per il
quale devono essere disponibili sufficienti dati sperimentali o derivati da calcoli per ricavare
parametri corretti. Altrimenti, occorre considerare l’energia libera di solvatazione.
I modelli empirici di solvatazione stimano l'energia libera di solvatazione in funzione dell'area
superficiale accessibile al solvente e producono dei risultati qualitativamente ragionevoli che incidono poco sui tempi di calcolo [21-23]. Tuttavia solo pochi programmi di docking utilizzano
attualmente tale approccio. Uno dei più recenti sistemi per valutare la solvatazione è costituito
da AMSOL [24], il quale combina un metodo semiempirico basato sugli orbitali molecolari
con un modello continuo di solvatazione comprendente 150 soluti neutri e 28 ioni. Tale meto6
do è in grado di commettere un errore assoluto di sole 0.7 Kcal/mol per i soluti neutri e di 2.6
Kcal/mol per gli ioni. Purtroppo AMSOL richiede un elevato tempo di calcolo ed è pertanto
applicabile solo a molecole piccole.
Il programma LUDI di H. J. Böhm [25], dedicato al “de novo design”, posiziona piccoli
gruppi funzionali presenti nel proprio archivio nel sito attivo dell'enzima. E' basato su una serie
di regole empiriche che considerano i tipi di atomi e le loro geometrie per i vari tipi d'interazione. Böhm ha definito una semplice funzione empirica per l'attribuzione del punteggio che utilizza il valore di 1.5 Kcal/mol proposto da Fesht [18] per il legame d'idrogeno e 25 cal/(mol.Å)
per l'interazione idrofobica. Lo scopo fondamentale è quello di determinare il peso relativo dei
legami d'idrogeno e dei termini idrofobici nell'interazione ligando - recettore:
Score =
∑
h − bonds
(100. f ( ∆R). f ( ∆α)) + 53 A
A rappresenta l'area di contatto idrofobico espressa in Ångstrom quadrati, i termini f ( ∆R) e
f ( ∆α) sono delle funzioni empiriche i cui valori sono compresi tra zero e uno, ed hanno il
compito di tener conto della distanza e dall’angolo di legame del legame d’idrogeno non ideale. Questa funzione è in grado di attribuire il punteggio alle singole situazioni spaziali in maniera
soddisfacente.
Bohacek e McMartin [26] hanno sviluppato una diversa funzione empirica per l'attribuzione
del punteggio che si basa sulla correlazione fra tipo di atomo del ligando, che interagisce con la
superficie del recettore accessibile al solvente, e tipo di superficie (idrofobica, accetrice o donatrice di legami d'idrogeno). Questa funzione può risultare molto utile per il docking.
7
Kuntz et al. [27] hanno sviluppato il programma DOCK, il quale è da considerarsi il primo
metodo di docking usato nella ricerca di nuovi farmaci. DOCK è nato allo scopo di realizzare
un algoritmo per il “design de novo” dei ligandi da inserire nel sito attivo di una macromolecola
di struttura nota. Questo programma ricerca gli orientamenti del ligando geometricamente
permessi, mantenendo rigida la struttura del ligando stesso e del sito recettoriale. Il metodo utilizzato per tale scopo comprende diverse fasi: la descrizione del ligando e della cavità recettoriale come insieme di sfere, l'orientamento del ligando e l'attribuzione del punteggio all'orientamento.
La cavità del recettore è descritta come un insieme di sfere in contatto con la superficie molecolare [28]. Ciascuna sfera è generata dalla tangente alla superficie molecolare in un punto a
(pertanto il vettore normale alla superficie comprendente questo punto passa per il centro della
sfera) che tocca la superficie in un altro punto b (Fig. 1).
Fig. 1 - Rappresentazione della superficie molecolare secondo Kuntz
8
Vari criteri sono utilizzati per ridurre il numero di sfere ad un valore accettabile. Solo le sfere
più piccole generate per ciascun punto della superficie vengono mantenute e viene considerata
solo una sfera per ciascun atomo del recettore. Viene calcolato l'angolo fra i due vettori passanti per il centro della sfera e i punti A e B della superficie al fine di scartare sfere che occupano “crepacci” di troppo limitate dimensioni. Solo le sfere che hanno il valore dell'angolo al
di sopra di una soglia di 90° sono tenute in considerazione. Infine, le sfere con raggio maggiore
di 5 Å, che si proiettano perlopiù nello spazio occupato dal solvente, vengono scartate. Le
sfere selezionate sono raggruppate in clusters (vedi oltre) caratterizzati dal fatto di essere costituiti da sfere sovrapposte le une alle altre. Generalmente il cluster più grande corrisponde al
sito di legame del recettore, mentre gli altri clusters descrivono le cavità più piccole presenti
nella superficie proteica. La forma del ligando è similmente descritta da sfere che sono generate posizionandosi all'interno della superficie molecolare del ligando. Le sfere sono prodotte
in maniera analoga alla precedente ad eccezione del fatto che la direzione della normale alla
superficie è definita da un punto interno alla struttura del ligando.
Dopo aver ottenuto una rappresentazione della forma del ligando e del sito recettoriale, è necessario confrontare le rispettive catteristiche spaziali e ciò è effettuato accoppiando i gruppi di
sfere del ligando con quelli della proteina. Per definire l'orientamento del ligando nel sito recettoriale, sono necessarie come minimo quattro coppie di sfere. Sfortunatamente non è possibile esaminare tutti gli accoppiamenti possibili fra le sfere poiché il loro numero è elevatissimo. Per Ns sfere del sito recettoriale e per Nl sfere del ligando si ha un numero di accoppiamenti pari ad N:
N=
Ns !
( N s − N l )!
9
Supponendo di avere Ns = 50 e Nl = 25, si hanno più di 1039 possibiltà e pertanto risulta impraticabile ogni forma di analisi sistematica. Allo scopo di superare questo problema, viene utilizzato un metodo di comparazione basato sul campionamento dei possibili confronti. In particolare, la comparazione delle sfere avviene per selezione di una coppia di sfere del sito recettoriale, misurando la distanza fra i due rispettivi centri. Generalmente le coppie sono scelte
partendo da quelle che hanno distanze maggiori. In seguito viene selezionata una coppia di
sfere del ligando, la cui distanza fra i centri deve essere compatibile con quella della coppia del
recettore considerando eventualmente una certa tolleranza. A questo punto viene selezionata
una terza sfera della proteina, vengono calcolate le distanze fra i tre centri ed infine le due distanze così ottenute vengono confrontate con quelle del ligando ottenute analogamente per
selezione di una nuova sfera. Questo processo continua fino a quando non è stato confrontato
un numero specificato di sfere o non è più possibile trovare nuove coppie che soddisfino la
tolleranza introdotta sulla distanza dei centri. Se si verificano le condizioni precedenti, il programma torna indietro e continua a cercare comparazioni alternative. Solitamente questo algoritmo seleziona molti dei possibili confronti.
DOCK [29] definisce le basi per un metodo di ricerca di composti sulla base della pura e
semplice complementarità della forma del sito recettoriale di struttura nota con le caratteristiche
geometriche del ligando. Questo programma è stato modificato per la ricerca di ligandi presenti nelle banche dati utilizzando un'appropriata funzione di assegnamento del punteggio. Ciò
consente di valutare l'interazione di moltissimi ligandi con una determinata macromolecola in
meno di un giorno di calcolo e di selezionare quel ligando che presenta la maggiore complementarietà di forma nei confronti del sito recettoriale. Ciascuna molecola specificata, presente
in archivio, viene confrontata con la proteina considerata, utilizzando il metodo precedente10
mente descritto. Ad ogni ligando è attribuito un punteggio e alla fine della ricerca viene evidenziata la molecola col miglior punteggio. Generalmente una semplice funzione di attribuzione
del punteggio, basata sul confronto della forma del ligando con quella del recettore ,è in grado
di discriminare i composti attivi da quelli inattivi in modo statisticamente significativo.
Da ciascun confronto delle sfere del ligando con quelle del sito recettoriale è possibile ricavare
un potenziale orientamento del ligando stesso. Tale orientamento è ottenuto dal computo della
rotazione e della traslazione del ligando che soddisfa la sovrapposizione delle coppie di sfere
del ligando con quelle del recettore. Il metodo originale utilizzato dal programma DOCK [27]
per attribuire il punteggio a ciascuna posizione trovata, tiene conto esclusivamente della sovrapposizione delle sfere di Van der Waals:
E overlap = ∑i ∑k ( ri + rk − d ik )
ri e rk sono i raggi di Van der Waals degli atomi del recettore e del ligando e dik è la distanza
interatomica fra gli atomi i e k.
Questa formula è stata modificata per meglio descrivere la complementarità delle forme del ligando con quelle del sito recettoriale [29]. Invece di minimizzare in modo indistinto tutte le sovrapposizioni interatomiche, la nuova funzione discrimina fra le sovrapposizioni costruttive ai fini dell'interazione, e quelle che, facendo prevalere l'effetto repulsivo, indeboliscono il legame ligando - proteina. Vengono misurate le distanze fra gli atomi del ligando e quelli del sito recettoriale e se queste ricadono all'interno della distanza ideale di Van der Waals più un certo valore di tollaranza, viene loro attribuito un punteggio positivo. Gli atomi,la cui distanza si allontana maggiormente da quella ideale, hanno un peso inferiore nel computo del punteggio. Se la
distanza è inferiore ad un valore minimo (generalmente alquanto inferiore rispetto alla distanza
11
ideale), viene attribuito un elevato punteggio negativo. Così facendo, molte deviazioni dalla distanza ideale sono comprese nella funzione allo scopo di compensare l'assenza di flessibilità nei
modelli del ligando e della proteina (Fig. 2).
La fase di attribuzione del punteggio di ciascun possibile orientamento del ligando è molto dispendiosa in termini di tempo di calcolo. Al fine di accelerare questa fase, viene precalcolato il
punteggio di ciascun punto di una griglia collocata nel sito recettoriale [30]. A ciascun atomo
del ligando viene attribuito il punteggio ottenuto considerando quello del punto della griglia più
vicino oppure interpolando i valori di più punti della griglia vicini all'atomo considerato. In questa maniera si ottiene un sensibile incremento della velocità di calcolo.
Fig. 2 - Funzione “soft” di attribuzione del punteggio
12
In molti casi può essere necessario considerare non solo la complementarietà geometrica, ma
anche quella elettrostatica e per questo motivo è stata sviluppata un'altra funzione per l'attribuzione del punteggio [31], che impiega i termini non - bonded del campo di forze AMBER
[32]. Per abbattere i tempi di calcolo, viene utilizzata una griglia precalcolata anche per le proprietà elettrostatiche.
Il metodo approssimato utilizzato dal programma non può garantire l’ottenimento di tutte le
comparazioni potenzialmente interessanti. Per questo motivo sono stati sviluppati sistemi alternativi di confronto [30], che si differenzino dal precedente principalmente per due aspetti.
La prima metodica prevede che le sfere vengano suddivise in più sottogruppi (bins) sulla base
della loro reciproca distanza. I vari sottogruppi hanno la funzione di restringere il campo di
scelta per la sfera successiva da aggiungere al confronto in corso. Le dimensioni del sottogruppo possono essere modificate dall'utente e questo permette un maggiore controllo della ricerca nel suo insieme. Un secondo metodo, molto differente rispetto all'originale, è stato sviluppato al fine di poter specificare l'ordine col quale vengono scelte le sfere. Ad esempio, selezionata una sfera di partenza, si può scegliere una sfera del sottogruppo che comprende le
sfere più lontane da quella di partenza, e così via.
Le strutture generate da DOCK non sono ottimizzate, benchè la funzione di assegnazione del
punteggio comprenda i termini non - bonded del force field AMBER, che sono altamente sensibili alle variazioni di distanze interatomiche [32]. Infatti, piccole variazioni di posizione del ligando possono portare ad ampie variazioni nel punteggio calcolato, evidenziando così una
cattiva correlazione fra accuratezza del docking e punteggio. Meng et al. [33] hanno messo a
punto un sistema semplice e veloce per l'ottimizzazione di strutture rigide, che migliora signifi-
13
cativamente il calcolo dell'energia e che orienta il ligando in una posizione comparabile a quella
ottenuta sperimentalmente.
Quando il ligando e il recettore sono entrambi proteine, il problema del numero di strutture generabili è veramente complesso. Per consentire il docking di sistemi proteina - proteina in tempi ragionevoli, è necessario sfruttare i metodi di confronto modificati descritti in precedenza.
Nella fase di generazione dei clusters, le sfere più grandi non sono considerate e ciò permette
di ottenere clusters di dimensioni generalmente più limitate dato che il numero delle comparazioni possibili cresce rapidamente all'aumentare del numero delle sfere. Clusters piccoli, cioè
costituiti da poche sfere, consentono un notevole risparmio sul tempo di calcolo. E' stato evidenziato con diverse prove che è possibile riprodurre gli orientamenti sperimentali per molti
complessi proteina - proteina ed inoltre sono state trovate posizioni di docking alternative.
Un particolare metodo di docking proteina - proteina, in prticolare antigene-anticorpo, è stato
messo a punto da Vakser [34] e si basa su un particolare algoritmo che cerca di superare
l’inaccuratezza dei dati sperimantali disponibili introducendo una certa “tolleranza” nella descrizione della superficie molecolare. Tale algoritmo prevede che la struttura sia del recettore,
sia del ligando venga proiettata su opportune griglie 3D, assegnando a ciascun punto un valore
proporzionale al volume adiacente libero o occupato dagli atomi, a seconda di punti interni o
esterni alla molecola in esame. Per il recettore, tale valore sarà negativo se il punto è interno
alla molecola, mentre sarà positivo se il punto è esterno. Nel caso della griglia del ligando, i
punti interni avranno un valore positivo, mentre quelli esterni ne avranno uno negativo. Moltiplicando punto per punto la griglia del ligando con quella del recettore, si ottengono valori nulli
per punti della griglia del ligando esterni ad esse, valori positivi qualora i punti interni al ligando
corrispondano a punti esterni del recettore ed infine valori negativi nel caso si verifichi una sovrapposizione delle due molecole considerate. Elaborando i valori ottenuti è possibile ottenere
14
un punteggio per ciascuna posizione del ligando ottenuta per rototraslazione di una griglia rispetto all’altra.
Tale tipo di approccio è stato modificato dallo stesso Vakser in funzione della precisione dei
dati cristallografici a disposizione. Infatti, l’utilizzo di griglie tridimensionali con una distanza fra
i punti (grid step) inferiore alla precisione dei dati sperimentali non è rigoroso ed è privo di significato. Un valore di grid step generalmente corretto in rapporto alle tecniche attuali potrebbe essere di circa 7 Å, ma purtroppo con una griglia così rarefatta l’algoritmo precedentemente proposto diventa inefficace e privo di significato fisico.
E’ stato proposto dallo stesso Autore di modificare l’algoritmo associando a ciascun punto
della griglia un valore che non risulti solo funzione del volume libero o occupato circostante al
punto stesso, ma anche del numero di atomi compresi nel cubo di 7 Å di lato delimitato dai
punti adiacenti a quello considerato.
Singh et al. [35] hanno sviluppato il programma GEMINI che, sfruttando le informazioni relative alle conformazioni delle catene laterali contenute in una banca dati di strutture cristallografiche, riesce a prevedere le conformazioni di un ligando proteico posizionato nel sito recettoriale. E’ stato dimostrato che GEMINI è in grado di generare un orientamento comparabile alle
rispettive strutture cristallografiche per una serie di peptidi legati agli enzimi endotiapepsina,
carbossipeptidasi A e termolisina.
Wang [36], prendendo in considerazione il problema assai complesso di ottimizzare il docking
molecolare al fine di poterlo applicare a complessi ligando - recettore completamente proteici,
evidenzia come un tale tipo di approccio possa essere semplificato considerando solo quelle
zone cave di una delle due proteine che consentano l'avvicinamento e l'interpenetrazione tra i
due partners in esame. In questo modo si riduce notevolmente l'estensione dello spazio conformazionale da analizzare. Per ottenere tutto questo, è stato ripreso l'approccio ideato da
15
Connolly [28] che si basa sulla ricerca dei motivi strutturali salienti dal punto di vista della
complementarità (per esempio cavità ed estroflessioni) presenti sulla superficie molecolare delle
due proteine allo scopo di definirne il reciproco orientamento. Una posizione ottimale può essere identificata evidenziando almeno quattro coppie di tali motivi complementari.
Wang, dopo aver constatato la difficoltà di individuare quattro coppie di zone d'interazione,
ipotizza come sia possibile limitare tale ricerca ad una sola coppia di zone, definendo poi
l'esatto orientamento reciproco attraverso procedure dinamiche di docking, come ad esempio
quelle basate sul confronto di griglie molecolari. Per evidenziare i punti critici di contatto, è necessario definire una superficie molecolare quale quella proposta da Richards [37], che si ottiene esplorando la superficie di Van der Waals con una sfera di raggio prestabilito. In questo
modo, si definisce la superficie accessibile a tale sonda sferica e, nel caso in cui il suo raggio
corrisponda a quello dell'ipotetica sfera che approssima le caratteristiche geometriche dell'acqua, viene caratterizzata la superficie accessibile al solvente. Un tale tipo di descrizione può rivelarsi inadatta per il docking proteina - proteina e pertanto Wang propone di utilizzare una
superficie mista, classificando gli atomi in "esterni" ed "interni" sulla base dell'accessibilità nei
confronti del solvente. Gli atomi esterni vengono descritti mediante la superficie accessibile al
solvente, mentre gli atomi interni sono rappresentati attraverso la semplice superficie di Van
der Waals, la quale permette di descrivere, in modo più preciso, le caratteristiche geometrico
- strutturali delle cavità interne che risultano decisive nell'interazione fra le due proteine. Entrambe le strutture proteiche vengono proiettate su una griglia tridimensionale e per ciascun
punto della superficie interna, viene costruita la sfera che interseca la parte interna della proteina, in modo del tutto analogo alla descrizione delle sfere negative di Kuntz [27]. Considerando
i raggi di tali sfere è possibile descrivere le caratteristiche geometriche della cavità (concavità,
convessità, ecc.), mentre considerando il volume di tali sfere è poi possibile risalire al volume
16
delle cavità intramolecolari e, in modo del tutto analogo, è possibile derivare il volume e la forma delle estroflessioni molecolari. La ricerca della zona d'interazione si realizza confrontando
le estroflessioni del ligando con le cavità del recettore sulla base della compatibilità dei volumi e
della complementarietà geometrica. Individuando un' accoppiamento fra zona recettoriale e ligando, si determinano cinque dei sei gradi di libertà conformazionale, ossia i tre gradi traslazionali e due gradi rotazionali. Il terzo grado di libertà rotazionale viene determinato con una procedura di docking sistematico attraverso 360/teta casi possibili (teta = "step size" della ricerca), valutando l'energia d'interazione e scartando i casi in cui si ha sovrapposizione molecolare.
L'orientamento migliore che si ottiene, viene ottimizzato col metodo Monte Carlo. L'intero
processo di ricerca richiede più tempo rispetto a quello proposto da Connolly, anche se risulta
applicabile ad un numero maggiore di complessi soprattutto grazie all'originale approccio per
la descrizione della superficie molecolare.
Da molti Autori [38-40] è stata evidenziata la possibilità di incorporare nella procedura di
comparazione la complementarità elettronica oltre la già citata complementarità sterica. Alla base di questo approccio, detto DOCK mirato, vi è l'idea di limitare il procedimento di confronto degli atomi (o delle sfere) del ligando con le sfere del recettore sfruttando delle proprietà
compatibili. Per esempio certe sfere recettoriali possono essere identificate come siti donatori
- accettori di legami d'idrogeno e pertanto solo gli atomi del ligando di tipo appropriato possono interagire con queste sfere. I vantaggi di questo tipo di approccio consistono sia in un
potenziale aumento della velocità di comparazione (superiore da uno a tre ordini di grandezza,
in funzione della ristrettività dei criteri applicati) sia in un miglioramento dell'orientamento del ligando rispetto a quanto si ottiene con la sola complementarità geometrica fra ligando e sito
recettoriale.
17
Una delle maggiori limitazioni del metodo impiegato dal programma DOCK consiste nel fatto
che le geometrie del ligando e del recettore sono mantenute rigide. Ciò è legato alla metodologia di confronto che si basa su distanze fisse. Infatti, se fosse permessa la variazione delle distanze stesse, il problema combinatoriale relativo ai confronti diverrebbe di complessità proibitiva. Fra i vari metodi sviluppati per superare almeno in parte questo tipo di ostacolo bisogna
ricordare quello proposto da R. L. DesJarlais et al. [41] che sfrutta la possibilità di suddividere
il ligando in più parti rigide ciascuna delle quali viene posizionata a livello del sito recettoriale
utilizzando l’algoritmo standard del programma DOCK. I vari orientamenti ottenuti dal docking dei singoli frammenti vengono analizzati sequenzialmente allo scopo di individuare il modo
di riassemblare le diverse parti nella struttura originaria del ligando, rispettando le appropriate
geometrie di legame. Purtroppo molti ligandi non possono essere suddivisi in un numero ridotto di parti e pertanto questo approccio non rappresenta la soluzione completa al problema
della flessibilità.
Un altro metodo [40] sfrutta il docking di una porzione iniziale del ligando (frammento
d’ancoraggio) nel sito recettoriale seguito da una ricerca conformazionale per la rimanente
parte del ligando. Il frammento d’ancoraggio è posizionato nella macromolecola utilizzando il
programma DOCK, ottimizzando in seguito la complementarietà dei legami d’idrogeno col sito
recettoriale. I vari orientamenti così ottenuti sono raggruppati in clusters; per ciascuno di essi
viene evidenziato un orientamento del punto d’ancoraggio che viene utilizzato come base per la
ricerca conformazionale. Quest’ultima fase consiste nel mantenere fisso il frammento
d’ancoraggio, mentre i legami presenti nella restante parte del ligando vengono ruotati sistematicamente per ottimizzare l’interazione del ligando col sito recettoriale.
I due metodi elencati precedentemente introducono flessibilità esclusivamente nella struttura del
ligando e non in quella del recettore. Leach prende in considerazione la possibilità di ottenere
18
una parziale flessibilità del recettore tramite il blocco conformazionale dello scheletro base
(backbone) della proteina. Solo alcune catene laterali, preventivamente selezionate, sono mobili secondo parametri predefiniti, mentre le restanti sono mantenute fisse. Nel caso del ligando
vengono generate una serie di conformazioni tramite analisi conformazionale, mentre le catene
laterali del recettore sono selezionate ricavandone la conformazione da strutture proteiche
sperimentali ad alta risoluzione[42]. Il ligando è posizionato nel sito recettoriale sfruttando un
metodo modificato del docking con strutture flessibili, descritto in precedenza. Poi vengono
analizzate varie combinazioni di conformeri del ligando con le posizioni delle catene laterali
della proteina utilizzando algoritmi di ricerca combinatoriale che sono in grado di individuare le
strutture caratterizzate da un valore che ricade in un certo intervallo del minimo energetico globale.
Un’accorgimento utile per tener conto della flessibilità dal ligando consiste nell’impiego di algoritmi genetici [43] (vedi oltre) che consentono di ottimizzare l’accoppiamento tra sfere del ligando e del recettore per il fitting. La flessibilità del ligando viene inclusa ammettendo la variabilità delle torsioni nell’ottimizzazione tramite algoritmi genetici. In questo modo, è possibile
eseguire contemporaneamente un’ottimizzazione dell’orientamento e della conformazione del
ligando variando gli angoli di torsione e selezionando le coppie di sfere per il confronto geometrico.
Il programma DOCK può essere utilizzato per la ricerca di composti “testa di serie” (lead
compounds) dai quali trarre spunto per modificazioni successive. Un esempio di tale applicazione è rappresentato dal docking di una serie di composti (circa 10000, provenienti dalla
banca dati CSD) con la proteasi HIV-1 [39]. I composti ai quali sono stati attribuiti i 200 migliori punteggi utilizzando il metodo della complementarietà geometrica, sono stati analizzati
tramite grafica molecolare allo scopo di selezionarli, utilizzando diversi criteri, per esempio
19
quello della semplicità di sintesi. Uno dei composti maggiormente interessanti si è dimostrato il
bromoperidolo, benchè fosse classificato solo cinquantunesimo nella graduatoria dei migliori ligandi. L’aloperidolo, già disponibile sul mercato e strutturalmente correlabile al bromoperidolo, è stato sperimentato e si è rivelato un inibitore debole ma specifico della proteasi HIV-1.
Successivamente sono state ottenute due strutture cristallografiche del complesso fra la proteasi e un derivato dell’aloperidolo e una di queste si è dimostrata comparabile con
l’orientamento generato da DOCK.
La metodica sfruttata da DOCK e dai programmi ad esso correlati ha consentito di ottenere
composti “lead” con modesta attività verso una serie di siti di legame. Tuttavia, in alcuni casi,
una successiva analisi cristallografica ha indicato modalità di legame piuttosto differenti da
quelle suggerite da DOCK. Ciò può essere interpretato innanzitutto con la moderata affinità
dei ligandi esaminati, che consente di ipotizzare modalità di legame non univoche. Inoltre occorre considerare che la conformazione ricavata dall’analisi cristallografica può essere influenzata dall’impaccamento delle molecole nel cristallo, dalla presenza di molecole del solvente di
cristallizzazione e dalla flessibilità del ligando. Nella maggior parte delle applicazioni, DOCK
non tiene conto di tutto ciò; pertanto non è sorprendente che i risultati talvolta non siano buoni.
Ciò non significa che la modalità di legame suggerita da DOCK non sia possibile, ma solo che
non è la più probabile.
Oltre ai sistemi di docking correlati al metodo introdotto dal programma DOCK, sono state
sviluppate altre metodiche che prendono in considerazione il problema del posizionamento di
piccole molecole in strutture macromolecolari tramite particolari funzioni di attribuzione del
punteggio che sintetizza l’interazione ligando - recettore.
Molti metodi di docking si basano sulla rappresentazione dell’ipersuperficie molecolare oppure
della superficie accessibile al solvente come quella proposta da Connolly [28]. Le tecniche
20
che considerano esplicitamente i punti che compongono la superficie molecolare richiedono
tempi di calcolo elevati, mentre i metodi più veloci impiegano una serie di sfere per rappresentare la superficie.
Una particolare classe di programmi di docking utilizza il metodo Monte Carlo per tener conto
della flessibilità conformazionale del ligando, come nel caso di AUTODOCK [44]. Questo
programma calcola l’energia d’interazione fra ligando e recettore utilizzando un campo di forze di meccanica molecolare (per esempio AMBER [32]). Il docking può essere eseguito a
partire da una posizione determinata manualmente oppure casualmente e sfrutta un algoritmo
di simulated annealing per ottimizzare l’energia totale del sistema. Le traslazioni, le rotazioni e
gli angoli di torsione del ligando sono modificati simultaneamente e ciascuna variazione viene
accettata o scartata utilizzando un algoritmo tipo Metropolis [45] il quale consente di evitare i
minimi energetici locali (vedi sezione sperimentale di questa tesi). All’inizio la temperatura è
impostata su valori estremamente alti, ma poi viene progrssivamente diminuita. Lo svantaggio
principale di questo approccio è costituito dalla lentezza di calcolo che limita l’applicazione
esclusivamente allo studio di ligandi di piccole dimensioni e non permette di eseguire ricerche
su un elevato numero di ligandi. Trattandosi di un metodo stocastico, è necessario disporre di
un metodo per discernere le situazioni spaziali geometricamente simili e quindi dotate di energia
d’interazione simile da quelle potenzialmente interessanti. In ogni caso il programma AUTODOCK è in grado di generare degli orientamenti spaziali corretti e plausibili dal punto di vista
sperimentale.
Hart e Read [46] hanno messo a punto un metodo di docking simile a quello sfruttato dal programma AUTODOCK ed hanno ottenuto dei risultati ragionevoli ma in tempi non propriamente brevi. Yue [47] ha applicato un metodo Monte Carlo - simulated annealing per rico-
21
struire il complesso proteasi serinica - inibitore utilizzando un insieme di misure ottenute sperimentalmente tramite raggi X come vincoli per restringere il campo di ricerca.
I metodi Monte Carlo richiedono tempi di calcolo piuttosto elevati. Malgrado ciò, essi rappresentano una buona scelta per il docking di un numero ristretto di molecole se si desidera un
campionamento molto fitto, ma risultano essere inadeguati per la ricerca dei ligandi a partire
dalle strutture presenti nelle banche dati.
Una particolare tecnica di docking consiste nell’effettuare una preselezione degli orientamenti
ligando - recettore tramite accoppiamento di atomi donatori e accettori di legami d’idrogeno e
nell’attribuzione di un certo grado di flessibilità della struttura del ligando. Smellie, Crippen e
Richards [48] hanno sviluppato un veloce sistema di docking grafico che coinvolge un ligando
conformazionalmente flessibile e un recettore rigido e che è stato applicato al docking di una
serie di triazine con la diidrofolato riduttasi di E. coli complessata col NADPH. Lo schema di
docking comprende coppie donatore - accettore di legami d’idrogeno vincolate da una serie
di distanze minime e massime tra l’atomo accettore e donatore. Per due interazioni simultanee
fra gli atomi recettoriali i,k e gli atomi del ligando j,l, devono essere soddisfatti i seguenti vincoli:
d ik ≤ U jl + δ
d ik ≥ L jl − δ
dik è la distanza fra i punti i e k del sito recettoriale, Ujl e Ljl descrivono le distanze minima e
massima tra gli atomi del ligando j e l che definiscono la flessibilità strutturale del ligando stesso
e δ è un parametro che riassume la flessibilità del sito recettoriale.
22
L’algoritmo utilizzato ricerca tutti i possibili accoppiamenti fra gli atomi donatori e accettori di
legami d’idrogeno e i relativi vincoli di distanza con lo scopo di trovare tutti gli orientamenti del
ligando che permettano almeno cinque legami d’idrogeno col sito recettoriale. In questo modo
vengono evidenziati una serie di orientamenti che soddisfano i vincoli di distanza senza però tener conto delle caratteristiche geometriche del sito di legame. Pertanto, molti orientamenti ottenuti sono stericamente impossibili per effetto della sovrapposizione dei raggi di Van der
Waals del ligando con quelli della proteina e, di conseguenza, vengono scartati tramite un adeguato sistema di controllo.
Questo tipo di approccio si dimostra utile soprattutto nelle fasi preliminari in cui si esegue una
valutazione generica dei possibili siti di interazione caratterizzati da un’alta densità di atomi donatori e/o accettori di legami d’idrogeno.
Kasinos et al. [49] ha messo a punto un metodo di docking simile al precedente in cui il ligando e il sito recettoriale sono rappresentati da una serie di punti ottenuti dalla proiezione di un
vettore della lunghezza di 1.6 Å passante per le coordinate di ciascun atomo del ligando e del
sito recettoriale capace di formare legami d’idrogeno. L’algoritmo di comparazione consite
nella valutazione diretta della complementarità fra i punti calcolati del ligando con quelli della
proteina. Il ligando è poi posizionato nel sito recettoriale sulla base di una sovrapposizione dei
punti d’interazione sede di legame d’idrogeno che sfrutta il metodo dei minimi quadrati per ottimizzare le distanze.
Yamada e Itai [50,51] hanno sviluppato un metodo di docking che considera flessibile la
struttura del ligando tramite una ricerca sistematica degli angoli di torsione, come pure la struttura del recettore tramite minimizzazione energetica d’interazione del complesso iniziale ligando
- proteina utilizzando il campo di forze empirico AMBER [32]. Tale algoritmo è alla base del
programma ADAM, la cui velocità di calcolo è dovuta anche al precalcolo degli orientamenti
23
ligando - sito recettoriale sfruttando la complementarità dei legami d’idrogeno. Tutto questo è
ottenuto tramite il calcolo di una griglia tridimensionale composta da atomi fittizi (dummy
atoms) accettori e donatori di legami d’idrogeno a livello del sito di legame della proteina.
Quindi vengono generati una serie di coppie di atomi accettore/donatore di legami d’idrogeno
ottenuti dalla comparazione degli atomi del ligando con quelli fittizi della griglia. Le varie coppie
vengono raggruppate in sottoinsiemi ciascuno dei quali descrive un particolare orientamento
ligando - macromolecola che viene successivamente ottimizzato variando sistematicamente i
legami mobili della molecola posizionata tramite confronto delle distanze fra gli atomi del ligando e quelli fittizi della griglia recettoriale. I migliori orientamenti ottenuti vengono ottimizzati con
il programma AMBER sia a livello del ligando, sia a livello recettoriale, introducendo così una
completa libertà conformazionale. ADAM è stato utilizzato per il docking del metotressato e
del trimetoprim nella diidrofolato reduttasi dimostrandosi efficace soprattutto perchè il sito attivo dell’enzima considerato è altamente polarizzato con cariche fortemente localizzate, rappresentando così un caso ideale. L’approccio che vede l’impiego della complementarità dei legami d’idrogeno può risultare carente nel caso di siti di legame altamente lipofili come nel caso
della fosfolipasi A2. In questo caso, che rappresenta l’estremo opposto di quello della diidrofolato riduttasi, la scarsa densità di atomi donatori/accettori di legami d’idrogeno nel sito recettoriale non consente di posizionare il ligando in modo adeguato poichè le interazioni idrofobiche non sono contemplate dal metodo applicato.
Un buon metodo dal punto di vista del computo dell'energia d'interazione ligando - macromolecola recettoriale è stato messo a punto da Vedani [52]. Infatti, viene data notevole importanza alla valutazione dei legami d'idrogeno, includendo anche i termini geometrici, nell'espressione dell'energia senza però trascurare le interazioni elettrostatiche, di Van der Waals e quelle
con centri metallici. Il metodo sfruttato per il docking comprende una serie di fasi successive. Il
24
file in formato PDB [53] delle coordinate di partenza del recettore, che deve essere di natura
proteica, viene convertito in un formato appropriato e contemporaneamente vengono individuate tutte le posizioni degli atomi d'idrogeno che presentino un legame polarizzato, compresi
anche quelli delle molecole d'acqua. Al fine di ridurre i tempi di calcolo, la struttura macromolecolare viene ridotta considerando soltanto i residui aminoacidici compresi in un certo intorno
del sito attivo (regione sito - attivo estesa). Alla struttura così semplificata vengono assegnate
le cariche parziali e i parametri del campo di forze (proprio dell'Autore, benché venga utilizzata
la stessa simbologia di AMBER per i tipi atomici) e viene generata la connettività. Il ligando
viene posizionato nel sito recettoriale adattando le catene laterali degli aminoacidi e mantenendo fisso lo scheletro principale della proteina. Questa approssimazione è giustificata dall'incremento di velocità che ne deriva e che permette di analizzare lo spazio conformazionale. Il programma di ottimizzazione del complesso che viene generato sfrutta una minimizzazione eseguita
col metodo dei gradienti coniugati. Ripetendo le fasi precedentemente descritte per più ligandi
è possibile non solo evidenziare quello dotato di migliore energia d'interazione, ma anche di
ottenere per quest'ultimo un modello ottimizzato di interazione grazie alla flessibilità introdotta
sulle catene laterali della proteina.
Nel tentativo di attribuire una giusta importanza all’ idrofobicità, Vakser e Aflalo [54] hanno
proposto un metodo per analizzare il riconoscimento del ligando da parte del recettore (molecular recognition) sulla base delle caratteristiche idrofobiche, data l’importanza che giocano
nell’interazione fra due molecole e nelle equazioni di QSAR classico. Vengono considerati
esclusivamente gli atomi idrofobici, quando si confrontano due superfici molecolari. Questa
approssimazione è giustificata dall’osservazione che nella maggior parte dei complessi ligando
- proteina cristallograficamente noti, il sito attivo presenta una densità di gruppi idrofobici più
alta che altrove. Una spiegazione razionale di tale scelta risiede nel fatto che i gruppi idrofilici
25
sono generalmente solvatati dall’acqua e relativamente liberi di muoversi e quindi meno adatti
all’interazione rispetto a quelli idrofobici. Questi ultimi hanno il vantaggio di non risultare solvatati dall’acqua e di avere un peso maggiore nella determinazione della struttura molecolare.
Prendendo in considerazione esclusivamente gli atomi idrofobici si riesce ad eliminare i falsi
positivi tipici dell’interazione idrofilica, a rendere superfluo lo studio di modelli flessibili nel caso
in cui si assuma che la componente idrofobica sia strutturalmente rigida, ed infine a ridurre il
numero di atomi da confrontare rendendo i calcoli meno onerosi.
E’ possibile combinare la ricerca geometrica in banche dati con il docking di strutture mantenute rigide, come nel caso del programma CLIX [55]. In questo caso particolare, viene utilizzato il programma GRID [4] per la preliminare identificazione delle potenziali regioni di interazione, sfruttando una sonda atomica e la visualizzazione grafica di griglie. La ricerca di strutture tridimensionali tramite banca dati fornisce una serie di ligandi che soddisfano le relazioni
geometriche fra punti farmacoforici che, successivamente, vengono posizionati all’interno del
sito recettoriale. Ciascun orientamento viene ottimizzato mediante meccanica molecolare allo
scopo di ridurre le interazioni repulsive. CLIX attribuisce un punteggio ad ogni orientamento
basato sull’energia di interazione tra punti GRID del recettore e atomi del ligando.
Bacon e Moult [56] hanno sviluppato un sistema basato sul metodo dei minimi quadrati per
ottimizzare la sovrapposizione delle superfici molecolari del sito recettoriale e del ligando. In
particolare, la superficie di Connolly [57] viene rappresentata come un insieme di punti connessi da vettori che costituiscono una struttura reticolare. La rappresentazione del sito recettoriale così semplificata viene sfruttata per un confronto con il ligando secondo il sistema proposto da McLachlan [58] che si basa sulla rototraslazione del ligando e sull’ottimizzazione in base al metodo dei minimi quadrati. L’interazione ligando - macromolecola viene valutata attri-
26
buendo un punteggio sulla base della complementarità geometrica e dell’interazione elettrostatica, nel computo della quale si tiene conto dell’effetto del dielettrico.
Fisher et al. [59] hanno ideato un metodo che esegue il docking fra strutture rigide sulla base
del confronto delle caratteristiche geometriche. La superficie molecolare di Connolly [28] viene rappresentata come una serie di sfere sfruttando lo stesso algoritmo del programma
DOCK. La superficie di ciascuna sfera viene “triangolata” ovvero viene descritta con delle
terne di punti che vengono memorizzate in una matrice per consentire successivamente un rapido accesso ai dati. Il confronto viene effettuato fra i triangoli che descrivono la superficie del
ligando e del sito recettoriale.
Walls e Sternberg [60] hanno sviluppato un metodo che consente di realizzare un docking fra
strutture rigide tra antigeni e anticorpo, sfruttando un algoritmo basato sulla complementarietà
di superficie descritta da un apposito “soft potential”. Essendo la struttura del sito anticorpale
non sperimentalmente nota, ma costruita tramite modelling, è presumibile che contenga errori
ed imprecisioni. Pertanto il metodo di docking impiegato doveva avere una certa tolleranza agli
errori e considerare la flessibilità del sito che ne deriva. Le superfici dell’antigene e
dell’anticorpo vengono rappresentate come una serie di porzioni di superficie poste a intervalli
di 0.5 Å. Le approssimazioni introdotte nell’elaborazione di tali superfici permettono di evidenziarne le principali caratteristiche, ignorando i dettagli considerati trascurabili. Una ricerca globale che implica i sei gradi di libertà rototraslazionali viene realizzata al fine di ottenere la reciproca orientazione tra antigene e anticorpo. Ogni orientazione ottenuta viene valutata con un
potenziale di Lennard - Jones modificato, in modo che la minimizzazione ottenuta con il nuovo
potenziale sia meno precisa della precedente e permetta una maggiore tolleranza
nell’ottimizzazione. Questo potenziale di superficie è stato modificato per favorire delle operazioni di docking che ignorino le ampie zone di superficie che risultano sovrapposte nel com27
plesso antigene - anticorpo. I migliori orientamenti ottenuti sono stati ulteriormente selezionati
in base alla componente elettrostatica. Ogni catena laterale è stata classificata come idrofobica, polare, positiva o negativa e rappresenta da una singola sfera. La componente elettrostatica viene calcolata sommando le interazioni delle porzioni sovrapposte, e quindi vengono scartate le strutture la cui sommatoria elettrostatica risulti maggiore di zero. L’accordo tra i dati
sperimentali e i risultati così ottenuti è generalmente buono anche se la scarsità delle applicazioni finora realizzate non consente una valutazione statistica di questo metodo.
Anche Jiang e Kim [61] hanno sviluppato un approccio di “soft docking” basato su strutture
rigide e ne hanno evidenziato l’efficacia su numerosi complessi proteina - proteina e sul complesso metotressato - diidrofolato reduttasi / NADPH. Tali Autori rappresentano la superficie
molecolare utilizzando vettori normali ai punti considerati su di essa e incorporando i punti
adiacenti in cubi capaci di approssimare la rappresentazione superficiale. Successivamente il ligando viene confrontato con i cubi rappresentanti e il sito attivo attraverso una ricerca rototraslazionale sistematica. I modelli ottenuti vengono valutati con un metodo del tutto analogo al
precedente. I risultati ottenuti per il complesso metotressato - diidrofolato reduttasi / NADPH
riproducono con buona approssimazione i risultati cristallografici. Tale metodo richiede un
tempo di calcolo proporzionale alla complessità del ligando e alla densità dei punti considerati
sulle superfici.
Billeter et al. [62] hanno applicato un metodo di ottimizzazione vincolata, noto come
“algoritmo ellissoidale” per realizzare procedure di docking su strutture sia rigide, sia flessibili.
La ricerca dei migliori orientamenti viene ottenuta tramite particolari tecniche di rototraslazione,
sfruttando una serie di vincoli sterici basati sul calcolo dei raggi di Van der Waals degli atomi
del ligando e del recettore. Il ligando è costretto a rimanere nel sito da tali vincoli, ma il suo
movimento è consentito da una certa flessibilità casuale. L’algoritmo ellissoidale permette di
28
ottenere la corretta orientazione del ligando all’interno di una cavità recettoriale, impedendone,
con i vincoli già citati, l’uscita, ma risulta incapace di indirizzare correttamente il ligando verso il
sito. A tal fine Billeter ha usato un semplice potenziale basato sulla superficie di Connolly capace di guidare il ligando nella giusta direzione. Il metodo è stato applicato al complesso metotressato - diidrofolato reduttasi / NADPH ottenendo risultati incoraggianti.
Una particolare nota di riguardo deve essere rivolta ai già citati algoritmi genetici per le loro
potenzialità e per la possibilità di includerli nelle procedure di ottimizzazione del docking. Gli algoritmi genetici, descritti per la prima volta nel 1975 da John H. Holland [63], possono costituire una nuova via per l'ottimizzazione di strutture molecolari. Fondamentalmente, si tratta di
algoritmi stocastici "intelligenti" sviluppati in analogia con la teoria evolutiva di Darwin e adattati
all'ottimizzazione.
Il problema principale dei calcoli di meccanica molecolare e di quantomeccanica è costituito
dal fatto che l'ipersuperficie di energia da esplorare comprende generalmente, oltre al minimo
assoluto, anche molti minimi locali. Applicando gli algoritmi classici (per esempio quelli a gradiente), è possibile raggiungere solo un minimo locale, e più precisamente quello più vicino alla
geometria di partenza [64]. Diverse strategie sono state proposte per superare questo problema e tra queste vanno menzionate il "simulated annealing", metodi di dinamica molecolare,
sistemi di analisi geometrica, la ricerca sistematica tramite una griglia di campionamento più o
meno fitta ed infine i metodi stocastici basati sulla generazione di configurazioni casuali. Fra
questi ultimi sono compresi gli algoritmi genetici. Essi operano su una popolazione di individui
ragionevolmente numerosa (da 10 a 100) chiamati cromosomi. Ciascun cromosoma è costituito da una stringa che codifica i parametri del sistema che si vuole ottimizzare. Un singolo
parametro codificato è chiamato gene. La stringa è generalmente formata da un insieme (array) binario, o di numeri reali, o di altri codici.
29
L'ottimizzazione con algoritmi genetici è un processo scandito dalle fasi successive dell'evoluzione segnate dal passaggio da una generazione a quella successiva. In primo luogo la popolazione attuale è valutata secondo un criterio prestabilito. Nel caso dell'ottimizzazione di strutture
molecolari, tale criterio è rappresentato ad esempio dall'energia del sistema. Sulla base dei risultati ricavati, viene creata la generazione successiva. Gli individui della generazione corrente
con le migliori caratteristiche hanno più probabilità di partecipare alla creazione della generazione successiva rispetto a quelli che presentano prestazioni peggiori. La valutazione delle caratteristiche di una popolazione è l'unica fase degli algoritmi genetici che non è casuale, ma richiede una "conoscenza" dominio - specifica, ovverosia una capacità di discriminare e quantizzare le capacità dei singoli individui.
La creazione di una nuova popolazione, ad eccezione della fase di inizializzazione, è basata
principalmente su tre operazioni: la semplice copia dei geni selezionati, il cross-over e la mutazione. Nell'operazione di cross-over, in particolare, due figli vengono prodotti da due genitori
selezionati. A livello delle coppie di cromosomi appartenenti rispettivamente ai due genitori
vengono parzialmente scambiati i geni (definiti da uno o più punti selezionati casualmente) in
maniera da ottenere delle nuove sequenze di codice. Nel caso in cui fosse utilizzata una codifica binaria, sarebbero scambiati i bit corrispondenti, come nel seguente esempio:
Genitori
Selezione del sito di cross-over
Figli
000000000
000000|000
111111000
111111111
111111|111
000000111
Per quanto riguarda l'operazione di mutazione, essa può essere ottenuta ad esempio per somma ad un gene di un incremento generato casualmente.
30
Gli algoritmi genetici, essendo metodi stocastici, non consentono di prevedere in anticipo il
successo di un singolo calcolo.
Oggi sono disponibili un elevato numero di algoritmi genetici atti a risolvere i più svariati problemi di simulazione ed ottimizzazione fra i quali vi è il docking molecolare.
Fra le tante applicazioni degli algoritmi genetici nel docking, bisogna ricordare il lavoro di Y. L.
Xiao e D. E. Williams [65], i quali hanno unito le potenzialità di questi algoritmi alla rapidità di
calcolo dell'energia d'interazione fra ligando e proteina tramite force field. Secondo questi autori, il metodo messo a punto da loro è molto più efficiente dei metodi puramente energetici
perché la costruzione di una griglia nell'iperspazio di docking, costituita da svariati milioni di
combinazioni geometriche, è fortemente limitata dalla potenza di calcolo dei computers attuali.
In particolare l'energia del sistema è stata valutata mediante un potenziale semplificato basato
su CVFF [73], includendo nel calcolo solo i termini non-bonded.
Il programma sviluppato da questi autori è chiamato GAME (Genetic Algorithm for Minimization of Energy) ed esegue la comparazione di due molecole A e B. La molecola A è mantenuta
fissa, mentre il centro di massa della molecola B è rototraslato liberamente usando sei gradi di
libertà. Questi sei gradi di libertà vengono codificati in una stringa di bit che costituisce un "individuo". Per esempio la stringa S
S = 100110101110021001010011
rappresenta l'insieme di parametri per la rototraslazione del tipo Pi(X, Y, Z, Ax, Ay, Az):
Pi = (4.5, 5.0, 9.0, 120, 100, 60)
31
Allo stadio iniziale viene creata una generazione di Ni individui tramite l'assegnamento di 0 o di
1 a ciascun bit della stringa in modo casuale. Ciò consente di ottenere una serie di individui che
rappresentano situazioni assai differenti del ligando nello spazio di ricerca. I dati immagazzinati
nella stringa vengono usati per rototraslare la molecola A. I vari individui vengono ordinati
sfruttando una funzione di valutazione basata sul calcolo dell'energia di interazione E della molecola A con B. La funzione di valutazione converte il valore di E nel fattore f che sintetizza le
prestazioni dell'individuo in maniera tale che valori negativi di E abbiano f positivo addizionato
ad un valore base F.
Gli algoritmi genetici consentono di evolvere dalla generazione attuale verso una nuova generazione con caratteristiche migliori tramite tre fasi: selezione, riproduzione e valutazione. Nella fase di riproduzione, l’algoritmo genetico utilizzato dal programma GAME sfrutta la generazione
di numeri casuali per selezionare il punto della stringa dove operare il cross-over e/o la mutazione. I migliori individui della nuova generazione vengono comparati con quelli della precedente e, nel caso in cui risultassero energeticamente peggiori dei precedenti, la nuova generazione verrebbe abbandonata. L'utente può definire il numero massimo di nuove generazioni
create il quale, una volta raggiunto, determina l'arresto del programma.
Oltre alla possibilità di sfruttare il docking per il confronto di ligandi depositati in banche dati
con una macromolecola a struttura tridimensionale nota, può essere assai interessante applicare gli algoritmi di docking a metodi di “de novo design” automatici al fine di valutare l’entità
dell’interazione fra recettore e ligando costruito sulla base della complementarità. Nishibata e
Itai [66] hanno realizzato un metodo per il “de novo design” basato sulla creazione di una serie
di strutture direttamente nella cavità del sito recettoriale, connettendo covalentemente degli
atomi, selezionati casualmente, ad una struttura di partenza. Le potenzialità del ligando vengono
valutate tramite il calcolo dell’energia d’interazione col recettore, sfruttando un campo di forze.
32
Il programma utilizzato (LEGEND) è capace di generare un gran numero di ligandi le proprietà
dei quali vengono valutate successivamente grazie al programma di docking LORE. LEGEND
possiede un archivio contenete informazioni per sedici tipi diversi di atomi come l’ibridazione, il
numero e il tipo di legame, la geometria, ecc. Per semplificare il calcolo dell’energia
d’interazione ligando - recettore e quindi per ridurre il tempo di calcolo, viene generata una
griglia tridimensionale a livello del sito di legame. Per ciascun punto della griglia viene calcolata
l’energia d’interazione fra una sonda atomica (C, N, O, H) e gli atomi del recettore, sfruttando
i termini non-bonded del potenziale MM2. Il potenziale elettrostatico di ciascun punto viene
determinato a partire dalle cariche atomiche degli atomi della macromolecola che si affacciano
verso la cavità del sito di legame. Tali cariche vengono attribuite tramite il programma AMBER. La generazione del ligando tramite il programma LEGEND si sviluppa in tre fasi: selezione del primo atomo, generazione degli atomi successivi e completamento della struttura
molecolare.
Il primo atomo selezionato (atomo d’ancoraggio) ha la funzione di fissre l’intera struttura del ligando a quella del recettore. La scelta di questo atomo avviene casualmente a partire dai numerosi eteroatomi della macromolecola che sono in grado di formare legami d’idrogeno. La
selezione del primo atomo del ligando è influenzata dalla sua capacità d’interagire con l’atomo
d’ancoraggio mediante un legame d’idrogeno. Successivamente, i restanti atomi del ligando
vengono generati uno per uno fino al raggiungimento del numero che viene preventivamente
impostato dall’utente. Per ogni atomo da aggiungere, viene scelto automaticamente, tramite
calcolo di un numero casuale, un atomo del ligando generato in precedenza, col quale formare
il nuovo legame covalente. Il tipo di atomo, il tipo di legame e l’orientamento spaziale del nuovo atomo sono selezionati casualmente. In particolare, la posizione dell’atomo è determinata
da un punto scelto casualmente compreso nel cerchio definito dal valore dalla lunghezza e
33
dell’angolo di legame formato con l’atomo d’attacco, entrambi assegnati in accordo con il relativo tipo di atomo descritto dal potenziale MM2. Nel caso in cui la posizione dell’atomo non
fosse accettabile o per sovrapposizione dei raggi di Van der Waals fra nuovo atomo e atomi
precedentemente generati o per uno sfavorevole contributo energetico, il programma scarta la
posizione d’attacco e tenta di trovarne un’altra sfruttando il metodo precedentemente descritto. Se LEGEND dovesse fallire dopo un certo numero di tentativi, l’atomo da aggiungere alla
struttura del ligando viene scartato e si tenta di generarne uno diverso. La struttura generata
tramite una serie di cicli come quello descritto viene successivamente completata aggiungendo
gli atomi di carbonio mancanti agli anelli aromatici frammentari e saturando le valenze rimaste
aperte con atomi d’idrogeno. La carica del ligando viene determinata col metodo di Del Re
[67] e alla fine la struttura viene ottimizzata. Siccome LEGEND è capace di generare un elevato numero di ligandi, le loro strutture vengono classificate sulla base di parametri energetici,
sfruttando il programma LORE.
In conclusione, i metodi di docking precedentemente descritti si basano essenzialmente o sulle
interazioni idrofiliche (legame d’idrogeno, legami ionici) o su interazioni idrofobiche (fattori sterici e interazioni di non legame). Entrambi gli approcci forniscono buoni risultati nei metodi citati. I metodi considerati si differenziano nettamente anche nella velocità di calcolo e, per alcune applicazioni, gli algoritmi più complessi risultano impraticabili. Le tecniche di docking dovrebbero essere supportate da metodiche di analisi di archivi di strutture capaci di realizzare
rapide ricerche delle molecole che soddisfino alcuni elementari requisiti di complementarità con
il recettore considerato. L’elemento fondamentale di tutti questi approcci risiede nel metodo
con cui vengono valutati gli orientamenti calcolati. Molto spesso le funzioni utilizzate risultano
imprecise, pochè molte componenti (in particolare quella entropica e quella di solvatazione)
vengono ignorate. Rimane da stabilire quanto rapidi possano essere metodi capaci sia di valu34
tare il corretto orientamento sia di calcolare un valore di energia libera accurato. Le componenti energetiche e strutturali del docking non si basano quasi mai sugli stessi algoritmi: infatti
strutture ragionevolmente accurate si possono ottenere con metodi standard di meccanica
molecolare incapaci tuttavia di fornire un’adeguata valutazione dell’energia di legame. Metodi
più precisi potrebbero richiedere un tempo di calcolo nettamente maggiore. D’altro canto, le
tecniche di docking costituiscono un tipico esempio di calcolo parallelo e potrebbero essere
velocizzate semplicemente distribuendo l’insieme dei possibili orientamenti da valutare su una
rete di processori. Le molecole d’acqua di cristallizzazione presenti nel sito di legame costituiscono un annoso problema: debbono essere tenute in considerazione durante il docking oppure devono essere ignorate ?
Considerare delle molecole di acqua legate da un’appropriata rete di legami d’idrogeno nel
sito attivo può aumentare l’energia d’interazione a causa della variazione d’entropia di desolvatazione. Questa strategia è stata spesso utilizzata nei metodi di “de novo drug design”, ma
ha avuto scarso successo. I migliori risultati si ottengono confrontando i modelli ottenuti con
entrambi gli approcci. Bacon e Moult [56] sono stati incapaci di posizionare l’ampicillina nel
sito della β-lattamasi senza includere una molecola d’acqua che probabilmente costituisce una
parte integrante del complesso enzima - inibitore.
Un’altra area di futuro interesse consiste nell’adeguato trattamento della flessibilità recettoriale.
Sono stati sviluppati diversi metodi che permettono almeno una certa libertà conformazionale
del ligando, mentre finora solo il lavoro di Leach [40] ha cercato di trattare la flessibilità del recettore, consentendo variazioni conformazionali a un sottoinsieme limitato di catene laterali, ritenute significative nell’interazione.
Un adeguato trattamento della libertà conformazionale del ligando dovrebbe consentire anche
ampi movimenti del recettore fino a una completa riorganizzazione della struttura proteica.
35
Questo è un problema di difficile soluzione, anche se è noto che in numerosi casi l’interazione
ligando - recettore viene accompagnata da ampie variazioni conformazionali della macromolecola. La capacità di predire tali cambiamenti risulterà di fondamentale importanza nella comprensione dei meccanismi di attivazione recettoriale e trasduzione del segnale.
Degno di nota è il tentativo di Vakser [54] di superare il problema della bassa risoluzione delle
tecniche cristallografiche oggi disponibili. Anche se il problema non è stato completamente risolto, in ogni caso il metodo proposto dall’Autore è in grado di fornire informazioni preliminari
riguardanti i siti di interazione probabili di alcune coppie di proteine. Tali informazioni sono
comunque suscettibili di ulteriori raffinamenti qualora risultino disponibili dati di maggiore precisione.
Un numero esiguo di metodi di docking è stato finora utilizzato a fini preditivi di strutture di
complessi ignote; nella maggior parte delle applicazioni è stata solamente valutata la capacità di
riprodurre strutture note cristallograficamente. Un’eccezione è rappresentata dal metodo di
Kuntz [27] che, nonostante la sua semplicità, ha sorprendentemente ottenuto una notevole serie di successi nell’individuare nuove molecole capaci di legarsi al sito recettoriale in esame. Il
successo dei metodi di docking non disgiunto dalla disponibilità di archivi strutturali sempre più
ricchi indica come queste procedure avranno in futuro un ruolo sempre maggiore nella progettazione di nuovi “lead compounds”.
36
Parte sperimentale
1. Introduzione al metodo
BioDock è un programma modulare creato per essere utilizzabile su piattaforme hardware
differenti apportando semplici modifiche. Questa possibilità è garantita dal fatto che il codice
base è stato sviluppato in linguaggio C [68-70]standard (ANSI). Inoltre, l'impiego nelle fasi
preliminari di più piattaforme hardware ha consentito di evidenziare gli aspetti critici che possono precludere o rendere difficoltoso il funzionamento del programma con altri computers. I
sistemi utilizzati nella stesura del codice sono stati: Silicon Graphics IRIS 4D/50, Silicon
Graphics INDIGO, Silicon Graphics INDIGO² [71], Commodore AMIGA (68030/68882) e
PC MS-DOS [72] (in quest'ultimo caso sono stati collaudati i soli programmi accessori e non
BioDock). Per ogni piattaforma sono stati impiegati gli opportuni strumenti per lo sviluppo di
software oltre a pacchetti software commerciali di modellistica molecolare: per il sistema Silicon Graphics si è fatto uso del Developer's Kit della stessa Silicon Graphics (nelle versioni per
IRIX V 3.3, 4.0.5, 5.3), del complesso di programmi della Biosym [73] (Insight II, Builder,
Biopolymer e Discover nelle versioni 2.0 e 2.1) e del programma quantomeccanico semiempirico MOPAC [74] (nella versione 5) per il calcolo delle cariche atomiche; per il sistema
AMIGA si è fatto uso del compilatore GCC (nelle versioni 2.6.3 e 2.7.0) [75] e del debugger
Enforcer [76] (nella versione 37); per il sistema MS-DOS si è fatto uso del compilatore Turbo
C/C++ [77] (nella versione 3.5).
37
2. Principi teorici.
Il programma BioDock è costituito attualmente da due sottoprogrammi: uno destinato al docking vero e proprio, l'altro destinato all'analisi geometrica e alla caratterizzazione dei clusters.
2.1 Algoritmi di docking.
Il principio generale su cui si basa la sezione di docking consiste nel considerare la proteina e il
ligando come strutture rigide e nel far compiere al ligando spostamenti casuali nello spazio. Di
volta in volta, viene valutata, con un potenziale empirico estremamente semplificato, l'energia
intermolecolare d'interazione tra il ligando nella nuova disposizione e il recettore. Il tipo di potenziale impiegato comprende i termini cosiddetti "non - bonded" (di non legame) e verrà discusso in dettaglio nella sezione di meccanica molecolare di questa tesi. Vengono scartate a
priori, senza valutare l'energia, e pertanto in modo molto rapido, le situazioni in cui si verifica
un'eccessiva sovrapposizione fra sfere di van der Waals degli atomi di ligando e recettore,
determinabile per semplice calcolo delle distanze interatomiche. Le situazioni favorevoli (frames), invece, vengono memorizzate in un opportuno file chiamato di traiettoria (trajectory),
che consente la successiva analisi e manipolazione delle strutture molecolari del complesso ligando - recettore ottenute. Come sistema di riferimento, BioDock utilizza una terna di assi
cartesiani, rispetto ai quali vengono eseguite tutte le operazioni di rototraslazione.
Il movimento della molecola da posizionare è realizzato attraverso la variazione dei sei gradi di
libertà del sistema: traslazione lungo gli assi X, Y ,Z e rotazione sugli assi X, Y, Z. (In teoria
sarebbero necessari ben nove assi di rotazione corrispondenti rispettivamente a ciascun direttore che determina la posizione reciproca nello spazio degli assi X, Y, Z. In questa modo si
avrebbero ben dodici gradi di libertà che, tuttavia, è possibile ridurre a sei imponendo le condi38
zioni di perpendicolarità degli assi e quindi anche dei piani XY, XZ e YZ. Tutto questo porta
ad una riduzione significativa del numero di parametri da considerare per generare un frame).
Per ciascun grado di libertà viene generato un numero casuale compreso in un intervallo definibile dall'utente e successivamente si procede alla rototraslazione secondo lo schema sotto riportato.
Supponendo di avere un ligando di N atomi, per ciascun atomo sarà necessario:
- portare il sistema di riferimento dall'origine degli assi al centro di rotazione (vedi oltre):
X = x - xm
Y = y - ym
Z = z - zm
- ruotare di un angolo γ attorno all'asse Z e cioè sul piano XY:
T = X cos γ - Y sin γ
Y = X sin γ + Y cos γ
X=T
- ruotare di un angolo α attorno all'asse X e cioè' sul piano YZ:
T = Y cos α - Z sin α
Z = Y sin α + Z cos α
Y=T
39
- ruotare di un angolo β attorno all'asse Y e cioè' sul piano XZ:
T = X cos β - Z sin β
Z = X sin β + Z cos β
- riportare le coordinate al sistema di riferimento originale ed eseguire la traslazione:
X = T + x m + ∆x
Y = Y + ym + ∆y
Z = Z + zm + ∆z
In queste relazioni x, y, z sono le coordinate originali dell'atomo generico; X, Y, Z sono le
nuove coordinate; x m, ym, zm sono le coordinate del centro di rotazione; α, β, γ sono le rotazioni attorno gli assi X, Y e Z; ∆x, ∆y, ∆z sono le traslazioni lungo gli assi X, Y, Z; T è una
variabile in cui sono riportati temporaneamente i risultati intermedi del calcolo (coordinata
temporanea). I valori di α, β, γ, ∆x, ∆y, ∆z possono essere impostati dall’utente o determinati
automaticamente dal programma a seconda delle modalità operative (vedi oltre).
Per quanto riguarda il centro di rotazione, esso corrisponde generalmente al centro di massa
del ligando (condizione di default di BioDock), anche se, in casi particolari, può essere definito
a piacimento da parte dell'utente.
Una volta eseguita la rototraslazione, è necessario controllare che non si abbia interazione
repulsiva tra gli atomi della proteina e quelli del ligando. Questa operazione consiste nel verificare che, dopo lo spostamento, la distanza fra gli atomi della molecola e quelli della proteina
40
non sia inferiore alla somma dei rispettivi raggi di Van der Waals. Considerando due atomi (i,
j), le rispettive coordinate (Xi, Yi, Zi e Xj, Yj, Zj) e i rispettivi raggi di Van der Waals (Ri e Rj),
si avrà:
Ri + R j ≤
(X
i
− Xj
) + (Y − Y ) + ( Z
2
2
i
j
i
− Zj
)
2
Se questa condizione non fosse soddisfatta, allora la situazione creatasi come conseguenza
dalla rototraslazione del ligando sarebbe scartata. In realtà, esiste la possibilità che i raggi di
Van der Waals si sovrappongano leggermente senza che si crei una situazione di repulsione
eccessiva. Si introduce quindi un certo grado di tolleranza nella precedente relazione, che diventa quindi:
Ri + R j − S ≤
(X
i
−Xj
) + (Y − Y ) + ( Z
2
2
i
j
i
−Zj
)
2
Il valore di S (grado di sovrapposizione permesso) è definibile dall'utente ed è compreso generalmente in un intervallo compreso fra 0.1 e 0.3 Å.
Superato il test sulla sovrapposizione dei raggi di Van der Waals, si procede al calcolo dell'energia di interazione fra ligando e proteina considerando esclusivamente l'energia di non legame. Per il lavoro di questa tesi, il force field utilizzato è stato CVFF [73]; tuttavia è allo studio la possibilità di rendere BioDock compatibile con l'uso di altri campi di forze standard o
con modificazioni introdotte dall'utente. I termini di non legame del CVFF considerati sono i
seguenti:
41
 Aij
E NB =
∑ ∑  r
Aij =
Ai . A j
Bij =
Bi . B j
i
j
12
ij
−
Bij
6
ij
r
+
qi .q j 

ε. rij 
i e j sono due generici atomi appartenenti rispettivamente al ligando e alla macromolecola; Ai,
Aj, Bi, Bj sono parametri empirici del force field; Rij è la distanza fra i due atomi i e j; ε è la costante dielettrica; qi e qj sono le cariche atomiche parziali di i e j.
I primi due termini dell’equazione precedente esprimono l'interazione di Van der Waals come
somma fra il termine repulsivo e quello attrattivo, mentre il terzo termine rappresenta l'interazione elettrostatica. Non compare esplicitamente un contributo del legame d'idrogeno al potenziale, in quanto, in prima approssimazione, si può considerare compreso nel termine coulombiano. E' allo studio la possibilità di introdurre anche un termine che esprima l'energia di
desolvatazione, che può influenzare in modo determinante l'energetica del sistema. Per quanto
riguarda il valore della costante dielettrica, è stato considerato valido il modello del dielettrico
continuo del CVFF [73].
La velocità di calcolo può essere sensibilmente aumentata introducendo un valore di "cutoff"
nel computo dell'energia. Tale valore soglia rappresenta la distanza interatomica oltre la quale
si possono considerare trascurabili, almeno ai fini del presente lavoro, le energie di interazione
tra i due atomi in esame. Un valore ottimale sembra essere 14 Å [78], ma è stato verificato
che, nel caso del docking fra la tacrina e la colinoesterasi (vedi oltre), un valore di cutoff pari a
10 Å può essere considerato accettabile. Infatti effettuando il calcolo rispettivamente con tutte
le possibili interazioni intramolecolari ed ignorando quelle oltre i 10 Å, la differenza di energia è
42
pari a sole 1-1.5 Kcal/mol, mentre il tempo di calcolo richiesto viene notevolmente diminuito.
Tale differenza energetica può essere considerata trascurabile date le dimensioni del sistema.
Sempre per incrementare la velocità di calcolo, è possibile realizzare il docking coinvolgendo
solo una specifica regione della macromolecola, trascurando tutti gli atomi posti oltre una certa
distanza dal centro di rotazione. Questa possibilità è utile quando i parametri impostati impongano che il ligando non si muova molto rispetto alla posizione di partenza, in particolare quando si desideri, ad esempio, non allontanarsi eccessivamente da una struttura determinata per
via sperimentale, o da un'ipotesi considerata sufficientemente affidabile.
BioDock consente, inoltre, di eseguire una simulazione completamente casuale senza approdare ad una situazione di minimo, oppure di convergere progressivamente verso un frame di
minima energia. Il primo procedimento è utile per ottenere un certo numero di punti di partenza
favorevoli e diversi fra loro per la successiva ottimizzazione. Il secondo modo di procedere
viene sfruttato per l'ottimizzazione stessa.
Tale possibilità è attuata grazie ad una variazione dinamica degli intervalli (inizialmente impostati
dall'utente) entro i quali possono oscillare i valori random utilizzati per la rototraslazione. In
particolare il programma genera una serie di frames che hanno soddisfatto il test di sovrapposizione dei raggi di Van der Waals (frames utili) e tra questi seleziona quello ad energia di interazione inferiore. Il frame così individuato viene utilizzato come nuova struttura di riferimento
rispetto alla quale vengono eseguite le successive rototraslazioni. Contemporaneamente gli intervalli vengono ridotti con un algoritmo selezionabile dall'utente. Nel caso in cui non venga
trovato alcun frame fra quelli della serie attuale che abbia un’energia inferiore rispetto a quella
del frame di riferimento, il programma continua senza convergere. Allo stato attuale, sono disponibili due modalità di intervento sull'ampiezza delle oscillazioni casuali: una caratterizzata da
43
una funzione di tipo lineare e una caratterizzata da una funzione esponenziale, con cui si ottiene
un decremento più graduale degli intervalli.
Durante l'ottimizzazione si possono introdurre dei fenomeni perturbativi che riducano il rischio
di approdare su un minimo locale. Nel presente lavoro è implementata la possibilità di applicare la seguente procedura tipo Metropolis [45] per la selezione dei frames di riferimento: se il
frame attuale risulta avere energia inferiore rispetto a quello precedente (∆E < 0) , viene selezionato e memorizzato come nuova situazione di riferimento. In caso contrario (∆E > 0), viene
computato il fattore F di Boltzmann:
F =e
−
∆E
RT
in cui R è la costante dei gas e T è la temperatura, contemporaneamente viene generato un
numero casuale compreso tra 0 e 1. Se tale numero è maggiore di F, il frame attuale è usato
come nuova struttura di riferimento, in caso contrario viene memorizzato, ma non utilizzato
come riferimento. Tale procedura è nota come “Boltzmann jump”.
Quando l’algoritmo di ottimizzazione converge verso una struttura molecolare, tutti gli intervalli
di rototraslazione vengono ridotti con la modalità lineare o esponenziale precedentemente impostata.
La procedura termina dopo il calcolo di un numero massimo di frames utili (stabilito
dall’utente), o per raggiungimento della convergenza. Quest'ultima situazione è valutata tramite
il calcolo della deviazione standard dell'energia. Vengono considerati i soli frames compresi
nell'ultima serie utilizzata per la ricerca del minimo. L’uso dell’ espressione
44
s=
∑ (X
i
i
− X )2
n −1
dove s è la stima della deviazione standard, n è il numero di osservazioni, Xi è il valore del
dato i, X è la media, non è ottimale (nel nostro caso Xi rappresenta l’energia del frame i e
X l’energia media degli n frames considerati). Infatti, per il calcolo della deviazione standard
con questo metodo, occorre che sia presente in memoria il valore di energia di ogni singolo
frame. Dal punto di vista dell’ ottimizzazione d’uso delle risorse di calcolo, e invece più vantaggiosa la seguente espressione:
s=
∑X
i
2
i
(∑ X )
−
i
n −1
2
i
n
che consente la memorizzazione di soli due valori (la sommatoria delle energie e la sommatoria
del quadrato delle energie).
2.2 Algoritmi di analisi.
Il grandissimo numero di strutture molecolari generate da BioDock rende necessaria la creazione di una sezione di analisi rapida dei frames, che consenta di eliminare in maniera veloce ed
efficace le configurazioni troppo simili fra loro.
Il numero di comparazioni n da eseguire in funzione del numero dei frames a , escludendo la
situazione di partenza, è dato dalla seguente relazione:
45
n=
1 2 1
a − a −1
2
2
Si supponga, in particolare, di aver simulato un docking completamente random che abbia generato una trajectory di 400,000 frames utili e di voler paragonare tutti i frames fra loro. Questo significa dover eseguire circa 80 miliardi di comparazioni. Se fosse possibile effettuare 100
comparazioni al secondo, sarebbero necessari circa 25 anni di calcolo.
Il primo algoritmo implementato nel corso di questa tesi prevedeva la comparazione diretta di
coppie di frames misurando le distanze fra gli atomi dell'uno e dell'altro, facendone poi la media. Se questo valore è inferiore ad una soglia definibile dall'utente, i due frames sono considerati simili. In particolare si ha:
D = ∑i
∆X i2 + ∆Yi 2 + ∆Zi2
n
D è la distanza media fra tutti gli atomi di due frames; ∆Xi, ∆Yi, ∆Zi sono le differenze fra le
coordinate atomiche di atomi corrispondenti nei due frames; n è il numero di atomi che compongono il ligando.
L' inefficienza di questo algoritmo è data dal fatto che i vari frames sono registrati nella trajectory soltanto riportando le rotazioni e le traslazioni relative rispetto al frame di partenza e
l'energia d'interazione. Tuttavia, per poter comparare i due frames, occorre generare le coordinate cartesiane di entrambi e successivamente valutare le distanze interatomiche tra atomi
corrispondenti come illustrato sopra.
In tal modo, riprendendo l’esempio sopra riportato, il calcolo potrebbe essere eseguito in tre
anni e mezzo da un computer come la Silicon Graphics INDIGO2, che è in grado di effettuare
circa 1000 comparazioni al secondo. E’ ovvio che in queste condizioni la comparazione sa46
rebbe comunque improponibile. E’ stato perciò sviluppato un algoritmo estremamente efficiente in grado di abbattere i tempi di calcolo sulla stessa piattaforma hardware e sullo stesso
esempio da tre anni e mezzo a 40 minuti.. Tale risultato è stato ottenuto eseguendo la rototraslazione dei frames prima di compararli tra loro e sintetizzando le posizioni degli atomi di ciascuna situazione spaziale in tre pseudocoordinate Xf, Yf e Zf, temporaneamente memorizzate in
un file. La successiva analisi non ha più come oggetto i singoli atomi di ciascun frame, ma utilizza solamente le terne di pseudocoordinate. Supponendo di avere un ligando costituito da 30
atomi, per ciascuna comparazione fatta col metodo precedentemente implementato, è necessario misurare la distanza fra trenta coppie di coordinate, mentre, utilizzando il nuovo algoritmo, basta confrontare solo tre coppie di valori.
Per discriminare il posizionamento e l’ orientamento reciproco del ligando nei vari frames che è
stato considerato accettabile per l’applicazione descritta in questa tesi, l’utente definisce la distanza media fra gli atomi di una generica coppia di frames al di sotto della quale i frames in
questione vengono considerati simili. Un valore di tale distanza è pari a 1 Å. Le situazioni considerate simili vengono memorizzate in una matrice binaria, generando così delle “famiglie” di
frames (clusters),per ciascuna delle quali viene salvato il frame dotato di minor energia
d’interazione con la proteina. Gli altri frames del cluster vengono scartati.
47
2.3 Descrizione del programma
2.3.1 Struttura di BioDock
Dal punto di vista della programmazione, l'intero codice del programma è suddivisibile in sezioni "fisiche" e in sezioni "logiche" (Fig. 3). Le sezioni "fisiche" corrispondono a vari files di
codice sorgente di dimensioni limitate che sono stati generati per rendere meno onerosa, in
termini di tempo, la fase di compilazione. Inoltre, la possibilità di disporre di una serie di listati
relativamente brevi consente un più facile debugging in virtù della maggiore leggibilità. Si tenga
presente che BioDock è costituito da più di 5000 righe di codice (60000 caratteri). Pertanto,
se non fosse stata applicata la tecnica che prevede la generazione e la compilazione di moduli
separati, ognuno deputato ad eseguire un numero limitato di operazioni, difficilmente si sarebbe
potuto modificare in maniera agevole il programma.
Le sezioni "fisiche" comprendono le parti di codice vero e proprio (files con estensione .c) e i
cosiddetti "handlers" (files con estensione .h). Questi ultimi, in particolare, contengono informazioni quali tabelle di dati, costanti, variabili preinizializzate, ecc.
Le parti di codice vero e proprio si distinguono in:
Main.c
- inizializzazione generale.
Parser.c
- interpretazione e controllo sintassi del file di input.
Dock.c
- modulo per il docking.
Loader.c
- caricamento ed assegnazione di carica, potenziale e raggio atomico agli
atomi contenuti nel file di coordinate.
Restart.c
- inizializzazione necessaria in caso di restart della simulazione.
Analyzer.c - modulo per l'analisi dei clusters.
48
Sub.c
- sottoprogrammi vari.
Per quanto riguarda gli "handlers" bisogna ricordare:
Ini.h
- strutture, variabili e costanti condivise con tutte le parti di codice.
Pm.h
- tabelle dei parametri del force field, dei raggi di Van der Waals e dei pesi
atomici.
Parser.h
- tabella delle parole chiave per l'interpretazione del file di input.
Analyzer.h - strutture e costanti per l'analisi dei clusters.
49
Fig. 3 - Diagramma di flusso del programma BioDock V1.1
50
Le sezioni "logiche" consentono di suddividere il flusso del programma in diversi passi (steps),
ciascuno dei quali presenta una serie di subroutines in grado di svolgere delle operazioni elementari. Le principali sezioni “logiche” che compongono BioDock sono le seguenti:
- Main
- Docking
- Lettura delle coordinate atomiche
- Assegnazione dei parametri di potenziale
- Generazione del file di restart
- Posizionamento ligando e calcolo energetico
- Memorizzazione su supporto magnetico
- Analyzer
- Lettura delle coordinate atomiche
- Rototraslazione dei frames (Fase 1)
- Analisi geometrica (Fase 2)
- Memorizzazione della nuova traiettoria (Fase 3)
2.3.2 Sezione "Main"
Costituisce il codice che viene eseguito per primo quando viene lanciato BioDock . Ha la funzione di interpretare le opzioni di linea come il nome del file di input e l'attivazione della modalità interattiva che vengono specificate al momento dell'invocazione del programma. Inoltre
questa sezione comprende l'interprete dei comandi che elencati nel file di input e impostati
51
dall’utente per generare la simulazione. In base a tali comandi viene selezionato il modulo opportuno che, allo stato attuale può essere o quello di docking , o quello di analisi.
2.3.3 Sezione “Docking”
Il modulo in questione è quello che in primo luogo carica e valuta la struttura sia del ligando,
sia della macromolecola e successivamente esegue il docking vero e proprio.
Le operazioni eseguite dal programma possono pertanto essere distinte in:
- Lettura delle coordinate atomiche.
Il file contenente le coordinate atomiche deve essere in formato Biosym [73] di tipo “nuovo”
(generato da Insight II versione 2.1 o superiore) o di tipo “vecchio”,generato da versioni di
Insight II precedenti, e deve contenere contemporaneamente, in forma di insieme molecolare
(assembly), le strutture del ligando e della proteina. In questo file vengono elencati per ciascun
atomo le coordinate cartesiane, il tipo atomico, la carica e il numero di residuo
d’appartenenza. In particolare tale numero è importante per individuare gli atomi del ligando.
E’ allo studio la possibilità di utilizzare i files PDB [53] per l’input delle coordinate e il relativo
programma atto a ricostruire la connettività (ossia la matrice che definisce le interconnessioni
atomiche) e ad assegnare i parametri di carica e potenziale.
- Assegnazione dei parametri di potenziale.
A ciascun atomo sono assegnati i valori di A e B che rientrano nell’espressione del calcolo
dell’energia d’interazione, come precedentemente descritto. Oltre a questi,a ciascun atomo
sono assegnati i valori del raggio di Van der Waals da utilizzare nel test di sovrapposizione
52
delle sfere atomiche; contemporaneamente è calcolata la carica totale e il peso molecolare sia
del ligando, sia della macromolecola. Inoltre, per il solo ligando, viene determinato il centro di
massa che eventualmente verrà utilizzato come centro di rotazione per generare nuovi frames.
La formula utilizzata per tale calcolo è la seguente:
Xm =
Ym =
Zm =
∑ M .X
i
i
i
PM
∑ M .Y
i
i
i
PM
∑ M .Z
i
i
i
PM
PM = ∑i M i
Xm, Ym, Zm sono le coordinate del centro di massa; Xi, Yi, Zi sono le coordinate cartesiane dell’
atomo i; Mi è la massa atomica dell’ atomo i; PM è il peso molecolare del ligando.
- Generazione del file di restart.
Il file di restart ha lo scopo di permettere il riavviamento di BioDock nel caso in cui la simulazione fosse stata interrotta da un evento imprevisto o dal raggiungimento della convergenza
(nel caso ad esempio che si desideri ridurre la soglia di convergenza) o del numero massimo di
frames specificati dall’utente.
Tale file viene creato all’inizio della simulazione e viene aggiornato al verificarsi di una serie di
situazioni quali il cambiamento del frame di riferimento, l’ottenimento di un numero di frames
utili multiplo di 100 e la sospensione esterna dell’esecuzione del programma decisa
dall’operatore.
53
- Posizionamento ligando e calcolo energetico.
Per ciascun frame vengono generati sei numeri casuali, compresi negli intervalli specificati
dall’utente, corrispondenti alle traslazioni e alle rotazioni rispetto al frame di riferimento e utilizzati per la rototraslazione del ligando. Con gli accorgimenti citati in una sezione precedente di
questa tesi, viene vagliato il frame e se accettabile viene registrato nel file di traiettoria. Sempre
in questa parte del programma vengono applicati gli algoritmi di convergenza, di “Boltzmann
jump” e di interruzione della simulazione secondo i criteri descritti in precedenza , salvo diversa indicazione da parte dell’utente nel file di input.
- Memorizzazione dei dati su supporto magnetico.
Questa sezione viene eseguita in parte contemporaneamente alle altre ed in parte solo al termine della simulazione. In particolare, in ogni istante viene aggiornato un file chiamato di “output”
in cui è riportata, in maniera leggibile da parte dell’utente, l’operazione in atto in quel preciso
istante.
In base agli eventi citati in precedenza, viene aggiornato continuamente il file di restart ed inoltre, quando è individuato un frame utile, i dati vengono registrati nel file di traiettoria. Al termine della simulazione, l’utente può imporre la registrazione del frame migliore direttamente in
forma di coordinate cartesiane.
54
2.3.4 Sezione “Analyzer”
Rappresenta quella parte di BioDock destinata all’analisi geometrica dei frames generati dalla
sezione di docking e all’eliminazione delle posizioni del ligando spazialmente simili.
Anche in questo caso le fasi dell’elaborazione possono essere distinte in:
- Lettura delle coordinate atomiche.
Valgono le stesse considerazioni fatte precedentemente per il modulo di docking, fatta eccezione per il fatto che viene considerato il solo ligando e non si ha alcun assegnazione dei parametri di force field e dei raggi di Van der Waals.
- Rototraslazione dei frames.
La geometria di ogni frame è memorizzata nel file di trajectory in forma binaria mediante i sei
valori corrispondenti alla traslazione e alla rotazione relative.
- Analisi geometrica.
Secondo la procedura sopra descritta, si calcolano le coordinate cartesiane di ciascun frame, e
si memorizzano in forma compatta trasformandole nelle pseudocoordinate Xf, Yf, Zf. Dalla
comparazione dei valori Xf, Yf e Zf di coppie di frames, viene costruita la matrice binaria che
identifica i frames geometricamente simili, e genera così una serie di clusters. I problemi concernenti questa fase sono già stati considerati nella sezione 2.2 della parte sperimentale di questa tesi.
55
- Memorizzazione della nuova traiettoria.
Alla fine dell’analisi, i frames di minima energia che caratterizzano ciascun cluster, vengono registrati in un nuovo file di traiettoria. La loro geometria è sempre espressa in termini di rotazione e traslazione relative rispetto alla posizione iniziale del ligando al fine di permettere una successiva manipolazione.
56
2.4 Programmi aggiuntivi.
BioDock genera i files di trajectory in formato binario (come descritto in precedenza) allo
scopo di immagazzinare il massimo numero di informazioni occupando il minimo spazio possibile sul supporto magnetico. Tuttavia, questo file non è visualizzabile direttamente da parte
dell’utente, sia per quanto riguarda i valori di rototraslazione di ciascun frame, sia per quanto
riguarda la struttura tridimensionale del ligando. Per questo motivo, sono stati creati i programmi ausiliari BtjView ed Extractor, i quali hanno rispettivamente la funzione di visualizzare i
valori α, β, γ, ∆x, ∆y, ∆z di ciascun frame e di produrre un file di coordinate cartesiane in
formato Biosym. Entrambi i programmi mettono a disposizione dell’utente la possibiltà di attivare selezioni multiple per rendere meno tedioso il lavoro di visualizzazione dei dati o di estrazione delle strutture dal file di traiettoria. Un esempio di trajctory resa leggibile da BtjView è
qui di seguito riportato:
Parametri assegnati:
~~~~~~~~~~~~~~~~~~~~
* Frames generati
:
* Ampiezza traslazione X:
* Ampiezza traslazione Y:
* Ampiezza traslazione Z:
* Centro rotazione X
:
* Centro rotazione Y
:
* Centro rotazione Z
:
* Overlap
:
* Cutoff
:
* Boltzmann Jump
:
42
30.0000
30.0000
30.0000
4.3726
70.2612
64.8854
0.1000
10.0000
Off
*
*
*
*
*
*
*
*
*
*
Versione trajectory
Ampiezza rotazione X
Ampiezza rotazione Y
Ampiezza rotazione Z
Assegnaz. centro rot.
Force check
Assegnazione autom.
BioCut
Costante dielettrica
Temperatura
:
1.02
: 360.0000
: 360.0000
: 360.0000
:
Off
:
On
:
Off
:
0.0000
:
1.0000
:
300.00
Frame
dX
dY
dZ
Ax
Ay
Az
Energia
========================================================================
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000 -49.1673
1 -9.6084
9.9907
5.1375 135.2189
-5.5189 109.0294 -17.4815
2 -12.4335
8.3784
2.2529
30.6895
69.8486 112.4279 -22.1076
3 -2.1839 12.2715 11.4432 -146.2777
38.4508
12.2491 -17.0867
4
1.2930 12.0070 11.2455 -11.2158
90.9734 140.5500 -17.9451
5
8.9967 10.1087 14.2233 -174.5302 -15.6555 -10.9913 -17.3939
6 -7.1054 -4.1132 13.3671
45.4772 -96.2285
-9.8499 -17.9763
57
7 -14.8060 -7.8680
8 -9.7164 10.2694
9 -7.6398 13.6298
10 -9.8928 -12.1042
11 -11.7675
5.7576
12 11.9433 14.6684
13 -13.2786 -9.0821
14
4.3127 14.9265
15 -11.5032 -1.0504
16 -13.6509 14.8600
17 -14.3728
9.4368
18 -1.6814 13.0577
19 -13.4946 14.8026
20 14.9711
8.7320
21 -1.0386 -14.7819
22 -1.1039 -7.1733
23 -0.6057
8.6807
24
9.5591 14.9407
25 14.6852 13.2575
26 -1.1684 13.2394
27 -4.7589 14.7903
28 -3.3733 -11.2427
...
14.9090
7.7967 -107.5647 -66.8206
1.2819
26.0660
83.9433 -61.5683
-1.8263 143.5696 166.8646
3.6894
11.8859 -133.9381 -93.2099 -47.9761
6.3808 -164.7298
3.9209 -59.4269
12.5527
47.4820 -150.6435 -65.5363
11.3931 -36.2546
75.2523 105.3169
14.7411
58.4103 -148.2335 -161.7120
14.6781
46.4945 -115.0212 -47.0665
14.1899
72.3756 -136.8579
98.5309
14.6894
44.8872 -101.6509
29.6334
14.8646
68.5401
61.1240 131.1476
1.9689
91.8059 -86.5172 -20.2886
9.3103 -135.9136 -137.4513 -23.1183
14.4537
10.0204 -29.5449 -165.3069
14.8766
35.9246
94.3144
86.8898
14.3003 137.6761 -78.2248
90.0376
6.8164
13.3329
58.6364 110.3728
3.5901 147.5473 -169.5997 -104.7204
1.9277
-6.5388 109.5520 -151.5140
-1.2908 -43.6929 -139.8753 -53.7135
8.2846 -42.8284 -144.5399 -75.1678
Minimo energetico (escluso frame di partenza):
28 -3.3733 -11.2427
8.2846 -42.8284 -144.5399
-75.1678
-3.5627
-19.7604
-18.1580
-18.7828
-23.6655
-3.0594
-19.3833
-2.7831
-12.1058
-0.3758
-2.8385
-5.6646
-4.1748
-11.6675
-11.7727
-15.6489
-21.1087
-13.6714
-17.1436
-18.9089
-14.5288
-27.7674
-27.7674
All’inizio della tabella sono riportati i parametri utilizzati nella simulazione (vedi oltre). Procedendo da sinistra verso destra, ciascuna colonna rappresenta: il numero del frame, le traslazioni rispetto agli assi X, Y, Z in Ångstrom, le rotazioni in gradi intorno gli assi X, Y, Z (corrispondenti agli angoli α, β, γ che compaiono nelle equazioni della sezione 2.1 della parte sperimentale di questa stessa tesi), l’energia di interazione fra ligando e proteina in Kcal/mol. Al
termine di questa stessa tabella è evidenziato il frame con la più bassa energia di interazione.
58
3. Applicazione del metodo BioDock
3.1 L’acetilcolinesterasi
I farmaci appartenenti alla classe degli stimolanti dei recettori colinergici comprendono sia
quelli ad azione diretta, sia quelli ad azione indiretta. I composti colinomimetici ad azione diretta, in base alle caratteristiche di selettività nei confronti dei recettori colinergici, possono
avere una prevalente attività muscarinica oppure nicotinica. I colinomimetici ad azione indiretta
inibiscono il processo d'idrolisi dell'acetilcolina eseguito dall'enzima acetilcolinesterasi che,
scindendola in acido acetico e colina, porta alla sua inattivazione. Riducendo l'entità del fenomeno di inattivazione, l'acetilcolina persiste negli spazi intersinaptici per un periodo più lungo
rispetto alla norma, esercitando così un'azione protratta sui recettori con conseguente potenziamento dell'attività.
L'acetilcolinesterasi (Fig. 10) è un enzima omodimerico caratterizzato dalla presenza di un sito
di legame per la testa cationica dell'acetilcolina e di un sito catalitico dal quale viene catalizzata
l'idrolisi. Quest'ultimo è composto dalla cosiddetta triade catalitica (his440, ser200, glu327) che,
prima dell'interazione col substrato, vede la formazione di un legame d'idrogeno fra un'azoto
imidazolico dell'istidina e l'idrogeno della serina. Tale interazione polarizza il legame idrogeno ossigeno (Fig. 4) della stessa serina al fine di rendere l'atomo d'ossigeno più suscettibile a formare un legame con il carbonio carbonilico dell'acetilcolina e a cedere l'idrogeno. Con un
meccanismo concertato della durata approssimativa di 150 microsecondi [79], l'idrogeno della
serina viene ceduto all'ossigeno della colina con contemporanea rottura del legame estereo co-
59
lina - acetato e formazione del legame acetato - serina. Successivamente, interviene l'acqua
che, grazie alle sue prorietà nucleofile, idrolizza il legame acetato - serina, riportando l'enzima
alla situazione di partenza e liberando una molecola di acido acetico. Il rilascio della colina dal
sito di legame avviene per riduzione dell'affinità con l'enzima e per spiazzamento da parte di
un'altra molecola di substrato.
I composti che sono dotati di attività inibente nei confronti della acetilcolinesterasi possono essere distinti in reversibili (edrofonio, carbammati) (Fig. 8) e in irreversibili (organofosfati). In
entrambi i casi, le molecole vanno ad occupare il sito attivo dell'enzima o una regione ad esso
adiacente, impedendo l'accesso al substrato naturale o la catalisi.
Fig. 4 - Rappresentazione schematica del sito attivo dell’acetilcolinesterasi
60
Fig. 5 - Struttura stereo degli aminoacidi del sito attivo dell’acetilcolinesterasi
Qualora l'interazione inibitore - enzima risulti reversibile, quindi non covalente, l'acetilcolina
pùo competere col farmaco eventualmente spiazzandolo se la sua costante di formazione del
complesso e la sua concentrazione lo consentono. Da ciò è evidente che, in questo caso, l'inibizione enzimatica non può essere completa ed è destinata a terminare dopo un breve periodo.
Qualora, invece, si abbia a che fare con un inibitore irreversibile, quest'ultimo forma legami covalenti con l'enzima rendendo vana qualsiasi forma di competizione da parte dell'acetilcolina. In
quest'ultimo caso, l'inibizione persiste per un periodo prolungato che è funzione del turn-over
dell'acetilcolinesterasi. I composti che sfruttano un tale meccanismo d'azione sono altamente
tossici (Sarin, DFP, Tabun, Soman, ecc). Benché agiscano con lo stesso meccanismo, gli insetticidi organofosforici come il parathion e il malathion sono meno tossici per l'uomo perché
61
inattivi nella forma nella quale vengono utilizzati. Alla base dell'azione di queste due molecole vi
è dapprima una loro attivazione per desolfurazione (sono tiofosfati) operata dal metabolismo
dell'organismo (Fig. 9) che li ha assunti, quindi la formazione del legame fosfoestereo con la
serina della triade catalitica. Quest'ultimo, essendo più stabile rispetto a un normale legame
estereo, non viene idrolizzato in tempi utili e pertanto il sito catalitico viene occupato in modo
irreversibile. Nel caso d'intossicazione da organofosforici, l'unico sistema per ripristinare la
funzionalità dell'acetilcolinesterasi è di somministrare un'ossima (pralidossima, propralidossima,
NIMA, trimedossima, obidossima, ecc) che, nel caso in cui il legame fosfoestereo non sia
troppo "invecchiato" (Fig. 6), può spostare la molecola che occupa l'ossidrile serinico, sfruttando una reazione di transesterificazione. I composti organofosforici trovano applicazione
come insetticidi e come armi chimiche, ma trovano rari impieghi nel settore medico (ecotiopato nella cura del glaucoma).
Il meccanismo d'azione degli inibitori reversibili può essere di due tipi: gli alcoli ammonici quaternari (edrofonio o tensilon) si legano reversibilmente al sito attivo o nei pressi di esso, mediante interazioni elettrostatiche che, essendo labili, permettono al complesso di sopravvivere
per un breve periodo (da 2 a 10 minuti), mentre i carbammati (fisostigmina o eserina, piridostigmina, neostigmina) carbamilano il sito esterasico dell'enzima. I carbammati, essendo più resistenti degli esteri, vengono idrolizzati lentamente e ciò consente di prolungare l'occupazione
del sito per un periodo compreso fra i 30 minuti e le 6 ore.
62
+N
CH3
CH N OH
N
CH N OH
CH3
Pralidossima
Pro-pralidossima
O
CH3 C CH N OH
NIMA
HO N CH
+
+
N CH2CH2CH2 N
CH N OH
Trimedossima
HO N CH
+
+
N CH2OCH2 N
CH N OH
Obidossima
Fig. 6 - Riattivatori dell’acetilcolinesterasi
63
Gli impieghi clinici degli inibitori dell'acetilcolinesterasi sono svariati e comprendono la cura del
glaucoma (fisostigmina, demecario, ecotiopato), dell'ileo paralitico post-operatorio, del megacolon congenito, dell'esofagite da riflusso, della ritenzione urinaria (come conseguenza di intervento chirugico, parto, lesione traumatica, malattia del midollo spinale), della miastenia grave,
di tachiaritmie sovraventricolari e delle intossicazioni da farmaci antimuscarinici. Recentemente
è stata evidenziata la possibilità applicare gli inibitori dell’acetilcolinesterasi, in particolare la tacrina, nella cura del morbo di Alzheimer [80], benchè già nel 1986 [81] ne siano stati documentati gli effetti farmacologici e la tossicologia in questo campo d’applicazione.
Grazie alla risoluzione della struttura cristallografica dell’acetilcolinesterasi, è stato possibile evidenziare un forte momento dipolare nella struttura dell’enzima [82] allineato con il canale imbutiforme (“gorge”) che pone in comunicazione il sito attivo con l’esterno. L’ingresso del substrato nella cavità è favorito dalla forza elettrostatica generata dal dipolo. All’interno del
“gorge”, le catene laterali dei residui aromatici schermano il ligando dall’interazione diretta con
residui anionici che generano il dipolo. L’affinità dell’ammonio quaternario dell’acetilcolina per i
residui aromatici e la forza elettrostatica generata dal dipolo agiscono in concerto per creare
un legame efficiente e selettivo fra substrato e acetilcolinesterasi.
Recentemente Harel et al. [83], hanno studiato il legame di una serie di inibitori
dell’acetilcolinesterasi cocristallizzati con l’enzima. Nel complesso con l’edrofonio (Fig. 7/B)
(non si tratta di un inibitore specifico dell’acetilcolinesterasi), il ligando interagisce con l’indolo
del Trp84 tramite l’azoto quaternario e con la Ser220 e la His440 tramite un doppio legame
d’idrogeno che coinvolge l’ossidrile in meta. La tacrina (Fig.7/C), invece, nel suo complesso
con l’enzima, interagisce con l’indolo del Trp84 sfruttando l’elevata densità elettronica del suo
sistema acridinico. Il decametonio (Fig. 7/A), la cui struttura si estende lungo il canale
d’accesso al sito attivo, utilizza i due azoti quaternari per l’interazione con l’indolo del Trp84 e
64
del Trp279. L’unica differenza importante riscontrata nella struttura enzimatica risiede
nell’orientamento dell’anello benzenico della Phe330, che nel caso del complesso col decametonio è disposto parallelamente alla superficie del “gorge”, mentre nel caso di quello con la tacrina e con l’edrofonio è posizionato in maniera da avere un interazione col ligando. In tutti i
casi presi in considerazione si evidenzia il ruolo importante svolto dai residui aromatici
nell’interazione di ligandi ammonici quaternari con l’enzima acetilcolinesterasi, mettendo in risalto che i residui Trp84 e Phe330 svolgono la funzione di “sottosito anionico” del sito attivo e il
residuo Trp279 rappresenta il sito anionico periferico.
65
Fig. 7
Complessi fra acetilcolinesterasi e decametonio (A), edrofonio (B), tacrina (C)
66
CH3
+ CH3
CH2 CH2 N CH3
CH3
C O
O
Acetilcolina
H3C
O
O
+ CH3
N CH3
CH3
N C O
H3C
H3C NH C O
Carbarile
Neostigmina
O
O
CH3
H3C NH C O
O C
CH3
N
CH3
N
N
CH3
CH3
+N
CH3
Fisostigmina
(eserina)
Piridostigmina
N
HO
+ NH3
+ CH3
N C2H5
CH3
Edrofonio
Tacrina
Cl
Cl
C2H5
O O
C2H5
+
+
CH2 N CH2 CH2 NH C C NH CH2 CH2 N CH2
C2H5
C2H5
Ambenonio
Fig. 8 - Inibitori reversibili dell’acetilcolinesterasi
67
CH3
H3C
H3C
CH3
CH O O
P
CH O
(H3C)3C CH O O
P
H3C
F
F
CH3
Isofluorpato
H5C2
H5C2
Soman
O O
P
O
S
+ CH3
CH2 CH2 N CH3
CH3
Ecotiopato
S
O2N
O
P
O C2H5
metab.
O C2H5
O
O2N
O
Parathion
H3C O
H3C O
S
P
O C2H5
O C2H5
Paraoxon
O
S
P
metab.
CH C
O C2H5
CH2 C
O C2H5
H3C O
H3C O
O
Malathion
O
P
O
S
CH C O C2H5
CH2 C O C2H5
O
Malaoxon
Fig. 9 - Inibitori irreversibili dell’acetilcolinesterasi
68
Fig. 10 - Struttura secondaria dell’acetilcolinesterasi
3.2 Docking del complesso tacrina - acetilcolinesterasi.
Per collaudare le potenzialità del programma BioDock, è stato preso in considerazione
l’enzima acetilcolinesterasi della Torpedo californica come struttura macromolecolare, e
l’inibitore tacrina come ligando da posizionare. Il file di coordinate spaziali utilizzato, reperibile
69
presso l’archivio Internet del Protein Data Bank [53] (file pdb1acj.ent), deriva dall’analisi ai
raggi X della struttura ottenuta dalla cocristallizzazione dellla proteina con l’inibitore [80]. Ciò
significa che si è volutamente partiti da una situazione in cui è nota sperimentalmente la posizione dell’inibitore nel complesso con l’enzima. Dato che i files presenti nell’archivio PDB [53]
sono di origine cristallografica, e pertanto non riportano le coordinate degli atomi di idrogeno,
si è proceduto alla loro modellazione utilizzando il software Insight II (modulo Builder). Sempre con tale programma (modulo Biopolymer), sono stati aggiunti alcuni segmenti di proteina
mancanti (da Asp1 a His3, da Ala486 a Ala489 , da Glu536 a Thr537) e sono state riparate le porzioni presumibilmente non corrette a causa della relativamente bassa risoluzione. Alla struttura
molecolare risultante sono state assegnate le cariche e i parametri di potenziale tramite il modulo Builder. Il campo di forze usato è stato il CVFF [73]. Le parti introdotte ex novo sono
state ottimizzate utilizzando il programma Discover e mantenendo fisse le restanti parti: la procedura di minimizzazione ha incluso 500 iterazioni sfruttando l’algoritmo di minimizzazione
Steepest descent (della discendente più ripida) trascurando le interazioni non-bonded oltre i
20 Å (cutoff ), poi 12500 iterazioni utilizzando l’algoritmo Conjugate gradients.
L’intera struttura, compreso il ligando, è stata poi ottimizata mediante 10500 iterazioni (Conjugate gradients e cutoff 20 Å) lasciando libere le catene laterali degli aminoacidi e mantenendo
fisso lo scheletro principale della proteina (backbone). Un’ulteriore raffinamento (20000 iterazioni, Conjugate gradients, cutoff 10 Å, con limite di convergenza pari a 0.001) è stato ottenuto mantenendo fissa l’intera struttura della tacrina e della acetilcolinesterasi ad esclusione dei
residui aminoacidici (backbone incluso) compresi in una sfera di raggio 15 Å centrata sul centro geometrico della tacrina (X = 4.025; Y = 69.943; Z = 63.296). La selezione dei residui da
muovere col criterio sopra citato è stata possibile grazie al programma BioLock appositamente creato per questa situazione.
70
Ovviamente, l’operazione di rendere mobile (“calda”) una zona interna di una proteina, al sciando “congelato” il resto della struttura vincolata ai dati sperimentali, non ha significato fisico e non sarebbe lecita. In particolare, il rischio è quello di creare artificialmente regioni a
gradiente d’energia relativamente elevato nella fascia di confine tre le zone calda e fredda della
proteina. Tuttavia, il raggio della sfera è sufficientemente ampio per considerare trascurabile, ai
fini del presente lavoro, un’eventuale influenza dei gradienti della zona di confine sul potenziale
del sito di legame dell’inibitore.
Il complesso così ottenuto è stato trattato con BioDock sfruttando il seguente file di input:
Project
1acj_min5
Module
Residue
Pri
Frames
Restart
Docking
999
19
20000
Off
Delta_X
Delta_Y
Delta_Z
1.0
1.0
1.0
Rot_X
Rot_Y
Rot_Z
5
5
5
AutoPar
Off
SaveMin
Off
Cutoff
10
OverLap
0.1
IntFrames
10
ConvMode
0
Confactor
0.95
BoltzJmp
Off
Envar
0.01
Dielectric
1
ForceChk
On
Ciascun comando è caratterizzato da una parola chiave seguita da un parametro. I comandi
attualmente implementati comprendono: AutoPar, BioCut, BoltzJmp, Cent_X, Cent_Y,
71
Cent_Z, Confactor, ConvMode, Cutoff, Delta_X, Delta_Y, Delta_Z, Dielectric, Envar, ForceChk, Frames, IntFrames, Module, Overlap, Pri, Project, Residue, Restart, Rot_X, Rot_X,
Rot_Z, Shell, Temp, Text.
Considerando i comandi nel file di input utilizzato per l’ottimizzazione, Project informa BioDock sul nome del progetto e quindi sul nome dei files che dovrà utilizzare e/o generare, Module seleziona il tipo di modulo (in questo caso quello di docking), Residue indica il numero
del residuo da muovere (ligando) presente nel file di coordinate, Pri inpone l’esecuzione del
programma con la priorità specificata (comando attivo solo sulla versione di BioDock per Silicon Graphics), Frames sono il numero di frames utili da generare, Restart attiva o meno la
modalità di restart (come descritto in precedenza), Delta_X, Delta_Y, Delta_Z, Rot_X,
Rot_X, Rot_Z definiscono gli intervalli delle rototraslazioni ammesse, AutoPar,se On, attiva
l’assegnazione automatica degli intervalli di rototraslazione, SaveMin, se On, impone la registrazione su disco del frame di minimo in forma di coordinate cartesianeoltre che nel normale
formato binario, CutOff definisce il valore in Ångstrom della distanza di cutoff, OverLap indica il valore di sovrapposizione dei raggi di Van der Waals oltre il quale la repulsione è da considerare eccessiva, IntFrames rappresenta il numero di frames utili dopo il quale applicare
l’algoritmo di convergenza, ConvMode è il tipo di algoritmo di convergenza (0 lineare, 1
esponenziale), Confactor è la pendenza con la quale avviene la convergenza in caso di applicazione dell’algoritmo a funzione lineare, BoltzJmp, se On, attiva la modalità Boltzmann jump,
Envar è il valore della deviazione standard calcolato sull’ultima serie di frames utili che consente a BioDock di stabilire il raggiungimento della convergenza (come citato in precedenza in
questa tesi), Dielectric è il valore della costante dielettrica, ed infine ForceChk, se On, forza il
superamento del controllo della sovrapposizione dei raggi di Van der Waals alla struttura di
partenza. L’output generato da questa simulazione è il seguente:
72
**********************
**** BioDock V1.1 ****
**********************
ATTENZIONE: BumpCheck sulla struttura di partenza non eseguito.
Ottimizzazione frame 1acj_plus_159406
**** Docking ****
Parametri assegnati:
~~~~~~~~~~~~~~~~~~~~
* Frames da generare
:
* Ampiezza traslazione X:
* Ampiezza traslazione Y:
* Ampiezza traslazione Z:
* Ampiezza rotazione X :
* Ampiezza rotazione Y :
* Ampiezza rotazione Z :
* Overlap
:
* Cutoff
:
* Costante dielettrica :
* Boltzmann Jump
:
10000
1.0000
1.0000
1.0000
5.0000
5.0000
5.0000
0.5000
10.0000
1.0000
Off
Caratteristiche ligando:
~~~~~~~~~~~~~~~~~~~~~~~~
Numero residuo:
Dimensione:
Peso molecolare:
Carica:
Coord. centro massa (X,Y,Z):
999
30 atomi
199.1235 Daltons
1.0002
4.372565 70.261192 64.885353
Caratteristiche proteina:
~~~~~~~~~~~~~~~~~~~~~~~~~
Dimensione reale:
Peso molecolare:
Carica:
8619 atomi
61712.4570 Daltons
-9.9999
Modalita' di convergenza 0 attivata.
* Numero interframes
:
10
* Fattore di convergenza: 0.950000
Nuovo frame di riferimento:
Nuovo frame di riferimento:
Nuovo frame di riferimento:
Continuo per altri 10 frames
Nuovo frame di riferimento:
...
Continuo per altri 10 frames
Continuo per altri 10 frames
Continuo per altri 10 frames
Continuo per altri 10 frames
Continuo per altri 10 frames
Continuo per altri 10 frames
Nuovo frame di riferimento:
Continuo per altri 10 frames
8
energia
-37.4559
11
energia
-42.2721
26
energia
-46.7950
riferendomi al frame 26
47
energia
-47.5598
riferendomi al frame 683
riferendomi al frame 683
riferendomi al frame 683
riferendomi al frame 683
riferendomi al frame 683
riferendomi al frame 683
789
energia
-49.7090
riferendomi al frame 789
73
File "1acj_THA.cor" salvato.
Minimo energetico: N.
Frames utili:
Numero collisioni:
Cicli utili:
Velocita' (Iter./min)
Lunghezza trajectory:
789 energia
-49.7090 Kcal/mol
810
701
53.61 %
21.26
22788 bytes
Tempo impiegato: 1 ora 11 minuti 5 secondi.
Successivamente la molecola della tacrina è stata isolata, protonata ed ottimizzata col programma quantomeccanico semiempirico MOPAC (con l’algoritmo AM1) al fine di ottenere
una distribuzione di carica migliore di quella empirica assegnata dal software Biosym [73], soprattutto per quanto riguarda il gruppo ammonico. Il ligando ottimizzato è stato reinserito
nell’enzima mantenendo la precedente posizione e il complesso così ottenuto è stato ottimizzato con BioDock. Rispetto al docking precedente, è stato eseguita una simulazione completamente casuale, senza cioè utilizzare gli algoritmi di convergenza. In questa maniera sono stati
ottenuti 600.000 frames tramite tre successive simulazioni da 200.000.
Il file di input utilizzato è stato il seguente:
Project
Module
Residue
Pri
Frames
Restart
1acj_plus
Docking
999
19
200000
On
Delta_X
Delta_Y
Delta_Z
30
30
30
Rot_X
Rot_Y
Rot_Z
360
360
360
AutoPar
Off
74
SaveMin
Off
Cutoff
10
OverLap
0.1
IntFrames
0
BoltzJmp
Off
Dielectric
1
ForceChk
On
I parametri rototraslazionali conferiscono al ligando la massima libertà di movimento, mentre
IntFrames è posto a zero per evitare che si converga su un frame specifico. In questa situazione, l’output generato è stato il seguente:
**********************
**** BioDock V1.1 ****
**********************
ATTENZIONE: BumpCheck sulla struttura di partenza non eseguito.
**** Docking ****
...
Minimo energetico: N.159406 energia
Frames utili:
Numero collisioni:
Cicli utili:
Velocita' (Iter./min)
Lunghezza trajectory:
-22.2366 Kcal/mol
200000
7168626
2.71 %
2102.96
5600108 bytes
Tempo impiegato: 2 giorni 10 ore 23 minuti 56 secondi.
In questo modo sono sati creati dei frames con lo scopo di trovare una serie di situazioni spaziali potenzialmente interessanti per l’interazione del ligando con la proteina, dalle quali partire
per successive ottimizzazioni. Il file di traiettoria risultante è stato sottoposto all’analisi geometrica tramite il modulo Analyzer di BioDock. Un tipico file di input per la comparazione dei
frames è il seguente:
75
Project
1acj_plus
Module
Dist
Analyzer
1
Il parametro Dist, non disponibile per il modulo Docking, rappresenta la distanza media fra gli
atomi al di sotto della quale i frames comparati si considerano simili (vedi sezione 2.2 della
parte sperimentale). Sono stati utilizzati tre valori di Dist, rispettivamente 1, 2 e 3 Å, che hanno
consentito di ottenere per ciascuna simulazione 755, 122 e 42 clusters. Questo fatto evidenzia
che allargando il criterio di similitudine, diminuisce il numero di famiglie di frames. Un tipico
esempio di output di analisi geometrica è il seguente:
**********************
**** BioDock V1.1 ****
**********************
**** Analyzer ****
Parametri assegnati:
~~~~~~~~~~~~~~~~~~~~
* Frames da analizzare
* Iterazioni previste
* Termine di paragone
* Numero ligando
* Atomi ligando
:
:
:
:
:
400000
8E10
3.0000
999
30
FASE 1: Rototraslo i frames.
FASE 2: Comparo i frames.
Lavoro Eseguito:
0
25
50
75
100
##################################################
FASE 3: Genero la nuova trajectory.
Numero clusters:
42
Tempo impiegato: 5 minuti 56 secondi.
76
Dalla traiettoria ottenuta con Dist pari ad un Ångstrom, è stato estratto il frame di minima
energia il quale rappresenta il cluster con i frame di miglior energia. I parametri di rototraslazione che caratterizzano questo frame sono:
Frame
dX
dY
dZ
Ax
Ay
Az
Energia
======================================================================
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000 -49.1673
Minimo energetico (escluso frame di partenza):
550 -3.3733 -11.2427
8.2846 -42.8284 -144.5399
-75.1678
-27.7674
Il frame 0 rappresenta la situazione di partenza, mentre il frame 550 corrisponde alla nuova
situazione spaziale individuata che successivamente è stata ottimizzata utilizzando il modulo di
docking di BioDock con il seguente file di input:
Project
1acj_min2
Module
Residue
Pri
Frames
Restart
Docking
999
17
20000
Off
Delta_X
Delta_Y
Delta_Z
1.0
1.0
1.0
Rot_X
5
Rot_Y
Rot_Z
5
5
AutoPar
SaveMin
Cutoff
OverLap
IntFrames
Confactor
ConvMode
BoltzJmp
Envar
Dielectric
ForceChk
BioCut
Text
Off
On
10
0.1
10
0.95
0
Off
0.01
1
Off
20
Ottimizzazione frame 550
77
Degno di nota è il parametro BioCut 20 che consente di accelerare notevolmente
l’ottimizzazione trascurando tutti gli atomi della proteina oltre i 20 Å di distanza dal centro di
rotazione del ligando (come già trattato in questa stessa tesi nella sezione 2.1).
78
4. Risultati e discussione
Il problema del posizionamento di antagonisti all’interno di una biomacromolecola non è di
semplice soluzione. Se infatti nel caso degli agonisti si suppone che essi debbano interagire
con il sito attivo della struttura proteica al fine di scatenare quella serie di reazioni che portano
all’attività biologica, gli antagonisti, o inibitori enzimatici possono:
• interagire anch’essi con il sito attivo, in modo transitorio o legandosi covalentemente;
• interagire con un sito adiacente al sito attivo, impedendo per ragioni conformazionali, elettrostatiche o steriche l’accesso dell’agonista al sito stesso;
• interagire con una porzione di macromolecola molto lontana dal sito attivo, impedendo ad
esempio l’apertura di “flaps” (come avviene per gli inibitori dell’ HIV proteasi) o ostruire
fisicamente l’ingresso di un canale che consente l’accesso dell’agonista al sito attivo.
E’ chiaro pertanto che i tentativi di docking di un antagonista debbano basarsi su una struttura
sperimentale, in quanto, in caso di molteplici possibilità di legame, i force fields attualmente disponibili ben difficilmente sono in grado di valutare con un buon grado di affidabilità quale di
esse sia la più probabile. Un metodo di docking può essere considerato valido se è in grado di
riprodurre, almeno come una tra le strutture probabili, quella ottenibile sperimentalmente.
Una fase successiva può essere quella predittiva, in cui di un ligando noto, correttamente posizionato, viene simulata la modificazione ai fini di migliorare la sua affinità per il sistema. Le
nuove strutture possono quindi essere sintetizzate, sottoposte all’analisi cristallografica e valutate dal punto di vista farmacologico. Tale processo iterativo purtroppo è raramente realizza-
79
bile allo stato attuale delle tecnologie, sia computazionali che strumentali, ma rappresenta un
ciclo ideale a cui tendere per lo sviluppo di farmaci mediante l’approccio diretto.
Il software sviluppato in questa tesi ha riguardato l’aspetto di docking privilegiando la velocità
di calcolo e la semplicità d’utilizzo. La tabella dei tempi di calcolo, riportata a pag. 72, evidenzia come le tre strutture considerate come probabili candidate siano state individuate in tempi
accettabili. In particolare il minimo più stabile è stato individuato rapidamente, come descritto
a pag. 70 (anche se ciò è scarsamente indicativo, in quanto tale minimo coincide con la struttura sperimentale che è stata scelta come punto di partenza).
Complesso
A
B
C
∆E (kcal mol -1)
0.0
14.94
25.51
Tabella 1
Il metodo BioDock ha consentito di individuare tre minimi (Fig. 12) in un intervallo di energia
pari a 25 Kcal/mol . Tale intervallo, piuttosto ampio, è stato scelto in modo da ovviare in parte
ai limiti dovuti all’approssimazione del potenziale ai soli due termini di Van der Waals e coulombiano.
Il minimo largamente più stabile (∆E = 0, vedi Tabella 1), coincide con la posizione sperimentale raffinata mediante il software Biosym. La tacrina è posizionata profondamente nel canale
che attraversa la proteina e in cui viene “aspirato” l’agonista secondo il meccanismo proposto
da Ripoll et al. [82].
L’antagonista si trova in una regione adiacente alla triade catalitica, anche se, dai risultati ottenuti nel presente lavoro, non sembra interagire direttamente con nessuno dei tre aminoacidi in80
teressati (Fig. 11/A e Fig. 13). L’interazione con la proteina è ampiamente stabilizzata da forti
interazioni π - π tra gli anelli insaturi dell’inibitore ed alcuni residui aminoacidici con catene laterali aromatiche (Fig. 14). Infatti si osserva una disposizione parallela ed una distanza di soli
3.70 Å tra il centro geometrico dell’arile della Phe330 e l’anello benzenico della tacrina.
Sull’altro lato della molecola, è osservabile un’analoga interazione tra il benzene del sistema
indolico del Trp84 e la sottostruttura piridinica della tacrina (distanza tra i centri geometrici 3.76
Å). L’interazione è ulteriormente rinforzata dal legame d’idrogeno di un protone della testa
ammonica con una delle molecole d’acqua di cristallizzazione presenti nella cavità in prossimità
del sito attivo. La stessa molecola d’acqua forma poi un secondo legame d’idrogeno con una
seconda molecola di solvente, che era già stato evidenziato nella struttura cristallografica [83].
81
Fig. 11
Non è stato invece osservato, a causa della diversa disposizione dell’ His440, il legame
d’idrogeno sperimentalmente trovato nel complesso cristallino, tra NH della His440 e l’azoto
del nucleo triciclico della tacrina. La disposizione dell’arile della Tyr442 rispetto alla tacrina
stessa evidenzia un’interazione, in questo caso perpendicolare, tra i due sistemi aromatici.
E’ interessante osservare come BioDock sia in grado di ipotizzare altre due posizioni di interazione alternative plausibili, anche se con energia d’interazione superiore, quindi meno probabili. Nel secondo minimo individuato (Fig. 11/B e Fig. 15), la tacrina risulta posizionata a qualche Ångstrom di distanza lungo il “canale d’accesso” al sito attivo. In questo caso il complesso
82
è stabilizzato dall’interazione π - π con il Trp279 (3.84 Å) tra il centroide del benzene indolico e
centroide dell’anello piridinico, da un legame d’idrogeno tra un protone ammonico e l’ossigeno
carbonilico della Tyr334 e da un debole legame d’idrogeno tra l’azoto del nucleo triciclico (raramente osservato nel caso di questo studio) e una molecola d’acqua di cristallizzazione. E’
interessante notare come in questo caso la funzione ammonica della tacrina si collochi in una
posizione analoga a quella occupata dalla testa cationica del decametonio nel complesso cristallizzato con lo stesso enzima (Fig. 7/A)
Il terzo minimo considerato (Fig. 11/C e Fig. 16), ulteriormente meno stabile di circa 10
Kcal/mol, si trova in una zona ancora più periferica, in pratica all’ingresso del canale imbutiforme dell’acetilcolinesterasi (Fig. 12). In questo caso, oltre a generiche interazioni di carattere
idrofobico non facilmente identificabili, si può osservare la formazione di un legame d’idrogeno
tra l’azoto del nucleo triciclico della tacrina e l’ossidrile della Tyr70.
L’andamento energetico nel confronto tra i tre complessi simulati A, B, C è ampiamente giustificato dal numero e dal tipo di interazioni favorevoli che la tacrina riesce ad instaurare con la
proteina. Guardando poi simultaneamente i tre minimi in rapporto alla struttura proteica, si ha
quasi l’impressione di un “cammino” percorso dall’inibitore dalla posizione C alla A lungo il
canale d’accesso al sito attivo in evoluzione verso la situazione più stabile.
83
5. Conclusioni
E’ evidente che le conclusioni del lavoro illustrato in questa tesi debbano essere considerate
ampiamente preliminari. Infatti, una sola applicazione è assolutamente insufficiente per valutare
appieno le potenzialità positive e/o le carenze del metodo BioDock. Inoltre, l’impiego di un
campo di forze limitato come quello impiegato in questa tesi fornisce un significato solo qualitativo ad esempio alle valutazioni energetiche dei complessi simulati. Sono in corso ricerche
per apportare diverse migliorie in questo senso alla metodica, ad esempio la possibilità di introdurre un termine esplicito per i legami d’idrogeno, o un termine che esprima il contributo
dell’energia di desolvatazione all’energia globale del sistema. Sarebbe opportuno anche introdurre una costante dielettrica dipendente dalla distanza al posto del dielettrico continuo standard adottato in quasta tesi.
Tuttavia sarà anche da valutare quanto l’introduzione di questi fattori sia vantaggiosa dal punto
di vista del rapporto tra la qualità dei risultati e l’efficienza di calcolo.
Sempre in rapporto al tipo di potenziale utilizzato, è allo studio l’idea di fornire all’utente la
possibilità di applicare campi di forze differenti dal CVFF; inoltre la possibilità da parte
dell’utente di fornire le coordinate atomiche mediante un file PDB anzichè in formato Biosym
renderebbe il software di uso più generalizzato.
In ogni caso si può senz’altro affermare che, nel caso del complesso acetilcolinesterasi - tacrina, la metodica BioDock si è rivelata in grado di operare una scansione estremamente dettagliata delle strutture di complessi possibili. Infatti sono stati simulati diversi milioni di strutture ,
da cui sono state estratte più di 600.000 modalità di interazione “utili”. Tra queste, attraverso
l’analisi dei clusters, l’eliminazione delle forme ridondanti e la successiva ottimizzazione, sono
84
state individuate come probabili le tre strutture A, B, e C (vedi sopra). Tra queste poi, la
struttura A è stata valutata come sostanzialmente migliore dal punto di vista dell’energia
d’interazione, in ottimo accordo con i dati cristallografici. I risultati sono stati ottenuti con tempi
di calcolo ragionevolmente brevi (la generazione e filtrazione di diversi milioni di strutture,
l’analisi di 600.000 strutture e l’ottimizzazione hanno richiesto circa 6 giorni di calcolo su una
INDIGO 2), e con un impiego efficiente dei supporti magnetici mediante memorizzazione compatta delle informazioni.
Si può pertanto concludere che, almeno nella presente applicazione, BioDock si è rivelato uno
strumento valido, efficiente, economico e di semplice utilizzo, ed è stato in grado di riprodurre
esattamente il dato cristallografico.
85
Fig. 12
Orientamenti A, B, C della tacrina rispetto alla struttura enzimatica
86
Fig. 13
Relazione spaziale fra orientamento A e triade catalitica
87
Fig. 14
Complesso A: interazione fra tacrina e residui aminoacidici dell’acetilcolinesterasi
88
Fig. 15
Complesso B: interazione fra tacrina e residui aminoacidici dell’acetilcolinesterasi
89
Fig. 16
Complesso C: interazione fra tacrina e il residuo Tyr70 dell’enzima
90
Bibliografia
1 - Y. C. Martin, M. G. Bures, P. Willet el al., Rewiews in Comput. Chem. 1,
pp 213-263, VCH, New York, NY (1990)
2 - F. C. Bernstein, T. F. Koetzle et al., J. Mol. Biol. 112, 535 (1977)
3 - W. C. Ripka, J. M. Blaney et al., “Topics in Stereochemistry” 20, pp 1-85, Wiley,
New York, NY (1991)
4 - P. J. Goodford, J. Med. Chem. 28, 849 (1985)
5 - M. Von Itzstein, W. Wu, Y. Kok et al., Nature 363, 418 (1993)
6 - N. Pattabiraman, M. Levitt, T. E. Ferrin, R. Cameron, J. Comput. Chem. 6, 432
(1985)
7 - J. M. Blaney, Ph.D. Thesis, University of California, San Francisco, CA (1982)
8 - C. R. Beddel, P. J. Goodford et al., Br. J. Pharmacol. 57, 201 (1976)
9 - S. H. Rotstein, M. A. Murcko, J. Comput.-Aided Mol. Design 7 ,23 (1993)
10 - S. H. Rotstein, M. A. Murcko, J. Med. Chem. 36, 1700 (1993)
11 - F. H. Allen, S. Bellard et al., Acta Crystallogr. B35, 2331 (1979)
12 - I. D. Kuntz, Science 257, 1078 (1992)
13 - J. T. Bolin, D. J. Filman et al., J. Biol. Chem. 257, 13650 (1982)
14 - P. M. Dean, “Molecular Fundations of Drug-Receptor Interaction”, Cambridge
versity Press, Cambridge (1987)
91
Uni-
15 - P. R. Andrews, M. Tintelnot, “Quantitative Drug Design” 4, pp 321-347, Pergamon
Press, New York, NY (1990)
16 - A. R. Fersht, Trends Biochem. Sci. 9, 145 (1984)
17 - D. Sali, M. Bycroft, A. Fersht, “Techniques in Protein Chemistry” 2, pp 295-303,
Academic Press, New York, NY (1991)
18 - A. R. Fersht, J. Shi, J. Knill-Jones et al., Nature 314, 235 (1985)
19 - J. M. Blaney, P. K. Weimer, A. Dearing, P. A. Kollman, E. C. Jorgensen, S. J.
Oatley, J. M. Burridge, C. C. F. Blake, J. Am. Chem. Soc. 104, 6424 (1982)
20 - G. Vistoli, Tesi di Laurea, Milano (1994)
21 - D. Eisenberg, A. D. McLachlan, Nature 319, 199 (1986)
22 - L. Wesson, D. Eisenberg, Protein Sci. 1, 227 (1992)
23 - T. Ooi, T. Oobatake, G. Nemethy, H. A. Scheraga, Proc. Natl. Acad. Sci. USA 84,
3086 (1987)
24 - C. J. Cramer, D. G. Truhlar, Science 256, 213 (1992)
25 - H. J. Böhm, J. Comput.- Aided Mol. Design 6, 593 (1992)
26 - R. S. BoHacek, C. McMartin, J. Med. Chem. 35, 1671 (1992)
27 - I. D. Kuntz, J. M. Blaney, S. J. Oatley, R. Langridge, T. E. Ferrin, J. Mol. Biol. 161,
269 (1982)
28 - M. L. Connoly, Science 221, 709 (1983)
29 - R. L. DesJarlais, R. P. Sheridan, G. L. Seibel, J. S. Dixon, I. D. Kuntz, R.
Venkataraghavan, J. Med Chem. 31, 722 (1988)
30 - B. K. Shoichet, D. L. Bodian, I. D. Kuntz, J. Comput. Chem. 13, 380 (1992)
31 - E. C. Meng, B. K. Shoichet, I. D. Kuntz, J. Comput. Chem. 13, 505 (1992)
32 - S. J. Weiner, P. A. Kollman, D. A. Case et al., J. Am. Chem. Soc. 106, 765 (1984)
92
33 - E. C. Meng, D. A. Geschwend, J. M. Blaney, I. D. Kuntz, Proteins (1993)
34 - I. A. Vakser, Protein Engeenering 8, 371 (1995)
35 - J. Singh, J. Saldanha, J. M. Thornton, Prot. Eng. 4, 251 (1991)
36 - H. Wang, J. Comp. Chem. 12, 746 (1991)
37 - B. Lee, F. M. Richards, J. Mol. Biol. 55, 379 (1991)
38 - B. K. Shoichet, I. D. Kuntz, Protein Eng. 6, 723 (1993)
39 - R. L. DesJarlais, G. L. Seibel, I.D. Kuntz, P. S. Furth et al., Proc. Natl. Acad. Sci.
USA 87, 6644 (1990)
40 - A. Leach, J. Mol. Biol. (1993)
41 - R. L. DesJarlais, R. L. Sheridan, J. S. Dixon, I. D. Kuntz, R. Venkataraghaven, J.
Med. Chem. 29, 2149 (1986)
42 - J. Desmet, M. Demaeyer, B. Hazes, I. Lasters, Nature 356, 539 (1992)
43 - J. S. Dixon, “Trends in QSAR and Molecular Modelling 92”, pp. 412-413,
ESCOM, Leiden (1993)
44 - D. S. Goodsell, A. J. Olson, Proteins 8, 195 (1990)
45 - N. Metropolis, A.W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, J. Chem. Phys.
21, 1087 (1953)
46 - T. N. Hart, R. J. Read, Proteins 13, 206 (1992)
47 - S.-Y. Yue, Protein Eng. 4, 177 (1990)
48 - A. S. Smeille, G. M. Crippen, W. G. Richards, J. Chem, Inf. Comput. Sci. 31, 386
(1991)
49 - N. Kasinos, G. A. Lilley, N. Subbarao, I. Hanfee, Protein Eng. 5, 69 (1992)
50 - M. Yamada, A. Itai, Chem. Farm. Bull. 41, 1200 (1993)
51 - M. Yamada, A. Itai, Chem. Farm. Bull. 41, 1203 (1993)
93
52 - A. Vedani, J. Comput. Chem. 9, 269 (1988)
53 - Protein Data Bank (PDB), Chemstry Department, Building 555 Brookhaven
National Laboratory Upton, NY (USA)
54 - I. A. Vakser, C. Aflalo, Proteins 20, 320 (1994)
55 - M. C. Lawrence, P. C. Davis, Proteins 12, 31 (1992)
56 - D. J. Bacon, J. Moult, J. Mol. Biol. 225, 849 (1992)
57 - M. L. Connolly, Biopolymers 25,1229 (1986)
58 - A. D. McLachlan, J. Mol. Biol. 128, 49 (1979)
59 - D. Fischer, R. Norel, H. Wolfson, R. Nussinov, Proteins 16, 278 (1993)
60 - P. H. Walls, M. J. E. Sternberg, J. Mol. Biol. 228, 27 (1992)
61 - F. Jiang, S.-H. Kim, J. Mol. Biol. 219, 79 (1991)
62 - M. Billeter, T. F. Havel, I. D. Kuntz, Biopolymers 26, 777 (1987)
63 - J. Holland, “Adaptation in Natural and Artificial Systems”, Ann Arbor, University of
Michigan Press (1975)
64 - T. Brodmeier, E. Pretsh, J. Comput. Chem. 15, 588 (1994)
65 - Y. L. Xiao, D. E. Williams, J. Phys. Chem. 98, 7191
66 - Y. Nishibata, A. Itai, Thetraedron 47, 8985 (1991)
67 - G. Del Re, J. Chem Soc., 4031 (1958)
68 - B. W. Kernighan, D. M. Ritchie, “The C Programming Language”, Prentice-Hall,
(1978)
69 - B. W. Kernighan, “ANSI X3.159-1989, the ANSI C Standard. Programming in
C-a tutorial”
70 - D. M. Ritchie, “C Reference Manual”
94
71 - Si ringrazia il Consorzio Interuniversitario di Biotecnologie (CIB) per aver concesso
l’uso del computer.
72 - Si ringrazia il Sig. Carlo Penati per aver concesso l’utilizzo del suo sistema
MS-DOS
73 - Biosym Technologies, 10065 Barnes Canyon Road, San Diego, CA
74 - J. P. Stewart, J. Comput.-Aided Mol. Design 4, 1 (1990)
75 - GNU - Free Software Fundation Inc., 675 Mass Ave, Cambridge MA (USA)
76 - Copyright 1992/94 Michael Sinz, 1992/93 Commodore - AMIGA Inc.
77 - Borland International Inc., 1800 Green Hills Road, Scotts Valley CA (USA)
78 - J. A. McCammon et al., J. Comput. Chem. 16, 454 (1995)
79 - B. G. Katzung, “Farmacologia generale e clinica” ed. italiana 1, pp 73-86, Piccin
Nuova Libraria, Padova (1989)
80 - D. Whinker, J. Am. Med. Ass. (JAMA) 6, 1023 (1994)
81 - M. V. Majovski, G. M. Marsh, T. K. Kling, N. Eng. Med. 315, 1241 (1986)
82 - D. R. Ripoll, C.H. Faerman, P. H. Axelsen, I. Silman, Proc. Natl. Acad. Soc. USA
90, 5128 (1993)
83 - M. Harel, I. Shalk, P.H. Axelsen et al., Proc. Natl. Acad. Soc. USA 90, 9031 (1993)
95
Sommario
Stato dell’arte
1. Introduzione al docking
pag. 1
2. Metodi proposti
pag. 4
Parte sperimentale
1. Introduzione al metodo
pag. 35
2. Principi teorici
2.1 Algoritmi di docking
pag. 36
2.2 Algoritmi di analisi
pag. 43
2.3 Descrizione del programma
2.3.1 Struttura di BioDock
pag. 46
2.3.2 Sezione “Main”
pag. 49
2.3.3 Sezione “Docking”
pag. 50
2.3.4 Sezione “Analyzer”
pag. 52
2.4 Programmi aggiuntivi
pag. 54
3. Applicazione del metodo BioDock
3.1 L’acetilcolinesterasi
pag. 56
3.2 Docking del complesso tacrina - acetilcolinesterasi
pag. 66
4. Risultati e discussione
pag. 76
5. Conclusioni
pag. 81
6. Bibliografia
pag. 88
Ringraziamenti
Desidero ringraziare il Prof. Luigi Villa, che ha contribuito a suscitare, sviluppare e discutere le
tematiche affrontate in questa tesi con grande esperienza e spirito innovativo.
Desidero inoltre ringraziare la Dott.ssa Anna Maria Villa per la pazienza e la competenza con
la quale mi ha guidato durante tutto il mio lavoro.
Desidero infine ringraziare il Dott. Giulio Vistoli per i preziosi consigli e il Dott. Patrizio Silvi
Antonini per le fotografie.
96