Modellistica molecolare
Davide Cittaro e Marco Stefani
Master bioinformatica 2003
Introduzione
La genomica ha mostrato come la sequenza di aminoacidi di una proteina non sia
sufficiente per determinare con precisione la sua struttura, e con essa la sua
funzionalità: si è visto infatti che sequenze molto simili, se non addirittura identiche,
possono generare proteine con strutture e quindi funzionalità diverse.
D’altro canto, la capacità di predire la struttura di una proteina, data la sequenza, è
importante per essere in grado di analizzare da un punto di vista teorico le
funzionalità di proteine sconosciute. Lo scopo è studiare la capacità di interazione
della proteina con l’ambiente circostante, cioè con le altre proteine e il solvente
presenti nella cellula, ma anche come la struttura della proteina cambia in seguito ai
cambiamenti ambientali o alle mutazioni della sua struttura.
La teoria
Uno degli strumenti che cercano di risolvere questo problema è la modellistica
molecolare, che, utilizzando le conoscenze sull’interazione tra atomi e molecole che
derivano dalla chimica e dalla fisica, cerca di determinare la struttura finale della
molecola. Questa metodologia è detta ab initio proprio perché si basa sulle
caratteristiche fondamentali degli atomi e sulle leggi che determinano le interazioni
tra atomi.
La modellistica molecolare studia lo stato energetico delle molecole partendo dalle
relazioni geometriche e spaziali tra gli atomi che le compongono.
Una molecola viene descritta attraverso tre valori, che determinano univocamente le
relazioni geometriche e spaziali di qualunque numero di atomi:
•
•
•
distanza di legame: è la distanza che separa due atomi
angolo di valenza: è l’angolo individuato da tre atomi
angolo diedro: è l’angolo di rotazione tra i due atomi centrali in un sistema di
quattro atomi. Viene anche detto angolo di torsione
L’energia di una molecola è definita come una ipersuperficie detta Superficie di
Energia Potenziale (PES), che ha (3N - 6) + 1 dimensioni, dove N è il numero di
atomi. Le 6 dimensioni sottratte sono i gradi di libertà legati ai moti esterni; ciò rende
il sistema leggermente meno complesso, ma soprattutto lo rende indipendente dalla
posizione della molecola rispetto al sistema di riferimento. La PES è una curva
continua e quindi in teoria ammette infinite configurazioni; di queste, però, viene
cercato solo il sottoinsieme delle configurazioni chimicamente possibili, e di queste,
quelle chimicamente più stabili, cioè con derivata prima nulla, hanno un interesse
biologico perché sono le configurazioni che la molecola assume con maggiore
probabilità.
Da un punto di vista teorico, utilizzando le leggi della meccanica quantistica, è
relativamente semplice costruire un modello teorico in grado di descrivere con una
buona approssimazione il sistema di atomi e la sua evoluzione nel tempo. Purtroppo,
quando si abbandona la descrizione analitica per calcolare numericamente lo stato del
sistema in un certo istante, la dimensione del problema rende sostanzialmente
intrattabile il problema per una tipica proteina di qualche centinaio di aminoacidi. La
complessità computazionale è dell’ordine di O(N3). Nella realtà, questo metodo è
praticabile solo con molecole di qualche decina di atomi.
La prima semplificazione è data dall’utilizzo della meccanica classica: descrivere gli
atomi come entità indivisibili riduce sia il numero di variabili sia la complessità
intrinseca del sistema, che scende a circa O(N2). Le relazioni tra gli atomi sono
espresse dalla legge di Hooke: gli atomi si comportano come se fossero legati tra loro
da molle.
L’energia interna della molecola è data da due tipologie di parametri: i cosiddetti
termini di legame e i termini di non legame. La quantità calcolata è sempre
considerata relativa all’energia della molecola allo stato naturale, cioè allo stato di
minima energia potenziale in condizioni di equilibrio energetico. I termini di legami
sono i legati alle relazioni spaziali tra gli atomi; i termini di non legame descrivono le
interazioni coulombiane e di Van der Waals, determinate cioè dalla loro carica.
Poiché l’energia potenziale della molecola in uno stato di non equilibrio non è nota,
ma, in condizioni normali, la molecola vibra in un intorno dell’equilibrio, è
ragionevole approssimare la funzione dell’energia intorno al valore di minimo con
una espansione in serie di Taylor. Sfruttando una serie di considerazioni matematiche,
è possibile semplificare ulteriormente l’equazione.
La PES di una molecola presenta un numero estremamente elevato di punti stazionari,
molti dei quali hanno un’energia troppo elevata per essere chimicamente verosimili.
Le conformazioni interessanti, quindi, sono quelle che hanno la maggiore probabilità
reale di verificarsi e in cui la molecola tende a stabilizzarsi. Questa probabilità è
definita dalla probabilità di Boltzman:
pi =
e − ∆Ei / KT
,
∑ e −∆Ei / KT
i
dove ∆Ei = E0 - Ei e E0 è il minimo assoluto.
La probabilità di Boltzman è funzione della temperatura dell’intero sistema: alcune
conformazioni sono più probabili a una certa temperatura T. Ma aumentando la
temperatura, cresce anche l’energia cinetica della molecola, la quale può superare le
eventuali barriere energetiche e cadere in un altro minimo.
La ricerca dei punti stazionari può utilizzare metodi di calcolo locali o globali, che
differiscono essenzialmente nella capacità di superare le barriere potenziali della PES.
I metodi locali si basano tutti sulla variazione del gradiente: dato un punto specifico
della PES, si cercano le variazioni della geometria della molecola fino a raggiungere
il minimo locale che rappresenta l’attrattore del punto di partenza. Poiché la derivata
prima della funzione dell’energia potenziale non è nota, occorre utilizzare delle
approssimazioni che danno il valore della derivata in un insieme discreto di punti,
partendo dal valore noto della derivata prima del punto di partenza. Il metodo del
gradiente segue la direzione dello spostamento dell’energia in modo da minimizzare il
gradiente, fino ad annullarlo. Esistono due tecniche fondamentali, che nella ricerca
vengono usate in successione:
• metodo Newton-Raphson
• metodo lineare.
I metodi globali si muovono su tutta la PES alla ricerca del minimo globale. Poiché
l’esplorazione sistematica della superficie è normalmente impensabile, i metodi
campionano un sottoinsieme delle possibili conformazioni secondo strategie
deterministiche o stocastiche; le conformazioni più promettenti (in genere quelle con
energia minima) vengono analizzate con un metodo locale alla ricerca del punto
stazionario di attrazione. Dopo un sufficiente numero di campionamenti e dei
rispettivi minimi locali, si considera il minimo globale come il punto stazionario di
minore energia tra tutti quelli trovati.
L’unico metodo in grado di garantire il minimo globale è il cosiddetto simulated
annealing, che però in pratica non viene usato perché converge troppo lentamente.
Tra le tecniche alternative, che non garantiscono il risultato ottimale, ci sono il
metodo stocastico Metropolis-Montecarlo: da un valore iniziale dato (seed),
l’algoritmo genera un insieme casuale di numeri sulla base dei quali si variano le
coordinate interne della molecola ottenendo una nuova conformazione; se l’energia
associata è minore dell’energia della configurazione di partenza, la configurazione
viene accettata e minimizzata; se l’energia è maggiore, si calcola la probabilità di
Boltzman p associata alla variazione dell’energia tra le due configurazioni e si genera
un nuovo numero casuale r: la nuova configurazione viene accettata se p > r. Le
configurazioni accettate sono la base per una nuova iterazione, fino a quando si
decide di avere trovato il minimo globale.
I tempi di convergenza non sono garantiti e c’è la possibilità di calcolare più volte la
stessa configurazione.
La ricerca sistematica è un metodo deterministico, la cui velocità di esecuzione è
inversamente proporzionale alla risoluzione della PES, cioè alla dimensione degli
incrementi con cui si esplora la PES. D’altro canto, è evidente che una bassa
risoluzione porta ad un campionamento approssimativo che potrebbe fare perdere
molti punti stazionari.
Un altro metodo deterministico modifica la funzione dell’energia in modo tale che, in
un dato intervallo, abbia un solo minimo, che viene trovato iterativamente dividendo
l’intervallo in sottointervalli sempre più piccoli. Questo metodo garantisce di trovare
il minimo, cosa impossibile per i metodi stocastici.
La dinamica molecolare analizza l’evoluzione della conformazione molecolare nel
tempo, cioè cerca di determinare le traiettorie degli atomi della molecola sottoposti
alle forze determinate dall’energia cinetica del sistema e dalle interazioni tra gli atomi
stessi.
Per ogni atomo occorre quindi risolvere l’equazione fondamentale del moto:
ρ
ρ
dr
f i = mi 2i .
dt
L’intero sistema di equazioni differenziali viene calcolato ad intervalli regolari, che
dipendono dalla frequenza del moto che si vuole studiare. Le vibrazioni molecolari,
per esempio, che sono il movimento fondamentale della molecola, sono dell’ordine
del femtosecondo. Ciò implica, che se si vuole ottenere una ‘istantanea’ dello stato
della molecola, è indispensabile definire gli intervalli di tempo della simulazione
dello stesso ordine di grandezza.
Le condizioni del sistema all’interno del quale si conduce la dinamica molecolare
sono fondamentali per determinare l’andamento della simulazione. Se si impone che
il sistema sia isolato, le forze interne devono rimanere costanti per tutta la dinamica;
se inoltre l’inizio della dinamica corrisponde allo stato “naturale” della molecola,
allora le forze interne devono essere nulle, cioè i momenti angolare e lineare devono
essere conservati.
Due parametri fondamentali per la definizione del sistema sono l’energia e la
temperatura. Poiché ad ogni grado di libertà della molecola corrisponde un’energia
termica, l’energia cinetica del sistema è data da
∑K
i
i
=
3N − 6
k bT .
2
Mantenendo costante uno dei due valori, si ottengono due tipici sistemi:
• NVE o microcanonico
• NVT o canonico
Entrambi hanno impongono che il numero di molecole e il volume all’interno del
quale si muovono gli atomi siano costanti. Nel modello microcanonico l’energia
totale del sistema deve restare costante, in quello canonico la temperatura.
Gli esperimenti
Utilizzando le tecniche della modellistica molecolare abbiamo cercato di verificare la
capacità dei simulatori di riprodurre il ripiegamento di una sequenza di aminoacidi,
come normalmente avviene per le proteine; gli esperimenti condotti sono stati
essenzialmente due: verificare la stabilità strutturale di una proteina nota, e la ricerca
del minimo globale nella PES di una breve sequenza di aminoacidi.
STABILITÀ CONFORMAZIONALE
Dalla banca dati PDB abbiamo scelto il peptide sintetico 1DNG; la sua struttura non è
determinata da un ambiente complesso come è il citoplasma; inoltre è breve (15
residui) e con una sola α-elica, che costituisce il 53,33% della struttura:
QAPAYEEAAE ELAKS
S HHHHH HHHTT
La struttura del polipeptide è stata determinata sperimentalmente tramite NMR, e
quindi è un buon punto di riferimento per la verifica dell’evoluzione della
simulazione.
Per questo esperimento abbiamo usato una serie di cicli di meccanica e dinamica
molecolare, per vedere l’evoluzione nel tempo della struttura e a quale conformazione
la ricerca della minima energia potenziale avrebbe portato la molecola.
Per la minimizzazione della molecola abbiamo usato il campo di forze Amber94, che
è ottimizzato per le sequenze aminoacidiche.
La dinamica si svolge in un ambiente canonico, dove cioè la temperatura è costante;
tenendo conto che i tempi delle vibrazioni molecolari sono dell’ordine del
femtosecondo, e che i tempi di folding sono perlomeno dell’ordine dei 100
picosecondi, abbiamo stabilito in 120.000 il numero di iterazioni per ogni ciclo di
dinamica.
Gli altri parametri delle dinamica non sono stati cambiati; in particolare la
temperatura è stata mantenuta a 300K, per conservare le tipiche condizioni
dell’ambiente cellulare.
Le configurazioni campionate sono state 1000; tutte sono state minimizzate, e la
migliore è stata utilizzata per il successivo ciclo di dinamica.
I valori iniziali della molecola sono riportati in Tabella 1. Il calcolo dei parametri
energetici della struttura 1DNG prima della minimizzazione rivela un valore di
energia interna positivo, probabilmente imputabile al fatto che la struttura depositata
in 1DNG è la struttura media di diverse soluzioni NMR.
E
str
ang
stb
tor
oop
Iniziale
182,676 35,025 15,473 0,060 82,593 0,060
Minimizzata -430,644 9,930 28,407 2,239 70,568 2,239
ele
vdw
0,060 44,070
-216,438 -12,223
sol
E - sol
5,335 177,341
-315,364 -115,280
Tabella 1 - Livelli di energia del peptide 1DNG
In soli tre cicli abbiamo notato che l’energia potenziale tende a stabilizzarsi, e quindi
abbiamo interrotto l’esperimento, ritenendo che in queste condizioni non ci sarebbero
stati miglioramenti sostanziali. La serie di grafici che seguono (Figura 1, Figura 2 e
Figura 3) mostrano l’andamento dell’energia nei tre cicli che abbiamo eseguito: per
ogni campione viene mostrato in blu il livello di energia potenziale assunto durante la
dinamica, e in giallo l’energia della configurazione dopo la minimizzazione. Si può
vedere come l’energia della molecola si stabilizzi intorno a –250 Kcal/mol già al
primo ciclo, e resti stabile sugli stessi valori anche nelle successive iterazioni.
Figura 1 – Energie nel primo ciclo di dinamica e minimizzazione
Figura 2 - Energie nel secondo ciclo di dinamica e minimizzazione
Figura 3 - Energie nel terzo ciclo di dinamica e minimizzazione
La seconda serie di grafici mostra la relazione tra l’energia potenziale e l’energia
minimizzata della molecola. Nel primo ciclo, Figura 4, è interessante notare come la
maggior parte delle conformazioni campionate durante la dinamica siano compresse
in un intervallo piuttosto ridotto di energia dopo la minimizzazione. Questa
rappresentazione rafforza l’idea che ci sia un grande bacino di attrazione in cui
cadono moltissime conformazioni. Nel secondo e terzo ciclo il risultato della
minimizzazione sembra essere indipendente dalla scelta del campione, e tutto ricade
nello stesso minimo di energia.
Figura 4 - Relazione tra energia potenziale e minimizzata al primo ciclo
Figura 5 - Relazione tra energia potenziale e minimizzata al secondo ciclo
Figura 6 - Relazione tra energia potenziale e minimizzata al terzo ciclo
La stabilità riscontrata potrebbe dipendere dalla scelta della temperatura a cui è stata
condotta la simulazione: i 300K utilizzati potrebbero non essere sufficienti a superare
una ipotetica barriera energetica della PES, e limitare il campionamento in un’area
diversa da quella del minimo globale. Per questo, abbiamo effettuato un ulteriore
ciclo di meccanica e dinamica a 800K, ma i migliori risultati delle energie potenziali
minimizzate sono risultate paragonabili ai precedenti, e quindi abbiamo abbandonato
la ricerca. Ancora una volta i grafici mostrano che il punto stazionario trovato ha
sempre lo stesso livello di energia, nonostante i livelli di energia potenziale dei
campioni siano notevolmente più alti, come ci si aspettava data la notevole
temperatura del sistema.
Figura 7 - Energie nel ciclo di dinamica e minimizzazione a 800K
Figura 8 - Relazione tra energia potenziale e minimizzata al ciclo a 800K
Anche se è difficile trarre delle conclusioni definitive usando queste metodologie, i
nostri esperimenti ci portano a pensare di avere effettivamente trovato un bacino di
attrazione molto forte, che forse non coincide con il minimo energetico assoluto da un
punto di vista matematico, ma che certamente stabilizza notevolmente la proteina,
rendendola insensibile anche a variazioni “traumatiche” del sistema.
Di seguito mostriamo l’evoluzione della struttura durante la sequenza di meccanica e
dinamica molecolare. L’apparente perdita di struttura visibile nelle immagini non
contraddice i risultati numerici visti sopra, perché se una piccola variazione di
posizione può influire sulla rappresentazione grafica della struttura eliminando la
freccia che indica l’α-elica, essa non determina variazioni fondamentali nella
molecola reale.
Figura 9 – Struttura tridimensionale originale
Figura 10 – Struttura tridimensionale dopo il primo ciclo
Figura 11 – Struttura tridimensionale dopo il secondo ciclo
Figura 12 – Struttura tridimensionale dopo il terzo ciclo
Figura 13 – Struttura tridimensionale dopo il ciclo a 800K
RICERCA DEI PUNTI STAZIONARI
Dalla sequenza della proteina 1DNG abbiamo selezionato i cinque aminoacidi dell’αelica
Ala – Glu – Glu – Leu – Ala.
Partendo dalla struttura lineare abbiamo usato la ricerca sistematica e la ricerca
stocastica alla ricerca dei punti di minimo, possibilmente globale.
La ricerca sistematica
La ricerca sistematica è un algoritmo deterministico, la cui capacità di garantire il
risultato ottimale dipende solo dalla risoluzione con cui si analizza la PES. Poiché la
PES è una superficie continua, la ricerca può essere potenzialmente estremamente
dispendiosa dal punto di vista del calcolo, quindi diventa indispensabile trovare il
compromesso migliore tra i tempi di esecuzione dell’algoritmo e la necessità di avere
una buona risoluzione della PES: intervalli troppo ampi nella discretizzazione
potrebbero fare mancare il minimo globale.
Le conformazioni vengono generate dalle rotazioni degli angoli di torsione.
Utilizzando dei passi di 60 o 120 gradi per i 17 angoli diedri della sequenza in esame,
il numero di conformazioni generabili supera i 16 miliardi. Dal grafico di
Ramachandran relativo ai cinque aminoacidi scelti, considerati all’interno di 1DNG
(Figura 14), si può vedere come gli angoli φ e ψ non sono multipli di 60º, e quindi
non verranno mai individuati dall’algoritmo con il passo dato.
Figura 14 – Grafico di Ramachadran
Inoltre, dati i limitati mezzi di calcolo, si è reso necessario un ulteriore aumento dei
vincoli di ricerca: abbiamo usato la funzione di minimizzazione del passo offerta da
MOE che incrementa gli angoli in modo da campionare il minimo locale dell’energia
potenziale per quel legame. Con questa euristica il numero di conformazioni trovate si
è ridotto a poco più di 400.000.
Di tutte abbiamo calcolato l’energia, e minimizzato il 10% con la minore energia. In
Tabella 2 riportiamo le 10 migliori conformazioni.
Struttura
1
2
3
4
5
6
7
8
9
10
U (Kcal/mol)
-135,28
-134,70
-134,00
-133,91
-133,87
-130,40
-128,39
-128,29
-127,52
-126,67
Tabella 2 - Energia potenziale delle strutture minimizzate
In Figura 15 è riportata la struttura della conformazione a minore energia potenziale,
dalla quale risulta immediatamente che il ripiegamento trovato non ha, purtroppo,
alcuna somiglianza con l’elica formata dalla sequenza da cui deriva il nostro
pentapeptide. Probabilmente, quindi, l’euristica utilizzata non permette la
soddisfacente esplorazione della PES.
Figura 15 – Configurazione a minore energia potenziale risultante dalla ricerca sistematica
La ricerca stocastica
Il campionamento casuale della PES permette di ottenere una “visione d’insieme”
dello spazio delle conformazioni della molecola molto più velocemente rispetto alla
ricerca sistematica; d’altro canto, l’algoritmo non è in grado di garantire che troverà la
soluzione ottimale, se non in tempi molto lunghi (tendenti all’infinito). A questo si
aggiunge, come nel caso della dinamica molecolare, il problema della terminazione,
perché non ci sono strumenti matematici in grado di valutare la probabilità di avere
raggiunto effettivamente il minimo globale.
Nella nostra ricerca abbiamo mantenuto i parametri proposti da MOE, escludendo
però l’inversione dei centri chirali nella generazione delle strutture, dal momento che
in natura la conformazione del carbonio α è sempre levogira. Le strutture analizzate
durante la ricerca sono state 10.000, da cui sono state selezionate sei configurazioni,
che abbiamo minimizzato, e i cui livelli di energia sono riportati in Tabella 3.
Struttura
1
2
3
4
5
6
U (Kcal/mol)
-72,37
-69,22
-68,90
-68,19
-67,93
-67,28
Tabella 3 – Energia potenziale delle strutture minimizzate
Della configurazione a minore energia potenziale è riportata la struttura in Figura 16,
dalla quale risulta che il ripiegamento trovato non ha neppure in questo caso alcuna
somiglianza con un α-elica. Il numero di strutture esplorate è decisamente troppo
piccolo per poter pensare di ottenere risultati interessanti, e confrontando questo
numero con quelli generati dalla ricerca sistematica l’esplorazione della PES è
decisamente superficiale. Ciononostante l’ordine di grandezza dell’energia delle
conformazioni trovate è paragonabile a quello delle conformazioni trovate dalla
ricerca sistematica. In più è interessante notare come le conformazioni generate con i
due differenti algoritmi siano simili, almeno nella disposizione dello scheletro, che si
ripiega a ferro di cavallo. Le differenti disposizioni delle catene laterali sono
evidentemente le principali responsabili delle differenze energetiche tra le due
molecole.
Figura 16 – Configurazione a minore energia potenziale risultante dalla ricerca stocastica
Dare una valutazione della qualità dei risultati non è semplice: in generale non si
hanno informazioni sulla struttura che si dovrebbe ottenere, e quindi l’esperienza del
ricercatore diventa un fattore fondamentale per capire se la ricerca sta dando dei frutti,
se si continuando l’esplorazione si possono ottenere miglioramenti, e, in definitiva,
quando interrompere la ricerca. Nel nostro caso constatiamo che due metodologie
diverse convergono a risultati simili, e questo ci porta a pensare che con questo livello
di dettaglio nella ricerca non si possono ottenere risultati migliori di quelli ottenuti.