Studio e sviluppo di metodi computazionali per l’analisi ... · e le tecniche dettagliate...

UNIVERSITA DEGLI STUDI DI PARMA

FACOLTA DI SCIENZE

MATEMATICHE, FISICHE e NATURALI

Corso di Laurea in Informatica

Tesi di Laurea

Studio e sviluppodi metodi computazionali

per l’analisi delle conformazionidi molecole

Relatore: Candidato:

Prof. Alessandro Dal Palu Tommaso Nanu

Correlatore:

Prof. Pietro Cozzini

Anno Accademico 2010/2011

A mannedda mea Mariantoniae a thiu Zoseppe.

Ringraziamenti

Sono passati ben cinque anni da quando ho fatto i bagagli e ho lasciato lamia famiglia, i miei amici e la mia terra di origine. Sono parecchi anni.Se in questo momento sto scrivendo queste poche righe, lo devo tutto ai mieigenitori, Lucia e Tonino. Nei limiti del possibile, mi e stato concesso tutto,non mi e mai stato fatto mancare niente.E grazie anche di aver concepito la mia sorellona, Susanna. Con lei ho passatoquasi tutti i momenti della mia vita, brutti e belli. Grazie Susa!Un sentito ringraziamento va al Professor Alessandro Dal Palu che nonostantemille impegni e riuscito a dedicare numerose ore alla buona riuscita del lavorodi tirocinio e di tesi. Ringrazio anche il Professore Pietro Cozzini per ladisponibilita concessa.Un ringraziamento va a tutti gli amici di Olbia e di Parma che, in tutti questianni, mi hanno sopportato. Un grazie particolare va a Michi: senza di lui, aquest’ora, potrei rotolare solo da una strada in discesa!

i

Indice

Introduzione v

1 Background 11.1 Proteine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 Caratteristiche strutturali . . . . . . . . . . . . . . . . 21.1.2 Sito attivo . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.2 Legame chimico e forze intermolecolari . . . . . . . . . . . . . 71.2.1 Forze intermolecolari . . . . . . . . . . . . . . . . . . . 8

1.3 Strutture dati fondamentali . . . . . . . . . . . . . . . . . . . 81.3.1 Grafo . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.4 Programmazione a vincoli . . . . . . . . . . . . . . . . . . . . 91.4.1 Alberi di ricerca . . . . . . . . . . . . . . . . . . . . . . 11

1.5 Drug Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.6 Algoritmi di Docking . . . . . . . . . . . . . . . . . . . . . . . 14

1.6.1 Algoritmi con corpi rigidi . . . . . . . . . . . . . . . . 141.6.2 Algoritmi con ligando flessibile . . . . . . . . . . . . . . 141.6.3 Docking attraverso simulazione . . . . . . . . . . . . . 151.6.4 Panorama del mercato attuale . . . . . . . . . . . . . . 16

2 Scopo della Tesi 192.1 Ligand-rotation . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.1.1 Analisi dei requisiti . . . . . . . . . . . . . . . . . . . . 21

3 Formalizzazione 233.1 Modello del problema a vincoli . . . . . . . . . . . . . . . . . . 24

3.1.1 Variabili . . . . . . . . . . . . . . . . . . . . . . . . . . 243.1.2 Dominio delle variabili . . . . . . . . . . . . . . . . . . 243.1.3 Vincoli implementati . . . . . . . . . . . . . . . . . . . 25

3.2 Visita dello spazio di ricerca . . . . . . . . . . . . . . . . . . . 253.2.1 Albero di ricerca . . . . . . . . . . . . . . . . . . . . . 26

3.3 Generazione delle conformazioni . . . . . . . . . . . . . . . . . 27

iii

iv INDICE

3.4 Rotazione di un legame . . . . . . . . . . . . . . . . . . . . . . 323.5 Esplorazione del sito attivo della proteina . . . . . . . . . . . . 34

3.5.1 Construzione delle basi ortonormali . . . . . . . . . . . 353.5.2 Generazione ligandi base . . . . . . . . . . . . . . . . . 353.5.3 Campionamento di S . . . . . . . . . . . . . . . . . . . 36

4 Dettagli implementativi 374.1 Inizializzazione di ligand-rotation . . . . . . . . . . . . . . . . 37

4.1.1 Acquisizione informazioni di configurazione . . . . . . . 384.1.2 Struttura file mol2 . . . . . . . . . . . . . . . . . . . . 40

4.2 Strutture dati del risolutore . . . . . . . . . . . . . . . . . . . 414.2.1 Matrice di Vicinanza . . . . . . . . . . . . . . . . . . . 424.2.2 Discretizzazione di un sottoinsieme limitato di R3 . . . 424.2.3 Celle contenenti atomi ‘vicini’ . . . . . . . . . . . . . . 45

4.3 Algoritmo di ligand-rotation . . . . . . . . . . . . . . . . . . . 464.3.1 Parte I: inizializzazione . . . . . . . . . . . . . . . . . . 474.3.2 Parte II: generazione conformazioni . . . . . . . . . . . 48

4.4 Consistenza del vincolo . . . . . . . . . . . . . . . . . . . . . . 504.4.1 Verifica dei vincoli . . . . . . . . . . . . . . . . . . . . 51

4.5 Accenni di Complessita . . . . . . . . . . . . . . . . . . . . . . 51

5 Risultati 535.1 Applicazioni possibili . . . . . . . . . . . . . . . . . . . . . . . 545.2 Efficienza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.2.1 Variazione parametri qualitativi . . . . . . . . . . . . . 575.3 Confronti su visite . . . . . . . . . . . . . . . . . . . . . . . . 60

5.3.1 DFS sul grafo della molecola . . . . . . . . . . . . . . . 605.3.2 Modifica della DFS: first fail . . . . . . . . . . . . . . . 62

6 Conclusioni e Sviluppi futuri 656.1 Molecole contenenti cicloesano . . . . . . . . . . . . . . . . . . 66

6.1.1 Cicloesano . . . . . . . . . . . . . . . . . . . . . . . . . 676.1.2 Nuovi gradi di liberta . . . . . . . . . . . . . . . . . . . 68

6.2 Espansione modello a vincoli . . . . . . . . . . . . . . . . . . . 686.2.1 Controllo lunghezza legame . . . . . . . . . . . . . . . 686.2.2 Implementazione propagazione dei vincoli . . . . . . . 68

Bibliografia 71

Introduzione

La costruzione di un farmaco capace di interagire correttamente con la causadi una malattia e un compito molto complesso. Durante le fasi iniziali, dopoaver individuato il processo biologico su cui intervenire per modificare il de-corso della malattia, ricercatori e medici selezionano alcuni composti guida,ossia potenziali candidati per il principio attivo del farmaco.

Un composto guida, precursore di un futuro farmaco, viene studiato etestato meticolosamente. Una prima fase prevede che venga sperimentato sucolture di cellule, in modo da attestarne l’efficacia e il grado di sicurezzaper l’organismo. Dopo aver superato i test della fase precedente, ha inizio lostudio su animali e uomo.

La ricerca del composto guida non e semplice; si possono seguire moltestrade. Per esempio, si possono analizzare le attivita farmacologiche possedu-te da estratti di piante, si puo partire dalla valutazione degli effetti collateralidi altri farmaci gia in commercio oppure si puo seguire un metodo sperimen-tale, usato in chimica farmaceutica, consistente nello studio della relazionestruttura-attivita della molecola in esame.

Un approccio utilizzato negli ultimi vent’anni consiste nella progettazio-ne di nuovi farmaci a partire dalle simulazioni delle interazioni tra farmacoe recettore. Cio e possibile grazie alle nuove tecnologie informatiche e alprogredire delle conoscenze di farmacologia molecolare.

Attualmente esistono numerosi pacchetti software in grado di progettaremolecole in modo efficiente. Tuttavia, esiste una lacuna insita nel loro algorit-mo di costruzione dovuta non ad un errore di programmazione o all’utilizzodi modelli scorretti (cioe le soluzioni ottenute sono, pur sempre, attendibi-li), bensı all’elevato numero di falsi positivi generati dai suddetti algoritmi,che potrebbero, quindi, causare la perdita di potenziali composti guida. Inbreve, un falso positivo e una molecola con struttura tridimensionale geome-tricamente corretta accettata da un pacchetto di docking e rifiutata da unaltro. Cio e causato dai diversi approcci su cui si basano i programmi, inparticolare dal differente algoritmo utilizzato per il calcolo dell’energia dellemolecole che verifichera se il composto guida, all’interno del target proteico,

v

vi INDICE

e energeticamente stabile.Nel presente lavoro di tesi, in collaborazione con il Dipartimento di Chimi-

ca e il Dipartimento di Matematica dell’Universita di Parma, verra presentatauna possibile soluzione al problema evidenziato esibendo un programma esen-te dal difetto che andra a costruire geometricamente dei potenziali principiattivi.

Nel capitolo 1 verranno introdotte e descritte le conoscenze di base, uti-li per capire al meglio le argomentazioni trattate in questa tesi, e i metodimatematico-informatici utilizzati nella formalizzazione e nella implementa-zione del pacchetto software. Nel capitolo 2 si definiranno i requisiti funzionalidel nuovo programma. Il capitoli 3 e 4 introducono formalmente i principie le tecniche dettagliate utilizzate affinche si generino delle molecole atten-dibili. Nel capitolo 5 si elencheranno tutti i test effettuati che andranno amisurare prestazioni del programma e qualita dei candidati generati. Infine,nel capitolo 6 verranno presentati dei possibili ampliamenti e suggerimentida adottare nelle versioni successive del software.

Capitolo 1

Background

Questo capitolo si propone l’obiettivo di illustrare le conoscenze di basetrattate in questo lavoro di tesi.

Nella sezione 1.1 verranno introdotte le proteine, sostanze di fondamentaleimportanza per gli esseri viventi; il tema principale e rivolto alla loro strutturae funzionalita.

Nella sezione 1.2 verra data la nozione di legame chimico. Verra, inoltre,presentata la struttura dati grafo e alcune sue proprieta nella sezione 1.3.1;grazie ad esso e possibile dare una rappresentazione astratta al concetto dimolecola.

Nella sezione 1.4 si descrivera la tecnica di programmazione a vincoli,dandone una definizione formale; parte di queste tecniche sono state utilizzatein questa tesi con lo scopo di poter esprimere vincoli da applicare agli atomidelle molecole che andremo a generare.

La parte finale di questo capitolo (sezione 1.5) affrontera le nuove tecnichee gli approcci utilizzati per la sintesi di nuovi farmaci, il Drug Design e il Doc-king. Faremo, in aggiunta, una carrellata dei principali programmi di Dockingin commercio, marcando le loro caratteristiche e i loro limiti funzionali.

1.1 Proteine

Le proteine sono polimeri naturali composti da unita di amminoacido legatefra loro da legami peptidici. Sono sostanze di primaria importanza per lastruttura, il funzionamento e la riproduzione della materia vivente.

Gli amminoacidi si concatenano nei peptidi e nelle proteine mediante laformazione di legami peptidici fra il gruppo carbossilico di un amminoacido eil gruppo amminico in α di un altro amminoacido, rispettivamente, il carbonio

1

2 CAPITOLO 1. BACKGROUND

Figura 1.1: Legame peptidico

legato ai due ossigeni in ‘Amino acid (1)’ e l’azoto legato a due idrogeni in‘Amino acid (2)’, figura 1.1.

Tranne che nella glicina, dove R = H, il carbonio in α e un centro ste-reogeno1. La figura 1.2 riporta i venti α-amminoacidi comunemente reperibilinelle proteine. Tutti hanno un nome comune, inoltre per ognuno c’e un codicedi tre lettere ed un codice di una lettera usata nella scrittura di formule dipeptidi o di proteine.

1.1.1 Caratteristiche strutturali

Si possono individuare quattro livelli di struttura di una proteina:

• Struttura primaria

• Struttura secondaria

• Struttura terziaria

• Struttura quaternaria

1Un atomo di carbonio legato a quattro gruppi diversi si definisce asimmetrico o chiralee costituisce un centro stereogeno

1.1. PROTEINE 3

Figura 1.2: Amminoacidi


Struttura primaria

La struttura primaria e formata dalla sequenza specifica degli amminoacidi.Lo scheletro peptidico di tale struttura e il risultato della regolare succes-sione di tre atomi −N − C − C− appartenenti, rispettivamente, al gruppoamminico, all’atomo di carbonio α e al gruppo carbossilico di ogni residuo.

I livelli superiori della struttura di una proteina, dati dalle modalita diripiegamento locale e dal ripiegamento dell’intera molecola, conferiscono allaproteina la forma finale biologicamente attiva; tuttavia questi ripiegamen-ti derivano dalla struttura primaria. Le differenti proprieta associate a unaprecisa sequenza amminoacidica determinano il modo in cui la proteina puoruotare o ripiegarsi assumendo una specifica e stabile struttura tridimensio-nale che la distingue da tutte le proteine. I legami coinvolti nella strutturaprimaria sono covalenti mentre nei livelli successivi sono presenti legami aidrogeno, piu deboli dei primi.

Struttura secondaria

Come descritta in [3], la struttura secondaria di una proteina e data dallemodalita di ripiegamento della catena polipeptidica. Esistono principalmentetre2 tipi diversi di struttura secondaria: α-elica, β-foglietti e loop.

L’α-elica e una spirale destrorsa nella quale i gruppi R si proiettano al-l’esterno dello scheletro peptidico perpendicolarmente all’asse dell’elica. Lastruttura a elica di un polipeptide e stabilizzata dalla formazione di legamia idrogeno tra gli atomi di idrogeno del gruppo amminico di un residuo am-minoacidico e gli atomi di idrogeno del gruppo carbonile di un altro residuo(figura 1.3). Quando questo modello di legami a idrogeno si ripete regolar-mente lungo un segmento di catena polipeptidica, viene a stabilizzarsi lastruttura ripiegata ad α-elica.

La struttura a foglietto β a pieghe si forma quando due o piu catenepolipeptidiche sono quasi completamente distese e giacciono l’una accantoall’altra. Il foglietto e stabilizzato da legami a idrogeno che si formano tra igruppi amminici di una catena e i gruppi carbonili dell’altra. Questa strut-tura puo essere formata da catene polipeptidiche diverse oppure da differentiregioni della stessa catena polipeptidica la quale si ripiega su se stessa.

Struttura terziaria

Affinche la molecola assuma la caratteristica struttura compatta e necessarioche la catena polipeptidica cambi direzione in corrispondenza di particolari

2Esistono altri tipi di strutture, ma comunque meno frequenti.

1.1. PROTEINE 5

Figura 1.3: Dettaglio di una α-elica.

punti, ripiegandosi in varie direzioni. La struttura terziaria rappresenta ladisposizione, l’organizzazione nello spazio che una proteina assume in dipen-denza della sua specifica struttura primaria. La stabilizzazione della strut-tura, secondo [3], e data dalle interazioni chimiche tra i gruppi R (le catenelaterali dei residui amminoacidici).

• Tra specifici residui di cisteina possono formarsi legami covalenti di-solfuro che contribuiscono a mantenere il corretto ripiegamento di unacatena polipeptidica.

• Le catene laterali idrofobiche possono aggregarsi all’interno della mo-lecola proteica venendo escluse dal contatto con le molecole di acqua econtribuendo al processo di ripiegamento.

• Le forze di van der Waals possono stabilizzare le strette interazioni traresidui idrofobici.

• Legami ionici tra catene laterali con carica positiva e negativa situateall’interno di una proteina, lontano dal contatto con le molecole diacqua, possono formare ponti salini.


Struttura quaternaria

Molte proteine nella loro forma funzionalmente attiva sono formate da dueo piu catene polipeptidiche, dette subunita, ognuna delle quali e ripiegatain modo da assumere la propria peculiare struttura terziaria. La strutturaquaternaria e il risultato del modo in cui le subunita proteiche si associanoe interagiscono nell’intera proteina. Come descritto in [3], l’emoglobina e unchiaro esempio di proteina con struttura quaternaria (figura 1.4). Interazioniidrofobiche, forze di van der Waals, legami a idrogeno ionici stabilizzano l’as-sociazione delle quattro catene polipeptidiche che costituiscono la molecoladell’emoglobina.

Figura 1.4: Esempio di struttura quaternaria: l’emoglobina. E possibileidentificare quattro sub-unita, due in rosso e due in blu.

La specifica forma delle proteine permette loro di legare non covalente-mente altre molecole e cio, a sua volta, e seguito da altri eventi biologicamenteimportanti come ad esempio:

• una sostanza puo penetrare all’interno di una cellula legandosi a unaproteina trasportatrice

1.2. LEGAME CHIMICO E FORZE INTERMOLECOLARI 7

• una reazione chimica puo essere accelerata quando una proteina enzi-matica lega uno dei reagenti

• segnali chimici come gli ormoni possono legarsi a proteine presenti sullasuperficie esterna di una cellula

L’acquisizione di informazioni riguardanti la struttura tridimensionale eil processo con cui una proteina si ripiega (protein folding) e di notevole im-portanza. Grazie a cio e possibile fare delle inferenze funzionali sulla proteinaa partire dal dogma fondamentale della biologia:

Struttura ⇐⇒ Funzione

nel senso che ad ogni diversa organizzazione strutturale della proteina eassociata una specifica funzione biologica.

La conoscenza dell’esatta forma di una molecola proteica e di cio che visi puo legare e fondamentale non solo per comprendere la biologia di base,ma anche in altri campi come la medicina. Un esempio fu la determinazionedella struttura tridimensionale di una proteina essenziale per la replicazionedel virus HIV che permise la progettazione di specifiche molecole capaci dilegarsi a questa e di bloccarne l’azione ([3]).

1.1.2 Sito attivo

Il sito attivo e una porzione di enzima3 implicata nella formazione di legamicon substrati4, che daranno luogo ad una reazione chimica. Il sito attivo diun enzima si trova solitamente in una tasca di quest’ultima ed e rivestitoda residui amminoacidici che partecipano al riconoscimento del substrato; eanche il punto su cui agiscono gli inibitori enzimatici. I substrati si legano alsito attivo per mezzo di legami chimici (si faccia riferimento alla sezione 1.2)come quello idrogeno, covalente o tramite interazioni idrofobiche, andando aformare un complesso.

1.2 Legame chimico e forze intermolecolari

Come definito in [4], il legame chimico e una connessione tra atomi. Si formatra due atomi se la risultante disposizione dei nuclei e degli elettroni possiedeenergia minore di quella totale corrispondente ai due atomi separati. Se la

3L’enzima e una proteina che catalizza reazioni chimiche4Molecole su cui agisce un enzima


minore energia si puo conseguire trasferendo completamente uno o piu elet-troni da un atomo all’altro, si formano ioni e il composto sara tenuto insiemedall’attrazione elettrostatica tra tali ioni, che viene definita legame ionico.Nel caso in cui l’energia minore si possa conseguire condividendo elettroni,gli atomi si congiungeranno tramite un legame covalente e si formeranno mo-lecole distinte. Un terzo tipo di legame e il legame metallico, che vede ungrande numero di cationi5 vincolati da una mare di elettroni.

Un legame tra due atomi puo essere di tre tipi: semplice, doppio o triplo.I legami che andremo a considerare nelle rotazioni saranno tutti legami cova-lenti del primo tipo, ovvero semplici. I restanti due, a causa della loro naturachimica, non possono ruotare su loro stessi senza rincorrere allo sconvenientedi spezzare il legame.

1.2.1 Forze intermolecolari

Le forze intermolecolari sono interazioni deboli di natura elettrostatica tramolecole neutre e ioni. Le energie coinvolte in questi tipi di interazione sonodi gran lunga inferiori rispetto al legame chimico intramolecolare. A diffe-renza dei legami intratomici, le forze intermolecolari tengono unite due o piumolecole in modo non covalente. In riferimento al lavoro svolto in questa tesi,queste forze rappresentano l’interazione che piu ci interessa. Grazie ad esse,un ligando ha la possibilita di entrare a contatto con il sito attivo di unaproteina, con la conseguente formazione di un complesso.

1.3 Strutture dati fondamentali

1.3.1 Grafo

Si definisce grafo una coppia ordinata G = (V,E) tale che:

• V e l’insieme dei nodi

• E e l’insieme (a, b) : a ∈ V ∧ b ∈ V , chiamato insieme degli archi.

Si possono distinguere due tipi di grafo: grafo orientato o diretto e grafonon orientato. Nel primo caso, ogni arco specifica la direzione della connes-sione, ovvero si puo distinguere il nodo di partenza e il nodo di arrivo; nelsecondo, non c’e distinzione tra i due nodi componenti l’arco.

Definizione 1.1. L’insieme degli archiEO di un grafo orientatoO = (VO, EO)e l’insieme delle coppie ordinate (a, b) con a, b ∈ VO

5Un catione e un atomo carico positivamente.

1.4. PROGRAMMAZIONE A VINCOLI 9

Definizione 1.2. L’insieme degli archi EN di un grafo non orientato N =(VN , EN) e l’insieme delle coppie non ordinate a, b con a, b ∈ VN

Si definisce grafo semplice un grafo che non contiene archi orientati.

Definizione 1.3. Si definisce percorso un insieme di vertici v0, v1, . . . , vne una sequenza di archi (vo, v1), (v1, v2), . . . , (vn−1, vn) che li collegano; v0e vn rappresentano gli estremi del cammino.

Un percorso che abbia gli archi distinti, viene definito cammino; se vo = vnil cammino si chiama circuito o ciclo.

Definizione 1.4. Sia G = (V,E) un generico grafo e siano u, v ∈ V duegenerici vertici di G. Se esiste un cammino con estremi u e v allora i duevertici sono connessi. Inoltre, la relazione di connessione e di equivalenza.

A partire dalla relazione di equivalenza precedente, si possono definire kclassi di equivalenza chiamati sottografi e definiti come

Gi = (Vi, Ei) per i = 0, . . . , k

dove Vi ⊆ V e Ei ⊆ E. Piu semplicemente, un generico Gi e un sottografomassimale che contiene tutti gli elementi connessi tra loro. L’insieme di questisottografi prende il nome di componenti connesse di G e la sua cardinalita siindica con γ(G); ne segue che se γ(G) = 1, allora il grafo e connesso.

1.4 Programmazione a vincoli

La programmazione a vincoli e un paradigma di programmazione dove lerelazioni tra variabili possono essere dichiarate in forma di vincoli. I vinco-li differiscono dalle primitive normalmente definite dagli altri linguaggi diprogrammazione per il fatto che non specificano azioni singole da eseguirepasso-passo, ma, piuttosto, si limitano a specificare le proprieta di cui deveessere dotata la soluzione da trovare.

La nozione centrale di questa tecnica e, appunto, il vincolo (per esem-pio X > 5, X + Y < 20); esso, definito su una sequenza di variabili, esemplicemente una relazione sul loro dominio. Piu formalmente,

Definizione 1.5. [Dominio] Sia y1, . . . , yk una sequenza di variabili. Diremoche D1 × · · · × Dk e il dominio loro associato se l’insieme di tutti i valoripossibili assunti da ogni yi e proprio l’insieme Di, per i = 1, . . . , k. In formule:

yi ∈ Di, ∀i = 1, . . . , k


Definizione 1.6. [Vincolo] Consideriamo una sequenza finita di variabiliY := y1, . . . , yk, con k > 0, e D1, . . . , Dk il dominio associato loro. Un vincoloC su Y e un sottoinsieme di D1 × · · · ×Dk.

A sua volta, si definisce Problema di soddisfacimento dei vincoli o Con-straint Satisfaction Problem, da cui l’acronimo (CSP), un insieme di vincoliapplicati a un insieme di variabili.

Definizione 1.7. Per Problema di soddisfacimento dei vincoli si intendeuna sequenza finita di variabili X := x1, . . . , xk, con dominio rispettiva-mente D1, . . . , Dk, insieme ad un insieme finito C di vincoli, ognuno su unsottoinsieme di X.

Denotiamo un CSP con la dicitura 〈C;DE〉, con:

DE := xi : xi ∈ Di ∀i = 1, . . . , n

dove, C e l’insieme dei vincoli definiti sul CSP e DE l’insieme dei valoridel dominio.

Una volta formulato il problema P con l’introduzione delle variabili (edel loro rispettivo dominio) e i vincoli definiti su di esse si procede con larisoluzione di P attraverso un Solver, composto principalmente da due par-ti: propagazione dei vincoli e ricerca delle soluzioni. La risoluzione di P cipermette di stabilire se:

• il problema e consistente, cioe se ammette soluzione;

• la soluzione o le soluzioni;

Definiamo formalmente il concetto di soluzione di un CSP.

Definizione 1.8. Sia 〈C;DE〉 un CSP con DE := x1 ∈ D1, . . . , xn ∈ Dn.Diciamo che (d1, . . . , dn) ∈ D1 × · · · × Dn soddisfa un vincolo C ∈ C sullevariabili xi1 , . . . , xim se (di1 , . . . , dim) ∈ C. Inoltre (d1, . . . , dn) ∈ D1×· · ·×Dn

e una soluzione per P se soddisfa ogni vincolo C ∈ C.

Esempio 1.1. Sia 〈x < y;x ∈ [0, 10] , y ∈ [5, 10]〉; allora tutte le soluzione delproblema appena esposto sono l’insieme delle coppie (a, b) con a ∈ [0, 10] e b ∈[5, 10] tali che a < b.

1.4. PROGRAMMAZIONE A VINCOLI 11

Ricerca delle soluzioni

La ricerca delle soluzioni di un CSP e effettuata attraverso l’esplorazionedi un albero, detto albero di ricerca, dove ad ogni livello compaiono tutte lepossibili scelte di una precisa variabile. Un percorso radice-foglia rappresentauna possibile soluzione e, se questa e consistente con l’insieme dei vincoliintrodotti nel problema, allora e effettivamente una soluzione.

Propagazione dei vincoli

L’esplorazione di tutto l’albero (il quale cresce esponenzialmente rispetto alnumero delle variabili) potenzialmente, potrebbe diventare un processo moltolento. Di conseguenza, si introducono le tecniche della propagazione, il cuiscopo e di riscrivere un vincolo C in uno equivalente, con l’applicazione dideterminate regole atte a soddisfare alcune proprieta di consistenza locale.L’obiettivo e quello di ridurre l’albero di ricerca con l’eliminazione dei suoisottoalberi che sicuramente non portano a soluzioni consistenti.

Nello specifico, si hanno tre tipi di regole:

• regole di Riduzione del Dominio, come conseguenza dell’applicazionedei vincoli alle variabili

• regole di Trasformazione, cioe semplificazione dei vincoli

• regole di Introduzione, grazie alle quali si aggiungono nuovi vincoliimpliciti

1.4.1 Alberi di ricerca

Affinche si trovino tutte le soluzioni consistenti del problema, si ricorre aglialberi di ricerca e agli algoritmi per esplorarlo, algoritmi di ricerca. E doverosoprecisare che l’algoritmo di ricerca non costruisce un albero di ricerca per poiesplorarlo successivamente; in realta, l’albero viene costruito ‘al volo’ durantel’algoritmo di esplorazione.

Definizione 1.9. Dati due CSP P1 e P2 e una sequenzaX di variabili comuniai due problemi, diremo che P1 e P2 sono equivalenti se:

• per ogni soluzione d per P1, esiste una soluzione per P2 e, quest’ultima,coincide proprio con d sulla sequenza di variabili X.

• per ogni soluzione e per P2, esiste una soluzione per P1 e, quest’ultima,coincide proprio con e sulla sequenza di variabili X.


Definizione 1.10. Siano P un CSP e x1, . . . , xn una sequenza di variabili.A e un albero di ricerca finito per P se rispetta le seguenti clausole:

• tutti i suoi nodi sono degli altri CSP

• la radice e proprio P

• i nodi di un livello pari hanno esattamente un discendente diretto

• se P1, . . . ,Pm con m ≥ 1, sono discendenti diretti di P0, allora l’unionedi P1, . . . ,Pm e equivalente a P0 (in relazione a X).

Figura 1.5: Esempio di un Albero di Ricerca.

Operativamente, l’idea alla base di una visita di un albero di ricerca ela seguente: si parte dalla radice e si prosegue visitando un figlio per tuttele sue scelte possibili. La visita e ricorsiva e nel momento in cui ci si trovaad avere un vincolo inconsistente si risale verso l’alto (backtracking) in cercadi un nodo con delle scelte ancora possibili; la ricerca di tutte le decisionitermina nel nodo radice.

1.5 Drug Design

Dal momento in cui strutture tridimensionali di proteine derivate da cri-stallografia a raggi X o spettroscopia NMR, divennero disponibili, si vide la

1.5. DRUG DESIGN 13

nascita di pacchetti software in grado di utilizzare queste informazioni perdar vita al progetto chiamato Drug Design. Il problema che questi softwaredevono affrontare e quello del Docking, ovvero di predire un complesso ener-geticamente favorevole composto da una proteina e una molecola (ipoteticofarmaco), chiamato in questo contesto, ligando.

Il Drug Design ha l’obiettivo di trovare una struttura guida, una piccolamolecola che si lega ad una precisa proteina target e che puo essere studiataa fondo per diventare un farmaco.

Dal punto di vista biologico, legare una piccola molecola ad una specificaproteina significa poter inibire la sua funzione, in modo tale da poter rende-re la sua azione verso l’organismo inoffensiva o, accelerare la sua funzione,simulando l’intervento del naturale ligando.

Il metodo che sta alla base del Drug Design e il concetto di chiave-toppa(ligando-proteina): quel che si vuole ottenere e una serie di chiavi valide perquella specifica toppa.

Figura 1.6: Figura rappresentante il docking tra una piccola molecola (inmarrone) e una proteina target per dar vita ad un complesso

Molti aspetti fanno del Docking un problema difficile da risolvere. Primotra tutti quello che riguarda il problema dello Scoring, ovvero calcolare eassegnare un punteggio al modo in cui ligando e proteina si legano tra loroper andare a formare il complesso. Al giorno d’oggi, non esiste una funzionedi scoring “d’uso generale” che permetta di predire in modo accurato cio. Insecondo luogo, e necessario considerare l’elevato numero di gradi di liberta:il piu importante e quello relativo al possibile orientamento spaziale dellaproteina e del ligando e la conformazione di quest’ultimo.Inoltre, puo variare la conformazione della proteina, si possono aggiungeremolecole d’acqua tra le molecole e puo cambiare lo stato di protonazione6.Cio implica che le funzioni di scoring contengano tipicamente molti minimilocali di energia difficili da ottimizzare.

6La protonazione e una reazione che consiste nell’addizione di un protone, cioe unidrogeno carico positivamente (H+), ad un atomo, ad una molecola o ad uno ione. Laspecie protonata subisce variazioni chimico-fisiche (idrofilia, proprieta ottiche, etc..).


Malgrado non esista una soluzione generale al problema del docking, sonostati sviluppati molti algoritmi specializzati in vari aspetti del problema eapplicati in modo soddisfacente.

1.6 Algoritmi di Docking

Il Docking e una tecnica computazionale atta a generare varie conformazionidi piccole molecole all’interno del sito attivo di una proteina. Successivamenteviene applicata una funzione di scoring ai risultati, in modo tale da poterottenere una lista dei vari ligandi in base ad un certo punteggio: maggiorpunteggio uguale maggiore stabilita.

Come descritto in [1], esistono diversi approcci per risolvere il problemadel Docking.

1.6.1 Algoritmi con corpi rigidi

Questi sono stati i primi algoritmi ad essere stati sviluppati. Sia la protei-na che il ligando vengono tenuti fissi nella loro conformazione spaziale e lacomplessita del problema si riduce alla ricerca dell’orientamento, con energiaminore, tra le due molecole. Sicuramente, il punto a favore di questa classedi algoritmi e la velocita di elaborazione: pochi gradi di liberta e vincoli assi-curano delle prestazioni notevoli. Tuttavia, maggiore e la velocita di ricerca,minore e la qualita dei ligandi trovati; infatti, a causa dei pochi movimentisimulati, non si esplorano tutte le possibili combinazioni di un ligando.

1.6.2 Algoritmi con ligando flessibile

La limitazione maggiore della precedente classe di algoritmi consiste nell’igno-rare completamente la flessibilita del ligando. Infatti, spesso piccole molecolehanno uno spazio conformazionale molto ampio con livelli energetici moltobassi, e quindi stabili.

Algoritmi a costruzione incrementale

Uno degli algoritmi usati frequentemente, facente parte di questa classe, e lacostruzione incrementale, che si basa sul concetto di frammentazione. Il ligan-do viene suddiviso in numerosi frammenti, i quali costituiscono una porzionedi ligando rigida, ovvero gli atomi che lo compongono non si muovono. Inbreve, l’algoritmo funziona nel modo seguente: viene piazzato il primo fram-mento del ligando (chiamato frammento ancora) nel sito attivo della protei-

1.6. ALGORITMI DI DOCKING 15

na; quindi si procede al piazzamento in successione dei restanti frammenti,collegandoli l’un l’altro in sequenza e validandoli con il calcolo dell’energia.

In altre parole, si possono riconoscere tre fasi:

• selezione del frammento base;

• piazzamento del frammento base;

• costruzione incrementale del ligando con i frammenti successivi.

Il grado di liberta di questo algoritmo e dato dalla rotazione dei legami singoli.

Algoritmi genetici

L’algoritmo genetico e un metodo euristico di ricerca e ottimizzazione cheimita il processo di evoluzione. L’idea generale e quella di partire da un certonumero di possibili soluzioni (individui) chiamate popolazione e, a ciascunaiterazione, operare una selezione di individui, impiegandoli per generare nuo-vi elementi della popolazione stessa, cosı da costituire una nuova popolazio-ne per l’iterazione (o generazione) seguente. Tale successione di generazionievolve verso una soluzione ottimale (locale o globale) del problema assegna-to. L’evoluzione viene ottenuta attraverso una parziale ricombinazione dellesoluzioni, e l’introduzione di mutazioni casuali nella popolazione di parten-za; sporadicamente nascono individui con caratteristiche non comprese traquelle presenti nei dati della specie originaria. Finita la fase di evoluzione lapopolazione risultante viene analizzata e vengono tenute solo le soluzioni chemeglio risolvono il problema: gli individui con le qualita piu adatte all’am-biente in cui si trovano hanno quindi maggiori possibilita di sopravvivere eriprodursi. Queste soluzioni subiranno una nuova fase di evoluzione e cosıvia. Questa classe di algoritmi viene utilizzata in vari campi; per esempio,in biologia molecolare e utilizzato per predire l’adattamento di un genomaall’ambiente, riconducendosi quindi all’evoluzione della specie introdotta daCharles Darwin.

1.6.3 Docking attraverso simulazione

Al contrario dei metodi appena menzionati, esistono algoritmi che affronta-no il problema attraverso tecniche di simulazione. Questi algoritmi partonoda una conformazione iniziale per poi passare ad altre ad energia minore,attraverso piccoli movimenti effettuati alla struttura e scartando quelle piuinstabili. Alcuni algoritmi facenti parte di questa categoria sono Simulazio-ni di Dinamica Molecolare, algoritmi di Monte-Carlo, Metodi ibridi ottenuticombinando due o piu algoritmi.


1.6.4 Panorama del mercato attuale

La scelta di un pacchetto di Docking e una operazione tutt’altro che semplice;prima di procedere all’acquisto, le case farmaceutiche eseguono particolaritest su uno spettro molto ampio di complessi. Nel 2006 e stato condottouno studio dalla GlaxosmithKline [7] volto a fare una indagine sullo statoattuale delle tecniche computazionali per il drug design ed in particolaresul docking e sulle funzioni di scoring. Da pochi anni a questa parte, sonostati pubblicati un numero sempre crescente di valutazioni su pacchetti didocking e funzioni di scoring, includendo recensioni sulle nuove tecniche,facendo confronti tra piu programmi di docking e studiando le correlazionitra i punteggi formulati dal docking e i punteggi calcolati dall’affinita delcomposto (quanto due molecole si legano bene). Questo, invece, differiscedalle solite valutazioni per due motivi: primo, si misurano le performace dimolti pacchetti di docking su numerosi target, secondo, il set dei composti perogni target e costituito da un gran numero di composti relativamente correlatitra loro, per cui le affinita sperimentali sono state misurate utilizzando unprotocollo standard, sviluppato solitamente dallo stesso gruppo di ricerca.

Sono stati presi in esame 10 differenti pacchetti; in particolare, alcuni diquesti offrono piu funzioni di scoring o algoritmi di docking per un totale di19 protocolli.Il focus della valutazione si basa su tre usi tipici di questi pacchetti:

1. predizione delle conformazioni di piccole molecole all’interno della pro-teina target;

2. virtual screening di banche dati volto all’identificazione di composti pervari target di proteine;

3. predizione dell’affinita dei composti

Risultati dello studio valutativo

Per quanto riguarda la predizione delle conformazioni di piccole molecole, siottengono dei buoni valori per tutti i tipi di target proteici; in particolare, pertutti i target (tranne uno), almeno un programma riesce a posizionare piu del40% dei ligandi entro i 2A rispetto la struttura cristallina7. Infatti, per moltitarget proteici, il 90% dei relativi ligandi potrebbe avere una orientazionecorretta; da cio si deduce che gli algoritmi di docking riescono ad esplorare lo

7Questo e un buon risultato: significa che la distanza media tra la molecola originata ela struttura cristallina e minore di 2A

1.6. ALGORITMI DI DOCKING 17

spazio conformazionale sufficientemente bene da ottenere delle buone confor-mazioni. Il problema e di natura diversa: non esiste nessun programma cheabbia valore generale, cioe che ottenga buoni risultati per un qualsiasi targetproteico.

Per quanto riguarda il Virtual Screening, questo ha successo quandosi utilizzano dei dati che simulano una tipica classe di composti farma-cueutici. Inoltre, in assenza di informazioni a priori sulla proteina target,le performance dei programmi sono in contrasto rispetto ai tipi di targetvalutati.

Capitolo 2

Scopo della Tesi

Uno degli aspetti comuni a tutti i pacchetti di Docking esistenti e che un unicoalgoritmo racchiude due aspetti completamente differenti: posizionamento delligando e calcolo dell’energia del sistema.

Una conseguenza che deriva dall’aspetto appena menzionato e che unchimico-farmaceutico, in possesso di una propria funzione energetica F , nelmomento in cui decidesse di applicarla, non avrebbe a disposizione un setdi ligandi “neutrali”: la sua funzione infatti verrebbe applicata ad un set dimolecole gia precomputate da un’altra funzione energetica F ′. In particolare,potrebbero mancare delle conformazioni che sarebbero state accettate da Fma rifiutate da F ′.

Un approccio completamente innovativo (e quindi sperimentale) potrebbeessere quello di separare in due parti un algoritmo di docking. L’idea e quelladi distinguere in modo netto la generazione di tutte le possibili conformazio-ni di un ligando all’interno del sito attivo della proteina target, ignorandocompletamente l’energia del sistema, dalla valutazione energetica.

2.1 Ligand-rotation

A fronte di tutto quello che si e detto nel capitolo precedente, quel che vo-gliamo e un pacchetto di Docking in grado di poter restituire delle confor-mazioni valide, cioe che rispettino la geometria della chimica, senza dovercipreoccupare dello stato energetico del complesso.

Lo scopo primo di questa tesi e dunque quello di studiare, analizzare,modellare e implementare un programma de novo di Docking, ligand-rotation,scritto in C++, in grado di poter generare delle conformazioni valide dalpunto di vista geometrico.

19

20 CAPITOLO 2. SCOPO DELLA TESI

Le tecniche alla base del funzionamento di ligand-rotation sono propriedella programmazione a vincoli. Nel capitolo 3 si andra a definire formalmen-te il problema CSP associato e quindi, l’insieme delle variabili e dei vincoli;si mostrera anche in che modo ligand-rotation assicura la consistenza deivincoli rispetto ai movimenti del ligando. Inoltre, l’esplorazione dello spaziodi ricerca viene eseguito efficacemente grazie ad algoritmi implementati ad-hoc per ligand-rotation, come la visita dell’albero di ricerca e l’algoritmo dibacktracking. Abbiamo ricorso ad una progettazione e realizzazione a manodel motore di ricerca, in quanto, gli oggetti del nostro algoritmo (i punti nelletre dimensioni) non sono trattati nativamente dai risolutori esistenti. Osser-viamo che, seppure il tutto si basi su tecniche di programmazione a vincoli,ligand-rotation non implementa la propagazione dei vincoli. La causa di cioe da imputarsi al tempo: avendone a disposizione poco, si e optato per unaimplementazione senza propagazione. Per ulteriori dettagli e approfondimen-ti si veda il paragrafo 6.2.2 nella quale si descrive una futura espansione diligand-rotation.

Il funzionamento si basa su un modello che prevede la proteina rigida e illigando flessibile, come la maggior parte dei pacchetti di Docking presenti sulmercato. Una motivazione valida di quanto appena detto e che la complessitadell’algoritmo cresce esponenzialmente all’aumentare del numero dei gradidi liberta (sezione 4.5): trattare anche i possibili movimenti della proteina,renderebbe ligand-rotation piu costoso.

L’algoritmo centrale di tutto ligand-rotation prevede di partire da unpotenziale candidato ligando e da una proteina con cui dovra interagire; siandra, quindi, a simulare i possibili movimenti che il ligando puo compiereall’interno del sito attivo della proteina. L’output sara composto da una seriedi molecole (in formato mol2); esse potranno essere studiate dal punto divista energetico, cioe valutate da una opportuna funzione di Scoring, la qualeandra a selezionare i ligandi che, inseriti nel contesto proteico, risulterannoenergeticamente favorevoli.

Evidenziamo il fatto che ligand-rotation e molto flessibile: i suoi parame-tri sono numerosi, molti dei quali definibili dall’utente. Attraverso un file diconfigurazione, l’utente e in grado di modificare il comportamento di ligand-rotation in base alle sue esigenze e al tipo di risultato che vuole ottenere.Ad esempio, si puo specificare il grado di qualita e precisione attraverso variparametri che andranno a modificare il campionamento delle strutture datiutilizzate o il numero dei gradi di liberta del sistema. A seconda delle scelteeffettuate, si andra ad incidere notevolmente sulla complessita computazio-nale e quindi sull’efficienza del software stesso, in termini di tempo impiegatodal pacchetto nel restituire una serie di ligandi.

Nei prossimi capitoli si valuteranno le prestazioni ottenute da ligand-

2.1. LIGAND-ROTATION 21

rotation in fase di test. A partire dall’insieme dei ligandi ottenuti, la valuta-zione vertera sui seguenti aspetti:

• complessita in tempo, in relazione al tipo di configurazione impostata

• qualita strutturale delle molecole ottenute

• generazione di un file di configurazione ottimale

Nel capitolo 5 descriveremo ampiamente i test effettuati su ligand-rotation.

2.1.1 Analisi dei requisiti

Il primo compito dell’analisi dei requisiti e identificare tutti i requisiti delsistema software. Questi vengono documentati in piu iterazioni con il com-mittente, in modo da chiarire eventuali dubbi e stendere gradualmente unadocumentazione. In questa fase si analizzano i punti di vista dell’utente fi-nale, esterno, senza occuparsi dei dettagli implementativi informatici, i qualirappresentato il punto di vista dell’utente interno. La difficolta maggiore diquesta fase e costituita dai problemi di comunicazione: lo scambio di infor-mazioni, pur su argomentazioni semplici, tra analista e cliente finale, si puorivelare complesso. Trovare un linguaggio comune tra i due e scambiare leproprie idee in modo comprensibile ad entrambi e di primaria importanza.

Ligand-rotation e un progetto software sviluppato in collaborazione conil Dipartimento di Chimica dell’Universita di Parma. L’idea e quella di rea-lizzare un pacchetto di Docking in grado di lavorare in maniera indipendentedall’energia del sistema ligando-proteina. In occasione di numerosi colloqui,abbiamo stilato i requisiti funzionali che ligand-rotation avrebbe dovuto ri-spettare. Innanzitutto, abbiamo definiti i gradi di liberta che il software deveimplementare. In primis, la molecola si muove grazie a degli step di rotazioneapplicati ai legami di tipo semplice. In secondo luogo, ligand-rotation deveessere in grado di esplorare il sito attivo della proteina.

Affinche il programma possa essere utilizzato da utenti non prettamenteinformatici, questi comunicheranno con ligand-rotation attraverso un file diconfigurazione. Si potranno specificare sia opzioni riguardati le funzionalitadi base, come input proteina, input ligando, etc.., sia parametri riguardantile prestazioni del software.

22 CAPITOLO 2. SCOPO DELLA TESI

Capitolo 3

Formalizzazione

In questo capitolo si esamineranno le tecniche e i principi utilizzati nell’im-plementazione di ligand-rotation, con lo scopo di conseguire gli obiettivi pre-fissati nei precedenti capitoli. Ligand-rotation iniziera la sua computazioneconsiderando il ligando base, ricevuto in input. La generazione delle confor-mazioni avra luogo effettuando sui legami tutte le rotazioni possibili, questeultime ottenute con una combinazione di rotazioni indipendenti individua-te da una visita opportuna della molecola. Si andra, quindi, a testare ognipossibilita di rotazione, verificando con i vincoli se queste realizzano un ligan-do valido geometricamente. Infine, sapendo generare tutte le conformazioni,queste saranno posizionate in uno spazio definito dall’utente in termini ditraslazioni e rotazioni. Questo affinche ligand-rotation possa analizzare lacavita del sito attivo. L’insieme di tutti questi gradi di liberta andranno adescrivere le diverse conformazioni ottenibili.

Innanzitutto, nella sezione 3.1 verra formalizzato il problema P , esplici-tando l’insieme delle variabili e dei vincoli rappresentanti il modello del CSPrelativo al problema P di ligand-rotation. La nozione di vincolo e intrinse-camente legata all’algoritmo che genera le conformazioni. Nella sezione 3.2verra presentato l’algoritmo di visita dello spazio di ricerca delle combina-zioni dove, per ogni atomo, si specificheranno l’insieme delle scelte possibilirappresentate dal numero di rotazioni che il legame relativo a quell’atomopuo compiere; l’idea e quella di trovare tutte le possibili conformazioni diun dato ligando. L’insieme delle sue soluzioni andra a costruire, a sua volta,l’insieme dei ligandi validi. Considerando le loro variabili, queste ultime do-vranno andare a verificare, per definizione di soluzione di CSP, ogni vincolodefinito per il problema P .

La sola generazione di queste ultime, nel sistema di riferimento del li-gando, e molto restrittivo nel senso che l’insieme delle soluzioni difficilmen-te ne conterra uno posizionato abbastanza correttamente rispetto a quello

23

24 CAPITOLO 3. FORMALIZZAZIONE

minimizzato, scaricato dalla banca dati PDB1.Una possibile soluzione presentata nella 3.5, implementata in ligand-

rotation, consiste nel posizionamento in sequenza dei ligandi ottenuti pre-cedentemente in tanti sistemi di riferimento specifici della proteina, andan-do a costruire ligandi con orientamento spaziale differente in modo tale daesplorare il sito attivo della proteina.

3.1 Modello del problema a vincoli

L’aspetto chiave di questo progetto e proprio quello di CSP, Constraint Sati-sfaction Problem. Definiamo quindi formalmente il modello utilizzato per larappresentazione del problema P riguardante ligand-rotation.

3.1.1 Variabili

Si possono individuare due tipi di variabili che compongono il problema P :

• variabili di tipo legame, ovvero le k rotazioni che caratterizzano unospecifico legame

• variabili di tipo punto, ovvero le coordinate corrispondenti agli atomidel ligando e della proteina

Le variabili di tipo legame sono associate ad ogni legame da ruotare nelligando di partenza. Si procede all’identificazione delle suddette variabili perpoi effettuare delle rotazioni successive che modificano la struttura della mole-cola. Intuitivamente, la rotazione di un legame influenza le due sottomolecoleseparate dal legame ruotato. Da cio, si deduce che non tutti i legami potran-no effettuare delle rotazioni, in quanto, se il legame in questione appartienead un anello, si va incontro alla rottura di quest’ultimo. Le variabili punto,invece, identificano tutti i possibili piazzamenti di un particolare atomo aseguito di rotazioni della sottomolecola a cui appartiene.

3.1.2 Dominio delle variabili

Il dominio di una variabile di tipo legame e un insieme di interi che identi-fica il numero di rotazioni che un legame puo effettuare. Questo valore puovariare a seconda della natura del legame chimico o delle impostazioni det-tate dall’utente. In generale, si possono riassumere i seguenti casi: sia l ungenerico legame,

1Il Protein Data Bank, o PDB, e un archivio per dati di struttura in 3D di proteine eacidi nucleici.

3.2. VISITA DELLO SPAZIO DI RICERCA 25

• se l compare all’interno di un ciclo aromatico, collega un atomo diidrogeno o e stato bloccato dall’utente, si definisce un’unica rotazione,ovvero l’identita

• in tutti gli altri casi, l effettua un numero di rotazioni pari a k, definitodall’utente; formalmente, l’insieme R degli angoli di rotazione e:

R =

2π

ki : i ∈ 0, . . . , k − 1

3.1.3 Vincoli implementati

Per meglio comprendere la natura chimica del vincolo che si vuole implemen-tare, consideriamo il legame piu debole esistente, legame di Van der Waalse immaginiamo una sfera di raggio pari alla lunghezza del legame con centrosull’atomo che stiamo esaminando. Affinche due atomi non entrino in colli-sione e necessario che ci sia una certa distanza tra i due, ovvero che le sfereidentificate dai due atomi non si intersechino mai o a meno di una certa sogliao tolleranza.

A partire da questa considerazione, nasce il vincolo di non sovrapposizio-ne, il quale verifichera la consistenza di un movimento del ligando. Siano N ,NP rispettivamente gli insiemi degli atomi del ligando e della proteina. Siaa ∈ N e sia ra il suo raggio atomico, allora ∀n ∈ (N r a) ∪NP

‖a− n‖ ≥ ra + rnj∀j = 1, . . . , |(N r a) ∪NP |

dove rnje il raggio atomico di nj.

L’esecuzione del test di sovrapposizione richiede un tempo quadratico sulnumero dei nodi. Osserviamo che, per motivi di efficienza, il controllo nonvertera su tutti gli atomi appartenenti all’insieme (Nra)∪NP ; nella sezione4.4 descriveremo l’implementazione lineare utilizzata da ligand-rotation perla verifica dell’unico vincolo esistente e si descriveranno le strutture datiutilizzate per poter ridurre l’insieme degli atomi da verificare.

3.2 Visita dello spazio di ricerca

Il ligando si presta bene ad essere rappresentato da una struttura dati ditipo grafo non orientato. Sia L = (N,A) il grafo in questione (illustratoin figura 3.4) dove N rappresenta l’insieme degli atomi e A l’insieme deilegami chimici esistenti tra atomi; e sia P = (NP , AP ) il grafo associato allaproteina target. Le informazioni relative ai nodi e agli archi possono essere


ricavate facilmente dai file di estensione mol2, formato standard utilizzatoper la rappresentazione di molecole (si faccia riferimento alla sezione 4.1.2per ulteriori approfondimenti).

Grazie ad una visita sul grafo L, e possibile costruire un albero Z (figura3.5) definito spanning tree (albero di copertura) per L, dove l’insieme deisuoi nodi coincide con N mentre l’insieme degli archi R dell’albero e unsottoinsieme di A. L’albero Z ha lo scopo di dettare un ordine di visita degliatomi della molecola al generatore di conformazioni. Cio da la possibilita dieffettuare il controllo di consistenza subito dopo la visita di un nodo, dato chei successivi non andranno a influire quello appena visitato (cio sta a significareche esso sicuramente non subira ulteriori rotazioni). Questa tecnica si traducequindi, in un risparmio in termini temporali: una visita topologica permette,infatti, di minimizzare il numero di volte che un atomo viene ruotato. Senzal’ausilio dell’albero, o di una opportuna visita della molecola, e visitando inodi in ordine casuale, gli atomi si stabilizzano piu tardi, dopo aver fattomolte piu rotazioni del necessario: cio a causa della dipendenza esistente tragli atomi. La verifica di consistenza si sarebbe dovuta applicare dopo averpiazzato tutti i nodi con un conseguente aumento dei tempi di esecuzionee rendendo ligand-rotation un programma inefficiente; inoltre, avrebbe resovano l’utilizzo di tecniche di programmazione a vincoli.

Grazie all’ordine stabilito dalla visita, e possibile costruire una lista D diatomi, la quale verra iterata dal generatore durante la costruzione delle mo-lecole. Ogni elemento della lista avra associato il suo sottoalbero, dipendentedallo spanning-tree effettuato su L. Da una rotazione dell’elemento Dk neconseguira una di tutto il suo sottoalbero.

3.2.1 Albero di ricerca

Per meglio comprendere in che modo vengono generate tutte le possibiliconformazioni, introduciamo il concetto di albero di ricerca. Daremo unadefinizione costruttiva a livelli dell’albero di ricerca T di ligand-rotation. Alivello 0 compare l’elemento D0. Al livello 1 compare D1 in tutte le possibiliscelte dettate dalla sua variabile di tipo legame, cioe in tutte le possibilirotazioni che, il legame caratterizzato dal nodo x ∈ Z corrispondente a D1 eil padre y ∈ Z di x, puo effettuare.

Il numero dei nodi che compaiono a livello k e uguale alle possibili sceltedella variabile legame dell’elemento Dk moltiplicata al numero di nodi pre-senti a livello k − 1 di T . Ne segue che ogni nodo del livello k − 1 ha unnumero di figli pari al numero di possibili scelte di Dk.

3.3. GENERAZIONE DELLE CONFORMAZIONI 27

Avremo, dunque, una corrispondenza tra gli alberi Z e T . Un nodo n ∈ Tcorrisponde ad una scelta di rotazione dell’arco di a ∈ Z, o legame chimico,che collega il nodo di Z corrispondente a n e suo padre.

Nelle figure 3.1 e 3.2 sono raffigurati due esempi di albero Z e T .

Figura 3.1: Albero Z

Conseguentemente, un cammino radice-foglia dell’albero T corrispondead una possibile soluzione e gli atomi compaiono nello stesso ordine dellalista D.

3.3 Generazione delle conformazioni

L’idea alla base della generazione delle conformazioni e quella di visitare l’al-bero di ricerca T andando ad esplorare tutti i suoi cammini, ovvero esaminaretutte le possibili scelte delle variabili di tipo legame di tutti gli elementi diD.

Il caso base e dato dall’esplorazione del nodo D0; sappiamo che, i no-di successivi non andranno a modificare la sua posizione in quanto radicedell’albero Z. Verifichiamo dunque la consistenza del vincolo su di esso: senon avviene violazione, si prosegue con il nodo successivo D1, altrimenti lacomputazione termina.

Il passo induttivo e caratterizzato dall’esplorazione del k-esimo elementodella lista, Dk, il quale si trova al livello k di un generico cammino dell’albe-ro T . Dalle considerazioni fatte precedentemente, sappiamo che gli elementidell’insieme Ik = D1, . . . , Dk, ruotati con un opportuno angolo delle ri-spettive variabili di tipo legame, sono stati validati e quindi, non subirannomodifiche; ne segue che Dk, essendo il figlio di un elemento di Ik gia piaz-zato, non subira ulteriori violazioni. Verifichiamo, dunque, la consistenza del


Figura 3.2: Albero T


Figura 3.3: Struttura molecolare del tamoxifene renderizzata con MolegroMolecular Viewer.


Figura 3.4: Grafo rappresentante la molecola in figura 3.3. I numeri al-l’interno dei nodi rappresentano l’id degli atomi mentre quelli negli archirappresentano l’id dei legami chimici. Grafo costruito con Graphviz.


2 8

2 6

1

2 7

1

2 5

5

2 4

5

2 3

5

2 2

5

2 1

1

2 9

2 0

1

2

1

1

1

3

5

1

4

1

1 3

5

5

5

7

5

1 4

1

1 9

6

1

8

1

1 2

9

1

1 0

1

1 1

1

1

1 5

1

1 6

1

1 7

1

1 8

1

1

Figura 3.5: Albero ottenuto da una visita sul grafo in figura 3.4. Sui nodisono visibili gli id degli atomi mentre agli archi sono associati il numero dellerotazioni che ogni legame puo compiere. Gli archi tratteggiati rappresentanoquelli effettivamente presenti sul grafo, ma eliminati per evitare cicli. Alberocostruito con Graphviz.


vincolo su Dk. Proseguendo con l’elemento Dk+1, del livello k + 1 di T , siaprono un numero di scelte pari al numero di rotazioni assegnate alla suavariabile legame, la quale corrisponde all’arco dell’albero Z che collega ilvertice x ∈ Z (identificato da Dk+1) e suo padre y ∈ Z. In questo puntodell’operazione, non e necessario conoscere la scelta esplorata sul nodo Dk,ovvero non e importante sapere quale particolare cammino stiamo seguendo;esploriamo, dunque, la prima scelta con conseguente verifica del vincolo.

Se ha esito positivo, si piazza l’atomo Dk+1, aggiornando l’insieme Ik, chediventa Ik+1 = Ik∪Dk+1 e proseguendo con Dk+2. In caso contrario, graziealla procedura di backtracking, si seleziona una scelta non ancora esplorataper Dk+1.

Giunti ad una foglia di T , ovvero all’elemento Dn, la quale si rivela con-sistente rispetto al vincolo implementato, il nuovo insieme In = In−1 ∪Dncontiene tutti gli atomi consistenti, posizionati rispetto alla particolare va-riabile di legame a loro corrispondente. Le loro coordinate attuali rappresen-tano una soluzione del problema P . A questo punto, si genera un file mol2,rappresentante una conformazione valida della molecola.

Si prosegue ricercando le altre soluzioni, attivando il backtracking cherisale lungo la lista e quindi lungo l’albero di ricerca T , per analizzare lescelte non ancora intraprese.

La procedura globale termina quando il backtracking arriva alla radice.

3.4 Rotazione di un legame

Le rotazioni dei legami rappresentano uno dei gradi di liberta di ligan-rotation. Attraverso di esse e possibile muovere strutturalmente il ligandoe, con l’intervento della consistenza dei vincoli, si andra a valutare geometri-camente se il movimento e accettabile o meno.

E doveroso precisare alcuni dettagli qui di seguito. Solitamente un ligandoha un numero elevato di legami; ne segue che non tutti ruoteranno. Cio edovuto al fatto che alcuni di essi vanno a formare particolari composti. Unesempio sono gli anelli aromatici, composti organici a struttura planare i cuilegami non sono singoli; cio implica che il composto e molto stabile e unarotazione dei suoi legami sarebbe molto improbabile.

Inoltre, eviteremo di far ruotare i legami nei quali in una delle due estre-mita compare un idrogeno. Si tratta di una scelta implementativa fatta so-stanzialmente per due motivi: in primis, la rotazione di un idrogeno none necessaria perche questa non apporta benefici all’energia del complessoproteina-ligando e, in secondo luogo, la rotazione di molti legami aumenta lacomplessita in tempo di ligand-rotation (si veda la sezione 4.5).

3.4. ROTAZIONE DI UN LEGAME 33

a

b d

c y

e x

l

f g

h

Figura 3.6: Esempio di albero Z. In rosso sono evidenziati i nodi che ca-ratterizzano il legame ruotante l e in blu gli elementi di B che subiranno larotazione.

Detto cio, possiamo ora dare la definizione formale di rotazione.Sia L l’insieme dei legami ruotabili, ovvero l’insieme dei legami tali che:

• non compaiano in composti aromatici;

• gli atomi caratterizzanti il legame non sono degli idrogeni

Poiche un legame identifica un particolare arco dell’albero Z, la sua ro-tazione andra a modificare tutto il sottoalbero con radice proprio il figliodell’arco. Sia l ∈ L un legame che ammette rotazione, x, y rispettivamenteil padre e il figlio del corrispondente arco di A e sia B l’insieme dei nodidel sottoalbero radicato in y (figura 3.6). L’obiettivo e quello di ruotare ogninodo di B. Ad esempio, data la rotazione del legame l di figura 3.6, l’insiemeB e composto dagli elementi f, g, h.

Siano X, Y i punti dello spazio associati ai nodi x, y. Sia inoltre SY

l’insieme dei punti dello spazio associati all’insieme dei nodi B. Sia P ∈ SY

generico e sia t la retta passante per i punti X e Y che identifica l’asse dirotazione.

Si definisce rotazione di Φ gradi su t la rotazione effettuata su tutti gliatomi appartenenti a SY .


Figura 3.7: Formula della Rotazione descritta in [8]

Operativamente, la rotazione viene effettuata nel seguente modo: sia P ∈SY come mostrato in figura 3.7 e sia s la retta passante per N ed O. Calcoloil versore n attraverso la formula:

n =N −O‖N −O‖

Applicando ora la formula 3.2 al punto P , otterremo la retta r′ e quindi ilpunto Q. Procedendo con questo algoritmo per tutti i punti di SY otterremola rotazione lungo s di angolo Φ di tutto il sottoalbero.

r′ =−−→ON +

−−→NV +

−−→V Q (3.1)

= r cos Φ + n(n · r)(1− cos Φ) + (r × n) sin Φ (3.2)

3.5 Esplorazione del sito attivo della proteina

La generazione delle combinazioni del ligando (presentata nella sezione 3.3)non e sufficiente per i nostri scopi. Essa ha luogo sul sistema di riferimento delligando, ignorando completamente la struttura della proteina. In altre parole,le semplici rotazioni lungo i legami del ligando non bastano, ma sono necessarialtri gradi di liberta che permettano a ligand-rotation di poter esplorare lacavita della proteina.

L’idea implementata in ligand-rotation consiste nella generazione di unaserie di basi ortonormali che andranno a campionare i tre gradi di liberta e,per ognuno di essi, si andra a testare diverse origini sul quale posizionare l’a-tomo D0. L’applicazione di queste basi al nostro ligando di partenza su tutte

3.5. ESPLORAZIONE DEL SITO ATTIVO DELLA PROTEINA 35

le origini definite dall’utente daranno luogo a una serie di ligandi, chiamatiligandi base, dall’orientamento spaziale molto vario. Avviando la ricerca delleconformazioni su ognuno di questi ligandi, ligand-rotation sara in grado diottenere risultati molto soddisfacenti.

Definiamo ora formalmente questi concetti.

3.5.1 Construzione delle basi ortonormali

Consideriamo l’insieme H ⊂ R3 che identifica n punti su una sfera di raggiounitario e su ognuno di essi costruiamo un sistema di riferimento. Effettuandok rotazioni di 2π/k lungo l’asse individuato dal versore uscente da ogni puntodi H si andranno a generare k sistemi di riferimento su ogni punto, con untotale di n · k.

Sia c il versore uscente da un generico punto h ∈ H e sia B = (0, 0, 1) unversore della base canonica di R3. Chiamiamo a e b rispettivamente i versoriottenuti dal prodotto vettoriale tra c e B e tra c e a, ovvero:

a = c× Bb = c× a

(3.3)

Nel caso in cui c sia proprio il versore B, i versori a e b avranno compo-nenti rispettivamente (1, 0, 0) e (0, 1, 0).

I versori a, b e c identificano il sistema di riferimento ortogonale R conorigine nel punto h. Eseguendo k step di rotazioni lungo l’asse c, ognuna di2π/k, otteniamo i k sistemi di riferimento.

Eseguendo questo procedimento per ogni punto della sfera otteniamo m =n · k sistemi di riferimento R1, . . . ,Rm.

3.5.2 Generazione ligandi base

Avendo ora a disposizione m sistemi di riferimento, andiamo ad applicarequesti ultimi al ligando iniziale, ottenendo appunto m ligandi base.

Per ogni sistema di riferimento Ri con i = 1, . . . ,m, definiamo le matriciMi ∈ R3 in modo tale che la prima riga corrisponda al versore a di Ri, laseconda e la terza riga rispettivamente b e c di Ri.

Mi =

ai

bi

ci

=

axi ayi azibxi byi bzicxi cyi czi

(3.4)

Chiamiamo queste matrici Mi, al variare di i tra 1 e m, matrice associataalla base ortonormale per Ri.


Il passo successivo consiste nell’applicazione del ligando su ogni sistemadi riferimento, con lo scopo di ottenerne un insieme da usare come input perligand-rotation. Questo processo puo essere descritto grazie alle matrici dellebasi ortonormali definite poco sopra dalla formula 3.4.

Sia L cosı definito

L = qj : j = 1, . . . t

l’insieme dei punti di R3, con |L| = t, tale che qj rappresenta la posizionenello spazio del j-esimo atomo del ligando di partenza. Quel che vogliamoottenere e un insieme di ligandi LR1 , . . . ,LRm tali che:

LRi =p(i)j : j = 1, . . . , t

i = 0, . . . ,m

con

p(i)j = Mi · qj j = 1, . . . , t

xp(i)j

yp(i)j

zp(i)j

=

axi bxi cxi

ayi byi cyiazi bzi czi

xqjyqjzqj

j = 1, . . . , t

i = 1, . . . ,m

(3.5)

dove le Mi rappresentano le matrici definite nella formula 3.4.

3.5.3 Campionamento di S

L’ultimo passo che resta da fare prima di far parire l’algoritmo vero e proprioe quello di traslare ogni ligando LRi su ogni punto di S in modo che D0

2

coincida con l’origine dei sistemi di riferimento; in questo modo partiranno|S| ·m istanze dell’algoritmo generatore di tutte le combinazioni valide.

L’insieme S e costruito dai dati forniti dall’utente: due punti e un passol. Si costruisce un cubo a partire dalle 2 coordinate e lo si partiziona in tanticubetti tali che il loro lato sia uguale a l e che il lato del cubo sia divisibileper l. L’insieme di tutti gli spigoli dei cubetti e proprio S.

Ligand-rotation non e un software statico, bensı dinamico e completamen-te configurabile dall’utente. Per quanto riguarda la costruzione dell’insiemeS, si dovranno specificare gli estremi dello spazio in cui si vuole che abbiainizio la costruzione del ligando e il passo di campionamento di questi punti;mentre, per quel che riguarda i sistemi di riferimento, le costanti k e n sonopersonalizzabili.

2D e l’elenco degli atomi del ligando; sono in ordine secondo la visita effettuata su G.Si faccia riferimento alla sezione 3.2

Capitolo 4

Dettagli implementativi

In questo capitolo verranno descritti gli algoritmi implementati e le strut-ture dati utilizzate in ligand-rotation. In particolare si mostrera (nella se-zione 4.1) come, a partire da un file di configurazione, si caricheranno tuttele informazioni relative al ligando ed alla proteina. Non soltanto, andremoanche a specificare tutte le opzioni disponibili per permettere una buonapersonalizzazione.

In aggiunta, andremo a definire (nella sezione 4.2) delle strutture dati chepermetteranno a ligand-rotation un incremento delle prestazioni, per quan-to riguarda la verifica della consistenza del vincolo (sezione 4.4). Questo eun aspetto molto importante, che ci permette di notare che la procedurache verifica i vincoli viene chiamata dopo ogni movimento del ligando (ossiaparecchie volte) e non e accettabile una inefficienza a questo riguardo.

Andremo ad esaminare nei dettagli (nella sezione 4.3) l’algoritmo prin-cipale del pacchetto. Descriveremo come le varie tecniche presentate fino adora interagiscono tra loro a partire dall’esplorazione del sito attivo fino allagenerazione di una conformazione consistente.

Come ultimo argomento (nella sezione 4.5) daremo alcuni accenni circala complessita di ligand-rotation, in termini di tempo e di spazio occupato.

4.1 Inizializzazione di ligand-rotation

In questa sezione verranno descritte le fasi iniziali in cui ligand-rotation vieneconfigurato, a partire dal file di input di configurazione, e inizializzate le suestrutture dati.

37

38 CAPITOLO 4. DETTAGLI IMPLEMENTATIVI

4.1.1 Acquisizione informazioni di configurazione

L’unico modo che ha un utente di comunicare con ligand-rotation e attraversoil file di configurazione; configurazione, in cui e possibile specificare tutti iparametri, sia funzionali che qualitativi, utilizzati dall’algoritmo di docking:

• ligando: specificare file mol21 contenente le coordinate del candidatoligando

• proteina: file mol2 contenente le posizioni degli atomi della proteinatarget

• step rotazioni: specifica il numero di rotazioni che un legame ruotantecompie

• legami bloccati: specifica, in sequenza, l’id2 dei legami che non devonoruotare

• coordinate sito attivo: specifica gli estremi dell’area in cui si trova ilsito attivo della proteina

• passo area di lavoro: indica il passo con cui l’area del sito attivo verrapartizionata

• set punti iniziali: specifica gli estremi dell’area in cui si piazzera il primoatomo del ligando

• passo punti iniziali: passo di partizionamento dell’area contenente ipunti iniziali

• punti sfera: numero dei punti sulla sfera di raggio unitario nei quali sigenereranno le basi ortonormali

• rotazioni su Z: step di rotazioni delle basi ortonormali sull’asse Z

Una volta in possesso di questi parametri, si andra a generare prima ditutto il grafo rappresentante il ligando, L = (N,A) e, in un secondo momento,l’albero T , grazie ad una visita V su L.

1mol2 rappresenta un formato standard di rappresentazione delle molecole.Approfondiremo l’argomento nella sezione 4.1.2

2L’Id corrisponde a quello specificato nel file mol2

4.1. INIZIALIZZAZIONE DI LIGAND-ROTATION 39

Implementazione

La realizzazione dei nodi e degli archi e stata resa possibile grazie all’ausiliodelle classi Nodo e Arco (nelle figure 4.1 e 4.2 viene presentata una possibileimplementazione delle due strutture). Il grafo che rappresenta la molecola estato implementato con delle liste di adiacenza (codice 4.3), con particolareriguardo all’efficienza, visti i numerosi accessi in lettura.

Codice 4.1: Struttura della classe Nodoclass Nodo

public :// . . .

private :int id atomo ;s t r i n g nome atomo ;MyReal x ;MyReal y ;MyReal z ;s t r i n g t ipo atomo ;MyReal ragg io atomico ;

; Codice 4.2: Struttura della classe Arco

class Arco public :

// . . .private :

int id l egame ;int atomo1 ;int atomo2 ;s t r i n g t ipo l egame ;int r o t a z i o n i l e g a m e ;

; Per motivi di efficienza, l’implementazione del grafo e stata realizzata per

mezzo di una matrice di vettori della libreria standard del C++. Nella primacolonna sono stati disposti tutti gli id dei nodi, mentre, su ogni riga, unalista degli id corrispondenti agli archi caratterizzanti il nodo stesso.


Codice 4.3: Struttura della classe Grafoclass Grafo

public :// . . .

private :s td : : vector<Arco> l egami ;std : : vector<Nodo> atomi ;std : : vector< std : : vector<int> > moleco la ;

// . . . ; Il passo successivo consiste nell’impostare, negli oggetti di tipo Arco, le

rotazioni possibili per ogni legame attraverso il metodo imposta legami(),specificate nel file mol2 con la seguente sintassi:

• ar: legame facente parte di un composto aromatico

• 1: legame singolo o semplice

• 2: legame doppio

Qui di seguito lo pseudocodice della procedura:

Codice 4.4: Pseudocodice del metodo imposta rotazioni()for arco in a r ch i m o l e c o l a ( )

i f arco . t ipo l egame i s ‘ ‘ ar ’ ’arco . r o t a z i o n i l e g a m e = 1

else i f arco . t ipo l egame i s not ‘ ‘ 1 ’ ’arco . r o t a z i o n i l e g a m e = 1

else i f arco . conta in s (H)arco . r o t a z i o n i l e g a m e = 1

elsearco . r o t a z i o n i l e g a m e = k

4.1.2 Struttura file mol2

Il file mol2 e uno dei sistemi per poter rappresentare la struttura di unamolecola attraverso file. Qui di seguito un esempio di file mol2:

4.2. STRUTTURE DATI DEL RISOLUTORE 41

@<TRIPOS>MOLECULE

4OHT

29 31 1 0 0

SMALL

USER_CHARGES

@<TRIPOS>ATOM

1 C1 31.1300 -1.5780 24.9020 C.ar 1 LIG600 0.00

2 O1 25.1930 -1.6480 26.9150 O.3 1 LIG600 0.00

3 C2 31.5550 -1.7460 26.2270 C.ar 1 LIG600 0.00

4 O2 33.4760 -3.6110 23.0430 O.3 1 LIG600 0.00

5 N1 34.3180 -4.9410 20.4290 N.3 1 LIG600 1.00

@<TRIPOS>BOND

1 1 2 ar

2 2 3 1

3 3 4 1

4 4 5 2

@<TRIPOS>SUBSTRUCTURE

1 LIG600 1 GROUP 1 LIG 1

Nella sezione MOLECULE vengono riassunte le informazioni principalidella molecola, come il suo nome, il numero degli atomi e dei legami chimici.Le sessioni ATOM e BOND sono le piu importanti per i nostri scopi. Nellaprima, infatti, vengono elencati tutti gli atomi della molecola, identificati daun Id, con le loro posizioni nello spazio e il tipo di ibridizzazione assunta. Peresempio, l’atomo indentificato dall’Id 2 e un ossigeno con coordinate x, y, zrispettivamente (25.1930,−1.6480, 26.9150), caratterizzato da una ibridizza-zione sp3. Nella seconda, si elencano i legami chimici esistenti tra atomi, incui la prima colonna contiene l’id del legame, la seconda e la terza, rispet-tivamente gli Id degli atomi coinvolti e la quarta il tipo di legame (ar perlegame facente parte di un composto aromatico, 1 legame singolo, etc..).

Per ulteriori approfondimenti si rimanda al sito webhttp://tripos.com/data/support/mol2.pdf

4.2 Strutture dati del risolutore

Qui di seguito verranno illustrate alcune strutture dati utilizzate dal risolu-tore dei vincoli di ligand-rotation.

http://tripos.com/data/support/mol2.pdf


4.2.1 Matrice di Vicinanza

La matrice di vicinanza, indicata conM∈ 0, 1n×n, n = |N | conN l’insiemedegli atomi del ligando, e una matrice quadrata di ordine n, booleana, dovegli indici di riga e di colonna corrispondono agli id degli atomi. Lo scopodi M e quello di specificare gli atomi adiacenti per ogni atomo del ligando.Nel momento in cui si posizionera un atomo a ∈ V e si andra a verificare laconsistenza del vincolo di distanza su di esso, si evitera di calcolare la distanzaeuclidea tra a e i suoi adiacenti. Questo perche, essendo i vicini una catenadi atomi legati tra loro covalentemente, non vale piu il concetto di sfera diVan der Waals. La lunghezza del legame covalente e inferiore rispetto allalunghezza di un legame di Van der Waals e, evitando questo controllo, siotterrebbero dei ‘falsi positivi ’.

Formalizziamo il concetto di adiacenza:

Definizione 4.1. Un atomo x ∈ N e adiacente ad un altro atomo y ∈ N see solo se il cammino tra x e y ha lunghezza massima 3, ovvero il camminonon deve essere composto da piu di tre legami.

Avremo quindi, per ogni atomo del ligando una serie di atomi adicenti.Indichiamo questo fatto con la seguente notazione: per ogni atomo a ∈ N

adj(a) = v ∈ N :M(a, v) =M(v, a) =M(a, a) = trueLa costruzione di M viene effettuata grazie a n visite in ampiezza sul

grafo L del ligando, dove |N | = |n|, specificando su ogni computazione illivello di profondita, cioe 3.

4.2.2 Discretizzazione di un sottoinsieme limitato diR3

Consideriamo una porzione cubica di spazio, V , identificato da due punti(due vertici del cubo) corrispondenti rispettivamente alle coordinate minimee massime delle tre dimensioni: lunghezza, larghezza e profondita. Siano Mine Max tali punti e sia p il lato del cubo.

Consideriamo una sua discretizzazione: dividiamo V in tante celle cubicheC di uguale dimensione, caratterizzate anch’esse da un minimo, MinC , unmassimo, MaxC , e un passo l ≤ p tale che p sia divisibile per l. Sia C l’insiemedi tutte le celle di V ; ogni cella rappresenta anch’essa un sottoinsieme limitatodi R3.

Ogni cella potra memorizzare al suo interno un insieme di atomi. Dato unatomo a, esso stazionera all’interno di una particolare cella Ca se le coordinatedi a appartengono all’insieme dei punti di R3, individuati da Ca.


Definizione 4.2. [Appartenenza]Siano (xa, ya, za) le coordinate di a e sia PCa l’insieme dei punti di R3 appar-tenenti alla cella Ca:

PCa =

(x, y, z) ∈ R3 : Minx ≤ x < Maxx ∧Miny ≤ y < Maxy ∧Minz ≤ z < Minz

dove (Minx,Miny,Minz) = MinC e (Maxx,Maxy,Maxz) = MaxC . Allora:

a ∈ Ca ⇐⇒ (xa, ya, za) ∈ PCa

Figura 4.1: Quattro celle cubiche con all’interno alcuni atomi.

L’insieme V appena introdotto avra un ruolo fondamentale nel risolutoredi vincoli da noi implementato. Le celle rappresentano la struttura base acui si accedera per la verifica della consistenza del vincolo di non sovrappo-sizione. Durante l’esecuzione dell’algoritmo di ligand-rotation, quando siamosicuri che un atomo, supponiamo Dk, non verra puo ruotato3, e necessarioverificare la consistenza su di esso. Per far cio, ci cerca di inserirlo all’internodi una cella; se nelle sue vicinanze ci sono altri atomi, l’inserimento fallisce,altrimenti ha successo. Viceversa, nel momento in cui si fa backtracking suDk, l’operazione da eseguire e quella inversa: eliminare l’atomo Dk dalla cellain cui si trova, in quanto si sta cercando una nuova configurazione e quindi su-bira una potenziale rotazione. In questo modo, la struttura dati conterra soloatomi consistenti e la verifica dei vicini sara effettuata molto efficientementegrazie ai metodi implementati su di essa.

3Nella sezione 4.3.2 descriveremo nei dettagli il momento esatto in cui si verifica laconsistenza.


Identificazione delle celle

Supponendo di dover inserire un atomo a all’interno di una cella, la suaidentificazione si ottiene ricavando tre indici, univoci per ogni cella, (X, Y, Z),a partire dalle coordinate (x, y, z) di a:

X = bx−Minx

lc Y = by −Miny

lc Z = bz −Minz

lc

Implementazione

Una simile struttura dati puo essere facilmente implementata con dei vectordella libreria standard del C++.

Codice 4.5: Implementazione di Vstd : : vector< std : : vector< std : : vector<Cel la> > >

Al contrario, una cella e stata implementata per mezzo di una classe,affinche possa memorizzare al suo interno l’insieme degli atomi.

Codice 4.6: Implementazione di una cellaclass Ce l l a

private :s td : : vector<Atomo> atomi ;

public :. . .

;

Insieme delle operazioni

L’insieme delle operazioni ammesse per la struttura V e per la Cella sonoelencate qui di seguito con il relativo pseudocodice.

• dato un atomo, identificazione della sua cella (codice 4.7)

• inserimento di un atomo in una cella (codice 4.8)

• eliminazione di un atomo da una cella (codice 4.9)


Codice 4.7: Pseudocodice per l’identificazione di una cellaV : : i d e n t i f i c a c e l l a ( atomo )

// . . .i d c e l l a [ 0 ] = ( atomo . x − Min . x )/p ;i d c e l l a [ 1 ] = ( atomo . y − Min . y )/p ;i d c e l l a [ 2 ] = ( atomo . z − Min . z )/p ;return i d c e l l a ;

Codice 4.8: Pseudocodice per l’inserimento di un atomo in una cella

V : : i n s e r i s c i ( atomo ) // . . .i d c e l l a = V : : i d e n t i f i c a c e l l a ( atomo ) ;i f ( i d c e l l a != NULL)

Ce l l a : : i n s e r i s c i ( i d c e l l a , atomo ) ;Ce l l a : : i n s e r i s c i ( i d c e l l a , atomo )

// . . .cubo . at ( i d c e l l a ) . push back ( atomo ) ;

Osserviamo che, durante l’identificazione della cella, nel caso in cui le

coordinate dell’atomo non appartengano alla porzione di spazio V riene re-stituito il valore NULL.

4.2.3 Celle contenenti atomi ‘vicini’

Lo scopo di questo paragrafo e illustrare in che modo vengono recuperatele celle contenenti gli atomi definiti ‘vicini’, ovvero quelli che subiranno ilcontrollo della distanza euclidea con l’atomo da inserire da parte del risolutoredi vincoli. Definiamo prima di tutto il concetto di atomo vicino e l’insiemeformato da questi ultimi.

Definizione 4.3. Sia a ∈ C un generico atomo appartente ad una cella,anch’essa generica. L’insieme G e composto dall’insieme degli atomi appar-tenenti ad una qualche cella di V tali che la loro distanza euclidea da a siainferiore o uguale alla somma del raggio atomico di a e del raggio atomicomassimale degli atomi del ligando.


Codice 4.9: Pseudocodice per l’eliminazione di un atomo da una cellaV : : e l im ina ( atomo )

// . . .i d c e l l a = V : : i d e n t i f i c a c e l l a ( atomo ) ;Ce l l a : : e l im ina ( i d c e l l a , atomo ) ;

Ce l l a : : i n s e r i s c i ( i d c e l l a , atomo )

// . . .cubo . at ( i d c e l l a ) . pop back ( atomo ) ;

Costruire l’insieme G non e affatto semplice, in quanto occorre esplorare

tutta la porzione di spazio che comprende l’atomo da inserire. L’approccioutilizzato da ligand-rotation consiste nell’individuazione di tutte le celle chepotenzialmente potrebbero contenere atomi vicini. Sia R il raggio massimodegli atomi componenti il complesso ligando-proteina, sia r il raggio dell’a-tomo da inserire e c la cella in cui inserirlo. Da queste due informazioni cicostruiamo il range delle celle da controllare, dato da:

range = bR + r

lc+ 1

L’insieme delle celle K i cui indici rispetto c spaziano della quantita rangesono proprio quelle che il risolutore andra ad esaminare per la verifica delvincolo di vicinanza.

4.3 Algoritmo di ligand-rotation

In questa sezione verra presentato nei dettagli il cuore dell’algoritmo di doc-king. Si parlera di consistenza in riferimento all’unico vincolo implementatosu ligand-rotation, ovvero quello di non sovrapposizione. La prima operazio-ne effettuata su V e quella relativa al piazzamento degli atomi della proteina.Questi, infatti, vengono inseriti immediatamente per il semplice motivo chela proteina e un corpo statico e quindi gli atomi non verranno mai sposta-ti di posizione. Notiamo che l’inserimento avviene senza il controllo dellaconsistenza.

4.3. ALGORITMO DI LIGAND-ROTATION 47

4.3.1 Parte I: inizializzazione

Come descritto formalmente nella sezione 3.5, questa fase si occupa dellagenerazione degli m sistemi di riferimento su cui applicare il ligando base,ottenendo quindi m ligandi di partenza pronti per essere computati, ognunocon un orientamento differente. Di seguito lo pseudocodice:

Codice 4.10: Pseudocodice relativo alla generazione dei ligandi baseP = p u n t i s f e r a ( k )S = c a m p i o n a p u n t i i n i z i a l i ( j )for p in P

Z = ass e zB = base (Z)X = p r o d o t t o V e t t o r i a l e (Z , B)Y = p r o d o t t o V e t t o r i a l e (Z , X)norma1 (X)norma1 (Y)T = numeroRotazioniZ ( )for t in T

ruota su Z ( angolo )t r a s l o m o l e c o l a (O(X, Y, Z) )ruota mo leco la (X, Y, Z)for s in S

t r a s l o m o l e c o l a ( s )gene ra con fo rmaz ion i ( )

angolo += step La funzione punti sfera(k) genera k punti disposti uniformemente su una

sfera di raggio unitario e, su ognuno di questi, si costruira un sistema diriferimento (come spiegato formalmente nella sezione 3.5.1). Successivamentelo si ruotera di T step lungo l’asse Z, ottenendo il terzo grado di liberta.Sull’origine dei sistemi di riferimento appena generati, si andra a posizionarel’atomo D0 del ligando e si ruotera la molecola in modo tale da applicarlela base appena generata. I ligandi che si otterrano andranno riposizionatisu ogni punto di S, insieme di coordinate definite dall’utente e campionatedal metodo campiona punti iniziali(j) (con passo j) facendo, come al solito,coincidere l’atomo D0 dei ligandi con i vari punti di S. In conclusione, perognuno dei ligandi applicati sui sistemi di riferimento costruiti, si fara partirel’algoritmo che generera tutte le strutture consistenti, attraverso la chiamataalla funzione genera conformazioni().


4.3.2 Parte II: generazione conformazioni

In seguito all’esecuzione della prima parte dell’algoritmo, ora abbiamo adisposizione un set di ligandi LRi , i = 1, . . . ,m da poter utilizzare comepartenza per la generazione delle conformazioni consistenti.

Sia fissato k ∈ 1, . . . ,m e consideriamo LRk il ligando su cui eseguirela computazione.

Sia inoltre D l’elenco ordinato degli atomi di LRk ottenuto da una visitaV sul grafo G, di cui ne abbiamo discusso nella sezione 3.3.

Qui di seguito lo pseudocodice dell’algoritmo utilizzato in ligand-rotation:genera con fo rmaz ion i ( i )

i f D[ i ] i s rooti f v e r i f i c a v i n c o l o (D[ i ] )

i++genera con fo rmaz ion i ( i )backtrack ing ( i )i−−backtrack ing ( i )

elsetermina computazione ( )

elses a l v a s o t t o a l b e r o (D[ i ] )a n g o l o r o t a z i o n e = 0for r o t a z i o n e in r o t a z i o n i (D[ i ] )

i f v e r i f i c a v i n c o l o (D[ i ] )r u o t a s o t t o a l b e r o (D[ i ] , a n g o l o r o t a z i o n e )i++i f not D[ i ] i s l a s t

gene ra con fo rmaz ion i ( i )backtrack ing ( i )

elses c r i v i c o n f o r m a z i o n e ( )

i−−r i p r i s t i n a s o t t o a l b e r o (D[ i ] )a n g o l o r o t a z i o n e += step

else :a n g o l o r o t a z i o n e += stepr i p r i s t i n a s o t t o a l b e r o (D[ i ] )r u o t a s o t t o a l b e r o (D[ i ] , a n g o l o r o t a z i o n e )

4.3. ALGORITMO DI LIGAND-ROTATION 49

Il primo passo dell’algoritmo sara quello di testare, sull’atomo radice, cioeD0, il vincolo di consistenza attraverso verifica vincolo(D0). Nel caso di falli-mento la computazione fallira immediatamente e si procedera con un nuovoligando LRh

con h 6= k; altrimenti, si memorizzera la posizione dell’atomo D0

sulla sua cella corrispondente, (cf. paragrafo 4.2.2) e si procedera con l’atomosuccessivo D1 con la chiamata ricorsiva genera conformazioni(1).

Supponiamo ora di essere arrivati ad una situazione di questo tipo: abbia-mo appena verificato e convalidato il vincolo sull’atomo Di−1 e ruotato il sot-toalbero con radice il nodo corrispondente a Di−1. Gli atomi attualmente pre-senti in V , oltre quelli appartenenti alla proteina, sono I = D0, . . . , Di−1.Siamo in procinto di avviare la computazione per l’atomo Di con la chiamataricorsiva genera conformazioni(i).

Una volta chiamata, la procedura salva sottoalbero(Di) salva la posizionedi tutti gli atomi del sottoalbero con radice l’atomo corrispondente a Di

(approfondiremo il perche di questa scelta nel paragrafo 4.3.2). Tra le possibiliscelte di rotazione su Di, si prende la prima e si verifica, quindi, il vincolo diconsistenza su Di.

Osserviamo come la verifica di consistenza viene effettuata sull’unico ato-mo Di. In questo istante, si e certi che Di non subira ulteriori rotazioni e sevalido, lo rimarra fino a quando non verra rimosso dalla struttura dati di cuisi serve il risolutore, tramite la procedura di backtracking. Inoltre, essendogli atomi dell’insieme I gia validati e quindi presenti all’interno delle celle diV si e certi che questi sono rimasti consistenti all’avanzare dell’algoritmo.

Nel caso di consistenza di Di, si procede alla rotazione del suo sottoalbero(come spiegato nella sezione 3.4). Fatto cio, si passa al nodo successivo e,se non e l’ultimo, si richiama ricorsivamente la procedura proprio su Di+1;in caso contrario, avendo superato con successo tutti i vincoli su tutti gliatomi, siamo arrivati alla fine della molecola. Creiamo, quindi, il file mol2 delligando appena generato. Nel caso in cui il vincolo di consistenza sull’atomoDi fallisce, si procede ad effettuare il backtracking per portarci all’atomoprecedente, si ripristina la situazione precedente alla rotazione e quindi, sitermina la computazione del nodo Di con la sua successiva eliminazione dallastruttura dati del risolutore attraverso la procedura backtracking(i). Se Di−1puo effettuare ulteriori rotazioni, si prosegue per questa scelta, altrimentitermina l’esplorazione del sottoalbero con radice Di−1.

Nel momento in cui si esplorano tutte le scelte dei nodi, l’ultima opera-zione da fare e eliminare l’atomo radice dalla sua cella e procedere con unnuovo ligando base.


Ripristino situazione precedenete alla rotazione

L’aspetto chiave della seconda fase dell’algoritmo e senza dubbio la visitadell’albero di ricerca. Qui vengono intraprese tutte le scelte possibili dei no-di e, di conseguenza, le rotazioni effettuate sul ligando base sono numerose;l’effetto negativo e la presenza di errori numerici. Infatti, sapendo di averea che fare con una mole di calcoli su R, dobbiamo preoccuparci dell’aritme-tica macchina e quindi del risultante accumulo di errore nelle operazioni dirotazione.

Si puo ricorrere a due tipi di soluzione:

• rotazione inversa

• memorizzazione informazioni prima della rotazione

Nel caso della rotazione inversa, prima di intraprendere le diverse sceltedi rotazione di un nodo, si effettua una rotazione inversa per portarci al casobase, quindi ruotare nuovamente per proseguire lungo una scelta. Lo svan-taggio di questa tecnica e che, seppur limitatamente, andiamo incontro adaccumulo di errore: la rotazione inversa non e la funzione inversa della ro-tazione. Per cui, l’applicazione di essa non ci portera alla situazione iniziale,bensı ad una simile che si riperquotera sulle rotazioni di tutto il suo sottoal-bero. Nella seconda soluzione, implementata su ligand-rotation, il ripristino sibasa su una apposita struttura dati. Prima della rotazione di un sottoalbero,cioe prima di effettuare una scelta sul livello k, si salva la situazione correntedei nodi, ovvero le loro coordinate spaziali, per poi venire ripristinata quandola procedura del backtracking ci riporta al livello k per percorrere un’altrascelta. Al contrario della prima soluzione, l’unico errore numerico commessoe nel calcolo della singola rotazione il quale, non si riperquotera ne sulle rota-zioni successive dello stesso nodo ne sulle rotazioni dei nodi del sottoalbero.L’unico svantaggio e nell’occupazione di memoria: per ogni livello e per ogniscelta si dovra memorizzare la situazione del sottoalbero. Comunque, avendoa che fare con ligandi composti da poche decine di atomi, al massimo uncentinaio, le informazioni da salvare sono contenute.

4.4 Consistenza del vincolo

La definizione di consistenza del vincolo data in precedenza, seppure esatta,e estremamente inefficiente. Infatti, e inutile effettuare il controllo di sovrap-posizione su atomi distanti o addiritura su atomi non ancora piazzati. Cioche formalizzeremo ora e l’algoritmo implementato da ligand-rotation per laverifica della consistenza, grazie alle strutture dati definite nella sezione 4.2.

4.5. ACCENNI DI COMPLESSITA 51

4.4.1 Verifica dei vincoli

Il controllo dei vincoli durante l’esecuzione dell’algoritmo di ligand-rotationfunziona in questo modo: una volta effettuata la rotazione di un atomo e delrelativo sottoalbero, si inserisce il primo atomo all’interno di una cella. Siaa l’atomo da inserire con coordinate spaziali (xa, ya, za) ∈ R e di raggio r.L’insieme C degli atomi da verificare e definito come:

C = G r adj(a)

dove:

• G e l’insieme degli atomi vicini, (cf. definizione 4.3)

• adj(a) e l’insieme degli atomi adiacenti, (cf. definizione 4.1)

Il vincolo di non sovrapposizione viene dunque verificato per ogni z ∈ C:

S(a, z) =

true, se ‖a− z‖ ≥ raggio(a) + raggio(z)

false, altrimenti

Qui di seguito lo pseudocodice implementato in ligand-rotazion.v e r i f i c a v i n c o l o ( a , Ce l l e , Adiacent i )

i f a not in Vreturn f a l s e

for c e l l a in Ce l l efor atomo in atomi . c e l l a

i f atomo i s not in a d i a c e n t ii f d(a , atomo ) < a . r agg i o + atomo . ragg i o

return f a l s eelse

return t rue

4.5 Accenni di Complessita

Il concetto di complessita computazionale e strettamente legato ai diversigradi di liberta che ligand-rotation possiede. Il grado di liberta che piu ‘ral-lenta’ l’algoritmo e la rotazione sul legame: avendo a disposizione k scelte per


ogni legame, si deduce abbastanza velocemente che ligand-rotation ha unacomplessita esponenziale sul numero dei legami n. In particolare, risulta:

CT (n) = O(kn)

dove k rappresenta il numero degli step effettuati per ogni legame ruota-bile.

Sebbene ligand-rotation cresca esponenzialmente, le tecniche della pro-grammazione a vincoli permettono di ridurre notevolmente la complessitaoperando dei tagli all’albero di ricerca. Infatti, intercettando l’inconsistenzadi un nodo, viene eliminato tutto il sottoalbero relativo e quindi parte dellesoluzioni totali, in modo tale da evitare la visita fino alle foglie dei percorsiinconsistenti.

L’altro grado di liberta e dato dal possibile posizionamento del ligandonelle tre dimensioni. Esso rappresenta una semplice costante moltiplicativaa CT , data dal numero di sistemi di riferimento totali e dal numero di puntiiniziali. Tuttavia, questa costante e abbastanza alta4 e, non avendo un numerodi legami ruotanti estremamente alto, il passo di discretizzazione puo influiresul costo complessivo dell’algoritmo.

4Nel caso di 210 diversi sistemi di riferimento e 27 punti iniziali, la costante e pari a5670, quindi da non sottovalutare.

Capitolo 5

Risultati

In questo capitolo si illustreranno i risultati ottenuti da ligand-rotation. De-scriveremo (nella sezione 5.1) il complesso preso in esame per la realizzazionedei test. Successivamente (nella sezione 5.2) analizzaremo in dettaglio i testeseguiti, mostrando il comportamento registrato da ligand-rotation al varia-re dei parametri di qualita. Presenteremo, in particolare, una configurazioneottimale. Infine (nella sezione 5.3) faremo dei confronti sulle possibili visiteimplementabili per generare l’albero e su come esse andranno ad incidere sultempo di esecuzione.

Come primo test di fattibilita non consideriamo una funzione energeticabensı confrontiamo i risultati ottenuti con il ligando minimizzato preso daPDB attraverso un test di tipo geometrico, calcolando il Root Mean SquareDeviation (RMSD) definito come segue:

Definizione 5.1. Siano L e C gli insiemi degli atomi rispettivamente del li-gando minimizzato e di quello appena generato. Sia N = |L| = |C|; definiamol’RMSD come:

RMSD =

√√√√ 1

N

N∑i=1

‖ni −mi‖ ni ∈ L,mi ∈ C (5.1)

Valori dell’RMSD pari a 0.5A stanno a significare che la distanza qua-dratica media tra gli atomi del ligando minimizzato e quello ottenuto daligand-rotation e di 0.5A. Le distanze di legame sono circa di 1.5A; com-mettere un errore sotto l’Angstrom caratterizza le posizioni degli atomi consufficiente precisione.

53

54 CAPITOLO 5. RISULTATI

5.1 Applicazioni possibili

Figura 5.1: Struttura molecolare del tamoxifene

Il complesso ligando-proteina studiato in questa tesi ha code number (re-lativo a PDB) 1XPC ; esso possiede numerosi ligandi, ma il nostro studio estato rivolto principalmente all’analisi del tamoxifene.

Sono stati effettuati i primi test di ligand-rotation con il tamoxifene(OHT). In figura 5.1 e illustrata la struttura del farmaco, mentre in 5.2 eillustrato il complesso ligando-proteina. Questo e un antitumorale apparte-nente alla famiglia dei farmaci interagenti con il recettore degli estrogeni eviene utilizzato nei casi di tumore mammario metastatico con buoni risultati.Successivamente sono stati scoperti gli ottimi benefici nella prevenzione dellaripresa della malattia in donne gia operate per tumore al seno.

Il tamoxifene presenta una struttura molecolare relativamente semplice;altri possibili ligandi non sono stati presi in considerazione, in quanto laloro struttura presentava dei cicli pentano e cicli esano che complicavanopesantemente l’algoritmo (cf. sezione 6.1). Avendo svolto questo lavoro intempi relativamente brevi, non si e riusciti nell’intento di gestire molecolecon le caratteristiche sopra citate.

5.2 Efficienza

Considerando il fatto che ligand-rotation e ancora un programma in evo-luzione e che necessita di notevoli estensioni per essere utilizzato in ambi-to professionale, i risultati ottenuti sono abbastanza soddisfacenti. Il nostroobiettivo e quello di generare conformazioni con un RMSD minore di 1A intempi ragionevoli. Come vedremo, e stato raggiunto con successo. La figura5.3 mostra uno tra i migliori risultati ottenuti con ligand-rotation.

5.2. EFFICIENZA 55

Figura 5.2: Immagine del complesso esaminato ottenuta effettuando un ren-dering con Molegro Molecular Viewer: in rosso e possibile vedere la strutturadella proteina, mentre in verde il ligando OHT


Figura 5.3: In verde e visibile il ligando ottenuto con ligand-rotation mentrein rosso e rappresentato il ligando minimizzato, scaricato da PDB. Il valoredi RMSD tra le due molecole e pari a 0.799538, ottenuto con la seguenteconfigurazione: 7 step di rotazioni, 210 sistemi di riferimento, passo cella di 1Ae punti iniziali campionati a 0.5A. Il tempo di esecuzione di ligand-rotation,con questa specifica configurazione, e stato di 172 minuti.

5.2. EFFICIENZA 57

5.2.1 Variazione parametri qualitativi

Questa classe di test ha l’obiettivo di stabilire il grado di efficienza di ligand-rotation al variare dei parametri di configurazione; si terra conto del valoredell’RMSD migliore e del tempo impiegato nel generare tutte le conforma-zioni.

Variazione degli step di rotazione

Esaminiamo in dettaglio il comportamento registrato dal pacchetto di doc-king a variazioni degli step di rotazione. Ricordiamo che lo step di rotazionerappresenta il valore delle variabili di tipo legame.

Lo scopo e quello di misurare l’efficienza e il corretto posizionamento delligando al variare degli step. Fissatone uno, i test verranno condotti su confi-gurazioni nelle quali variera il numero dei sistemi di riferimento. Il numero distep di rotazione presi in considerazione sono 2, 3, 4, 5, 6, 7. Come visibilenel grafico 5.4, l’esecuzione con un numero di step inferiore a 5, sebbene im-pieghi molto meno, in particolare non si superano i 5 minuti di calcolo, nonproduce risultati significativi sotto l’Angstrom. Dal grafico si puo notare ilmiglioramento delle conformazioni ottenute, in termini di RMSD, all’aumen-tare degli step di rotazione, a discapito del tempo impiegato, il quale restacomunque sempre ragionevole, ovvero entro le 6 ore di calcolo.

Variazione campionamento delle basi ortonormali

Questa classe di test e volta a misurare la reazione di ligand-rotation al variaredei sistemi di riferimento. Dopo averne fissato uno, i test eseguiti varierannoper il numero di step di rotazione. Nel grafico 5.5 sono visibili i valori ottenuticampionando le basi ortonormali: sull’asse delle ascisse sono riportati i tempiin scala logaritmica, mentre sull’asse delle ordinate i valori di RMSD miglioriper ogni istanza di ligand-rotation. I risultati, limitatamente al ligando testa-to, mostrano che anche con soli 541 sistemi di riferimento ligand-rotation siastato in grado di ottenere strutture molecolari molto valide, entro l’Angstromdi RMSD, in tempi molto brevi, precisamente in un’ora e mezza di calcoliper una conformazione con RMSD pari a 0,85534A.

Campionamento delle celle di V

Il seguente test e rivolto a identificare la dimensione ottimale delle celle diV . Nel grafico 5.6, l’asse delle ascisse individua il numero degli step di rota-zione rispetto ai quali si sono effettuati i test, mentre l’asse delle ordinate, in

1Esattamente 9 punti sulla sfera e 6 step di rotazione sull’asse Z


10−1

100

101

102

103

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

tempo espresso in minuti

rmsd

Prestazioni−qualità al variare degli step di rotazione

7 step6 step5 step4 step

Figura 5.4: Grafico rappresentante qualita in termini di RMSD, rispetto alligando minimizzato, e prestazioni, in termini di tempo di esecuzione delprogramma. Le molecole accettabili si trovano sotto l’Angstrom.

5.2. EFFICIENZA 59

10−3

10−2

10−1

100

101

102

103

0.5

1

1.5

2

2.5

3

3.5

tempo espresso in minuti

rmsd

Prestazioni−qualità al variare del passo di campionamentodella base ortonormale

2101207254

Figura 5.5: Grafico rappresentante qualita in termini di RMSD, rispetto alligando minimizzato, e prestazioni, in termini di tempo di esecuzione delprogramma. Le molecole accettabili si trovano sotto l’Angstrom.


scala logaritmica, identifica il tempo in minuti dell’esecuzione di ogni istanzadi ligand-rotation. La dimensione del lato delle celle e un parametro moltoimportante, che andra ad influire sul tempo totale della computazione manon sulla qualita dei risultati. Questo perche, nel momento della verifica deivincoli, si andranno a selezionare e controllare tutte le celle che potenzial-mente possono contenere atomi sovrapposti. Una dimensione troppo piccolacomporterebbe un numero di celle da visitare elevato, con conseguente per-dita di tempo; su molte di esse potrebbe non esserci nemmeno un atomo.Al contrario, una dimensione troppo grande comporterebbe una verifica deivincoli anche su atomi ‘lontani’ poiche una cella potrebbe contenere atomidistanti tra loro, rappresentando una porzione di spazio relativamente gran-de. Dai risultati, si evince che la misura ottimale e data da 2A: considerandol’esecuzione con 6 step, il tempo impiegato da quest’ultima rispetto ad uncampionamento a 0.5A si dimezza.

5.3 Confronti su visite

Un aspetto molto importante e la generazione dell’albero a partire dal grafo,dato che sta alla base della visita sullo spazio di ricerca delle combinazioni delligando. Modificando l’albero si influenzera l’ordine di esplorazione delle va-riabili (e di conseguenza l’albero di ricerca) e si potra avere un miglioramentoo un peggioramento delle prestazioni, in modo dipendente a seconda dallaparticolare posizione assunta degli atomi nella lista D e da quanto spessofalliranno.

5.3.1 DFS sul grafo della molecola

Il primo algoritmo, implementato su ligand-rotation con il compito di visitareil grafo G = (V,E) rappresentante il ligando, e la DFS, visita in profondita.La strategia di ricerca esplora il grafo andando, in ogni istante dell’esecuzio-ne dell’algoritmo, il piu possibile in profondita: gli archi del grafo vengonoesplorati a partire dall’ultimo vertice scoperto v che abbia ancora degli archinon esplorati uscenti da esso. Una volta terminata l’esplorazione di tutti gliarchi non esplorati del vertice v, si ritorna indietro per esplorare tutti gli ar-chi uscenti a partire dal vertice da cui v′ era stato precedentemente scoperto.Il processo di esplorazione continua fin quando tutti i vertici del grafo nonsiano stati esplorati.

Grazie al fatto che G e connesso, ovvero che la cardinalita delle suecomponenti connesse e 1, la ricerca DFS produce un unico albero.

5.3. CONFRONTI SU VISITE 61

4 5 610

0

101

102

103

104

105

step rotazioni

tem

po

esp

ress

o in

min

uti

Prestazioni al variare del campionamento delle celle

0.5Å1Å2Å3Å

Figura 5.6: Grafico rappresentante il tempo di esecuzione al variare del cam-pionamento delle celle di V . Raggruppati sulle ascisse, i test svolti con 4,5, 6 step di rotazione, mentre sull’asse delle ordinate il tempo impiegato daligand-rotation (in scala logaritmica).


La ricerca DFS, inoltre, marca ogni vertice con ben precise informazionitemporali, in particolar modo aggiorna due etichette per ogni vertice: d[v]che registra quando il generico vertice v e stato scoperto ed f [v] che registraquando e stata esplorata l’intera lista di adiacenza di v.

La complessita computazionale di una DFS e pari a Ω(V+E). In ogni caso,rispetto al tempo di esecuzione dell’intero programma, e assai irrilevante: laDFS fa parte della fase di pre processing di ligand-rotation, per cui vieneeffettuata una sola volta.

5.3.2 Modifica della DFS: first fail

Una piccola, ma sostanziale, modifica all’algoritmo della DFS consiste nell’or-dinare tutti i nodi adiacenti a quello attualmente visitato in base al numerodi rotazioni che possono effettuare: dal piu basso al piu alto. La strategia chesi vuole implementare e quella simile ad una first fail, ovvero si fa in modoche durante l’iterazione della lista D si fallisca il prima possibile. Cosı fa-cendo, poiche interviene il backtracking per portarci allo stato precedente, sielimina tutto il sottoalbero relativo al vertice in cui abbiamo fallito, andandoa ridurre il numero di soluzioni finali e risparmiando il tempo impiegato aduna esplorazione di un sottoalbero privo di soluzioni.

L’implementazione su ligand-rotation, in pratica, e una first fail statica.La scelta delle variabili viene effettuata nella fase di preprocessing, durante laDFS, e l’ordinamento riguarda gli archi del singolo nodo, indipendentementedagli altri.

La modifica della DFS puo trovare giustificazione nel fatto che avendo daattuare due scelte, una dipendente dall’altra, conviene intraprendere primaquella con meno possibilita. Cosicche, in caso di fallimento, il sottoalberoeliminato conterra piu nodi.

Nel grafico 5.7 sono rappresentati i risultati ottenuti sia da istanze diligand-rotation su cui e stata implementata la first fail che da altre senza. Daitest effettuati si puo concludere che la first fail ha avuto successo, riducendoi tempi di esecuzione.

5.3. CONFRONTI SU VISITE 63

4 5 610

0

101

102

103

104

step rotazioni

tem

po

esp

ress

o in

min

uti

Prestazioni al variare della visita sul grafo

DFS first failDFS

Figura 5.7: Grafico rappresentante il tempo di esecuzione di istanze di ligand-rotation con first fail e senza.

Capitolo 6

Conclusioni e Sviluppi futuri

In questo lavoro di Tesi si e affrontato il problema del docking. In particolare,lo scopo era volto alla costruzione di un nuovo pacchetto in grado di farinteragire tra loro due molecole, senza che cio fosse influenzato da funzionidi scoring che valutano la stabilita energetica.

La soluzione offerta a questo problema e, appunto, ligand-rotation; pro-gramma da me studiato e implementato, in grado di generare strutture mo-lecolari geometricamente corrette e, oltretutto, in maniera efficiente, per quelche riguarda il tempo impiegato nei calcoli.

Grazie all’approccio interamente geometrico, i test di validita eseguitisu ligand-rotation hanno dimostrato che e possibile raggiungere una qualitaaccettabile rispetto alle richieste reali da parte delle case farmaceutiche.

Con una configurazione standard, che porta il programma a conclusionenel giro di 6 ore al massimo1, si riesce a produrre posizionamenti di molecoleentro l’Angstrom di RMSD, cioe si ottengono ligandi dalla struttura moltosimile a quella minimizzata, scaricata dal database PDB.

Una passo successivo a questa Tesi consiste nell’applicare una funzionedi scoring per valutare l’energia complessiva del complesso ligando-proteinadeterminando la stabilita del sistema; in particolare, ci interessa quanto siacorrelata all’RMSD ottenuto dai nostri test.

Verranno ora mostrati alcuni riferimenti utili all’evoluzione di ligand-rotation. L’ingegneria del software insegna che un programma informaticoe in costante evoluzione, sia per quel che riguarda la risoluzione di bug esi-stenti che per l’aggiunta di nuove funzionalita. Essendo un progetto giovane,l’intento e quello di poter sviluppare nuove metodologie che constentiranno aligand-rotation l’elaborazione di differenti tipi di molecole, anche complesse,

1Il tempo impiegato da ligand-rotation dipende molto dal ligando scelto e dal numerodi legami ruotabili.

65

66 CAPITOLO 6. CONCLUSIONI E SVILUPPI FUTURI

garantendo tempi di esecuzione ragionevoli e piu precisione nella generazionedelle strutture dei ligandi.

Nella sezione 6.1 cercheremo di dare una possibile estensione all’algorit-mo esistente, in modo tale da poter trattare molecole con caratteristiche distruttura particolari.

Successivamente (nella sezione 6.2) si descrivera l’idea di un nuovo po-tenziale modello a vincoli che prevede l’inserimento di ulteriori vincoli el’implementazione della propagazione.

6.1 Molecole contenenti cicloesano

Un problema incontrato nello sviluppo di ligand-rotation e stato che, perparticolari ligandi, non si e riusciti ad ottenere un loro buon posizionamentoin termini di RMSD con il ligando minimizzato. Sono stati condotti numerositest e le migliori conformazioni si aggiravano intorno a 1.5A - contro il nostrolimite di accettazione di 1A. Il problema non stava nemmeno alla base delleopzioni di campionamento; infatti, all’aumentare dei parametri di qualita,non si apprezzavano miglioramenti da parte delle strutture generate.

Figura 6.1: A sinistra il ligando ottenuto con ligand-rotation, a destra illigando minimizzato scaricato da PDB

La strada da percorrere e completamente diversa. Alla base di cio ci sonoparticolari strutture molecolari denominate Cicloesano e Ciclopentano, chedovranno essere trattate in modo completamente differente.

6.1. MOLECOLE CONTENENTI CICLOESANO 67

6.1.1 Cicloesano

Il cicloesano fa parte della famiglia dei cicloalcani; e composto da sei atomidi carbonio e 12 atomi di idrogeno, e privo di doppi e tripli legami e la suastruttura e chiusa ad anello.

Gli atomi di carbonio hanno ibridizzazione sp3, andando a formare dellestrutture tetraedriche con angoli di legame di 109.5. La struttura quindinon sara planare, la quale ha angoli di 120, bensı tendera ad assumere unaconformazione a sedia, visibile in figura 6.2.

Figura 6.2: La figura mostra la struttura del cicloesano con a sinistra unaconformazione a sedia e a destra una conformazione a barca

Oltre quella a sedia, esistono altre conformazioni, sebbene con maggioreenergia e quindi meno stabili: conformazione a barca (in figura 6.2) e a mezzasedia. In figura 6.3 si esprime l’energia del cicloesano per ogni conformazione.

Figura 6.3: La figura rappresenta le possibili energie delle strutture delcicloesano in funzione alle conformazioni assunte


6.1.2 Nuovi gradi di liberta

Tutte le conformazioni del cicloesano elencate nella precedente sezione nonpossono essere ottenute a partire da semplici rotazioni. E necessario imple-mentare un nuovo algoritmo da inglobare a quello esistente, in grado di gestirepiu gradi di liberta. L’idea e quella di individuare queste strutture particola-ri attraverso la visita del grafo del ligando e, successivamente, nel momentoin cui si piazzano i suoi atomi, e necessario capire la conformazione attualeper andare a simulare tutti i suoi possibili stati conformazionali prima dellarotazione del sottoalbero.

6.2 Espansione modello a vincoli

Un primo ampliamento da realizzare su ligand-rotation potrebbe essere quellodi integrare il modello a vincoli, introducendono di ulteriori, come la verificadella lunghezza di un legame chimico, oppure implementare la propagazionedei vincoli. Nei paragrafi 6.2.1 e 6.2.2 svilupperemo questi due argomenti.

6.2.1 Controllo lunghezza legame

Con il modello a vincoli attualmente implementato in ligand-rotation, la ve-rifica della lunghezza dei legami chimici sarebbe completamente inutile, anzisarebbe un aggravio alla complessita totale del sistema. Al contrario, riuscen-do a trattare molecole come il cicloesano o il ciclopentano, si introdurrebbeun nuovo grado di liberta e il controllo diventerebbe sensato. Effettuandodelle rotazioni anche sui legami delle precedenti strutture si potrebbe andareincontro ad una rottura del ciclo.

6.2.2 Implementazione propagazione dei vincoli

Il modello a vincoli qui presentato non implementa su di se la propagazio-ne dei vincoli; questione non affrontata a causa dello scarso tempo avuto adisposizione. Dalla sezione 1.4 sappiamo che la propagazione ha il compi-to di ridurre lo spazio di ricerca, eliminando eventuali elementi del domi-nio. Nel nostro caso, l’eliminazione riguarda le possibili rotazioni. L’idea equella di evitare di visitare porzioni di albero inconsistenti. Nel contesto diligand-rotation, una possibile definizione di propagazione sarebbe questa:

nel momento in cui si identifica una rotazione che sicuramente portaa fallimento, la si elimina dal dominio della sua variabile legame. Comeconseguenza di cio, si potrebbe avere il risveglio a catena di altre propagazionie quindi l’eliminazione di altre rotazioni.

6.2. ESPANSIONE MODELLO A VINCOLI 69

L’interazione tra vincoli di non sovrapposizione e rotazione di un legamepossono fornire deduzioni sulla effettiva validita dei valori del dominio dellevariabili. Supponiamo di compiere un movimento alla struttura molecolare.Attraverso un algoritmo, possiamo essere in grado di stabilire a priori chesicuramente la rotazione appena effettuata non porta a nessuna soluzionevalida. In virtu di cio, si deduce che i valori ammissibili di alcune variabilidi tipo legame sono incongruenti con i vincoli implementati e si effettua,quindi, una riduzione del dominio delle variabili implicate. Osserviamo comeil costo computazionale degli algoritmi di propagazione debba essere moltobasso affinche il loro utilizzo sia giustificato.

Bibliografia

[1] Thomas Lengauer (Ed.). Bioinformatics - From Genomes to Drugs,Volume I. WILEY-VCH, 2001.

[2] Peter Atkins, Loretta Jones. Principi di chimica. Zanichelli, 2005.

[3] William K. Purves, David Sadava, Gordon H. Orians and H. CraigHeller. Elementi di Biologia e Genetica. Zanichelli, 2005.

[4] Harold Hart, David J. Hart, Leslie E. Craine. Chimica organica, quartaedizione. Zanichelli, 1998.

[5] Krzysztof R. Apt. Principles of Constraint Programming. CambridgeUniversity Press, 2003.

[6] Bjarne Stroustrup. C++: linguaggio, libreria standard, principi diprogrammazione. Addison-Wesley, 2000.

[7] Gregory L. Warren, C. Webster Andrews, Anna-Maria Capelli, BrianClarke, Judith LaLonde, Millard H. Lambert, Mika Lindvall, Neysa Ne-vins, Simon F. Semus, Stefan Senger, Giovanna Tedesco, Ian D. Wall,James M. Woolven, Catherine E. Peishoff and Martha S. Head. A criti-cal Assessment of Docking Programs and Scoring Function. Journal ofMedicinal Chemistry, 2006, Vol. 49, No. 20.

[8] http: // mathworld. wolfram. com/ RotationFormula. html .

71

http://mathworld.wolfram.com/RotationFormula.html

Studio e sviluppo di metodi computazionali per l’analisi ... · e le tecniche dettagliate...

Documents

Transcript of Studio e sviluppo di metodi computazionali per l’analisi ... · e le tecniche dettagliate...