Tesi 2
Tesi 2
Tesi 2
Introduzione I
1 Le proteine 2
1.1 La struttura delle proteine . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Folding vs Misfolding . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Ruolo biologico dei metalli . . . . . . . . . . . . . . . . . . . . . 15
2 La malattia di Alzheimer 20
2.1 Fisiologia e patologia . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2 La diagnosi della malattia . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Il peptide β -amiloide . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4 Ioni metallici e Aβ peptide . . . . . . . . . . . . . . . . . . . . . 36
3 Dinamica molecolare 40
3.1 Teoria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Force field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.3 Gli ensemble statistici . . . . . . . . . . . . . . . . . . . . . . . . 52
3.4 Schematizzazione dell’algoritmo di dinamica
molecolare classica . . . . . . . . . . . . . . . . . . . . . . . . . 60
i
INDICE
ii
INDICE
G Pseudopotenziali 182
Bibliografia 200
iii
Introduzione
Negli ultimi anni, è stato dimostrato che numerose malattie dipendono dall’al-
terazione del folding proteico. Tali malattie, raggruppate sotto il nome di Protein
Conformational Disorders (PCD), comprendono la malattia di Alzheimer (AD),
le encefalopatie spongiformi, la corea di Huntington, il morbo di Parkinson, il dia-
bete di tipo II, l’amiloidosi da dialisi, la sclerosi laterale amiotrofica (SLA) e più
di altre quindici malattie meno note. L’evento chiave per lo scatenarsi delle PCD
è un cambiamento nella struttura secondaria e/o terziaria di una proteina o di un
suo frammento, non accompagnato da alterazioni della struttura primaria. In tutti
i casi, c’è una progressiva transizione dalla proteina correttamente avvolta verso
uno stato aggregato, ricco di conformazioni beta.
I
INDICE
la principale e irrisolta questione. Tre differenti ipotesi sono state proposte per
descrivere le relazioni tra stati conformazionali e aggregazione. L’ipotesi del-
la polimerizzazione (A) afferma che è l’aggregazione ad indurre cambiamen-
ti nella conformazione proteica; quella conformazionale (B) che il misfolding
proteico è indipendente dall’aggregazione, che non rappresenta un necessario
punto d’arrivo del cambiamento conformazionale; infine l’ipotesi conformazio-
ne / oligomerizzazione (C), rappresenta una visione intermedia, nella quale pic-
coli cambiamenti conformazionali innescano l’oligomerizzazione che è essenziale
per la stabilizzazione del misfolding proteico.
II
INDICE
Nell’Aβ si possono osservare due regioni: (i) la regione compresa tra i residui
29 e 40, corrispondente alla regione transmembrana dell’APP, altamente idrofo-
bica, scarsamente solubile, e tendente ad aggregare in strutture dissimili da quelle
viste nei pazienti di AD e differenti a seconda del pH, del solvente e della tempe-
ratura; (ii) il frammento 1-28, che occupa la regione extracellulare dell’APP e può
produrre sia strutture ad α-elica che oligomeri β -sheet simili a quelli caratteristici
delle placche amiloidee e quindi responsabile della variabilità conformazionale
dell’Aβ . E’ stato osservato che le placche contengono quantità superiori alla
media di ioni metallici appartenenti al gruppo d, come lo zinco (Zn2+ ), il rame
(Cu2+ ) e il ferro (Fe3+ ).
Si è osservato che lo Zn2+ più degli altri ioni metallici favorisce l’aggregazione
di Aβ . Tuttavia il ruolo dello Zn2+ è soggetto a continui dibattiti: (i) alcune
osservazioni sperimentali sembrano indicare un possibile ruolo dello Zn2+ nel-
la formazione degli aggregati iniziali, in cui lo ione metallico fa da “ponte” fra
due o più peptidi differenti formando dei legami cosiddetti “inter-molecolari”; (ii)
mentre altri ipotizzano un ruolo centrale dello Zn2+ nel processo del misfolding.
Per quanto riguarda la coordinazione del metallo, alcuni studi hanno mostrato
che la regione che lega il metallo è compresa tra gli amminoacidi 6 e 28 del pep-
tide Aβ ed è stato suggerito il coinvolgimento dei due amminoacidi istidina-13 e
istidina-14. In recenti studi NMR sono stati proposti vari modi di coordinazione
intra-peptide in cui lo ione metallico si lega con tre istidine (6-13-14) e con la
regione ammino-terminale oppure con l’acido glutammico-11. Inoltre un recente
lavoro ha suggerito che il frammento peptidico minimo coinvolto nell’interazione
con il metallo sia Aβ1−16 .
Per tentare di contribuire a chiarire il possibile ruolo svolto dallo Zn2+ è sta-
III
INDICE
IV
INDICE
V
Capitolo 1
Le proteine
2
1.1 La struttura delle proteine
Nelle proteine sono presenti 20 tipi di catene laterali che variano per dimen-
sioni, carica, capacità di formare legami idrogeno e reattività chimica. Tutte le
proteine in tutte le specie - batteri, archeobatteri ed eucarioti - vengono costruite a
partire dagli stessi 20 amminoacidi [1]. Le catene laterali degli amminoacidi han-
no differenti tendenze a partecipare alle interazioni l’una con l’altra e con l’acqua.
Queste differenze influenzano profondamente i loro contributi alla stabilità e alle
funzioni delle proteine.
3
Le proteine
Figura 1.2: Struttura e caratteristiche chimiche delle catene laterali degli amminoacidi
4
1.1 La struttura delle proteine
Gli amminoacidi anfipatici (fondo viola in Fig.1.2) possono essere sia idrofo-
bici sia idrofilici a seconda delle condizioni del microambiente in cui si trovano.
Questa loro caratteristica è ideale per la formazione delle interfacce [2].
5
Le proteine
Figura 1.4: Catena polipeptidica estesa in cui sono mostrate le tipiche lunghezze di legame, gli
angoli di legame e i diedri
6
1.1 La struttura delle proteine
che molte combinazioni sono proibite a causa di impedimenti sterici tra gli atomi.
I valori permessi possono essere visualizzati attraverso un grafico bidimensionale
detto appunto grafico di Ramachandran. Nella Fig.1.5 è mostrato un grafico di
Ramachandran in cui con il colore rosso indichiamo quelle combinazioni degli
angoli torsionali che sono “permesse”, nelle quali cioè non ci sono impedimenti
sterici, con il rosa le regioni permesse solo se alcuni impedimenti sterici sono ri-
lassati. Come si vede dal grafico, circa tre quarti delle possibili combinazioni tra
gli angoli torsionali vengono escluse da impedimenti sterici.
7
Le proteine
1.5Å lungo l’asse dell’elica e forma con esso un angolo di 100◦ il che significa
che vi sono 3.6 residui amminoacidici per ogni giro di elica. Il passo dell’α-elica,
8
1.1 La struttura delle proteine
che corrisponde al prodotto dello spostamento (1.5Å) per il numero di residui per
giro (3.6Å), è di 5.4Å. Le α-eliche sono una struttura compatta [6], con valori
approssimativi degli angoli φ e ψ di -60◦ e -50◦ rispettivamente (vedi Fig.1.5).
Per buona parte delle proteine, la struttura terziaria rappresenta l’ultimo livello
di organizzazione strutturale. E’ il caso delle proteine cosiddette monomeriche,
costituite cioè da un’unica unità funzionale, biologicamente attiva. Molte altre
proteine (ad esempio, un gran numero di enzimi), nella loro forma attiva sono
invece costituite dall’associazione di due o più unità di struttura terziaria (dette
monomeri o subunità), uguali (proteine omo-oligomeriche) o diverse (proteine
etero-oligomeriche). Si parla in tal caso di struttura quaternaria, per riferirsi all’or-
ganizzazione multimerica della proteina. Nella struttura quaternaria, le subunità
sono tenute insieme da interazioni generalmente non covalenti. Raramente, più
catene peptidiche sono unite da legami covalenti, come accade ad esempio nelle
immunoglobuline, in cui le catene leggere e pesanti sono tenute insieme da ponti
disolfuro [8].
9
Le proteine
10
1.2 Folding vs Misfolding
11
Le proteine
12
1.2 Folding vs Misfolding
13
Le proteine
14
1.3 Ruolo biologico dei metalli
15
Le proteine
16
1.3 Ruolo biologico dei metalli
17
Le proteine
18
Capitolo 2
La malattia di Alzheimer
20
2.1 Fisiologia e patologia
• Demenza vascolare.
• Demenze da prioni.
Si stima che circa il 50-60% dei casi di demenza sia dovuto alla malattia di
Alzheimer (Alzheimer’s disease - AD), il 20% sia associato a demenze vascolari e
a corpi di Lewy, il 15% a malattie neurodegenerative quali la corea di Huntington
e alla malattia di Parkinson ed il restante sia dovuto ad altre cause.
L’AD è un disordine neurodegenerativo eterogeneo con progressione irreversi-
bile che interessa primariamente le regioni ippocampali e neocorticali del cervello
[22]. Questa malattia fu descritta per la prima volta nel 1907 dal neuropsichiatra
tedesco Alois Alzheimer. L’AD è la quarta causa di morte in età adulta - dopo le
malattie cardiovascolari, il cancro e l’ictus - e colpisce circa 20 milioni di persone
nel mondo (oltre 450,000 in Italia). Per ragioni facilmente intuibili, legate alla
crescente aspettativa di vita, le previsioni per il futuro sono catastrofiche. Questa
patologia è destinata ad aumentare e il numero di pazienti tra vent’anni potrebbe
21
La malattia di Alzheimer
più che raddoppiare. I costi per lo Stato sono già enormi e la scoperta di farmaci in
grado di ritardare anche soltanto di cinque anni l’esordio della malattia potrebbe
far risparmiare al servizio sanitario nazionale un numero cospicuo di milioni di
euro e/o dollari.
Al momento non esiste alcuna cura efficace contro di essa e la patogenesi della
malattia resta ancora oggetto di molte teorie che coinvolgono fattori genetici ed
epigenetici (come fattori metabolici od ambientali).
22
2.2 La diagnosi della malattia
La diagnosi avviene solo con la congiuntura dei suddetti esiti, non è mai sicura e
certa, ma solo possibile e probabile. La diagnosi sicura avviene solo dall’autopsia
e dalle analisi istologiche del tessuto cerebrale.
Dall’esame autoptico emerge una riduzione nelle dimensioni delle regioni del-
l’encefalo coinvolte nei processi implicati nella memorizzazione, nell’apprendi-
mento e nei comportamenti emotivi. Queste zone comprendono la corteccia en-
torinale, l’ippocampo, il prosencefalo e l’amigdala (Fig.2.1).
La patologia
23
La malattia di Alzheimer
24
2.2 La diagnosi della malattia
Il fattore genetico
25
La malattia di Alzheimer
26
2.2 La diagnosi della malattia
La terapia
Figura 2.5: Sinapsi colinergica con i principali siti d’azione delle molecole attive sulla
trasmissione mediata dall’acetilcolina
27
La malattia di Alzheimer
dei pazienti. Tre studi longitudinali1 hanno evidenziato un effetto protettivo della
terapia di sostituzione con estrogeni sul rischio di Alzheimer [30]. Gli estrogeni
modulano il metabolismo secretorio dell’APP e un polimorfismo del gene alfa
del recettore degli estrogeni sembra interagire con l’Apo ε-4 nel determinismo di
forme sporadiche di Alzheimer.
Inoltre ai malati di AD vengono somministrati anche farmaci antidepressivi,
antipsicotici, ansiolitici e sedativi, soprattutto nelle fasi terminali della malattia.
Negli ultimi anni numerosi studi sono stati indirizzati alla ricerca di molecole
naturali o di nuova sintesi in grado di inibire l’aggregazione dell’Aβ e ridurne
gli effetti citotossici. Ricerche di questo tipo potrebbero rappresentare il punto di
partenza per la progettazione di terapie farmacologiche di prevenzione e/o cura
e, inoltre, contribuire alla comprensione dei dettagli molecolari che regolano il
processo di fibrillogenesi.
Molti studi dimostrano che piccole molecole con anelli aromatici sono in gra-
do di influenzare significativamente la cinetica di aggregazione dei peptidi Aβ .
Queste includono la curcumina, l’acido rosmarinico, il fullerene, le tetracicline,
la melatonina, i polifenoli del vino, ecc.
In prospettiva, le strategie più promettenti per il trattamento dell’AD sembrano
essere tre:
28
2.2 La diagnosi della malattia
Per ultimo vogliamo citare il Progetto NAD (Nanoparticles for therapy and
diagnosis of Alzheimer Disease), partito nel settembre 2008. Un progetto di ricer-
ca multidisciplinare con l’obiettivo di diagnosticare precocemente e contrastare
la malattia di Alzheimer. L’obiettivo è quello di realizzare nanoparticelle (NPs)
in grado di attraversare la barriera emato-encefalica per raggiungere il cervello,
sede principale della malattia di Alzheimer. Alle nanoparticelle verranno legate
molecole in grado di riconoscere e distruggere le placche amiloidi che si deposi-
tano nel cervello in tale malattia. Il progetto NAD coinvolge diciannove partner
tra centri di ricerca e piccole e medie imprese provenienti da Italia, Spagna, Por-
togallo, Francia, Slovacchia, Svezia, Olanda, Ungheria, Finlandia, Grecia, Belgio
e Inghilterra, impegnati in un connubio multidisciplinare
29
La malattia di Alzheimer
30
2.3 Il peptide β -amiloide
31
La malattia di Alzheimer
Figura 2.6: Schema della struttura di APP e della sua maturazione proteolitica
32
2.3 Il peptide β -amiloide
Notiamo che il peptide contiene 3 istidine (residui 6, 13 e 14), che sono am-
minoacidi spesso coinvolti nel legame con i metalli, come evidenziato nel §1.3.
Il metallo si lega generalmente all’istidina attraverso uno dei due atomi di azoto
che fanno parte dell’anello imidazolico. Studiando il profilo di idropaticità del
33
La malattia di Alzheimer
34
2.3 Il peptide β -amiloide
35
La malattia di Alzheimer
Figura 2.11: Schema del processo che porta alla formazione delle fibrille
che metalli come zinco, rame e ferro sono presenti nei depositi amiloidi della
malattia di Alzheimer, e che essi sono in grado di indurre l’aggregazione in vitro
del peptide Aβ e la conseguente formazione di fibrille.
Vi sono almeno due generiche reazioni con i metalli che devono essere prese
in considerazione. La prima è relativa all’associazione di un metallo ad una pro-
teina che porta a misfolding e/o ad aggregazione proteica. Questa reazione può
coinvolgere ioni metallici non redox come lo Zn2+ , o metallo redox attivi come il
Cu2+ e il Fe3+ [40]. La seconda è relativa all’ossidazione proteica che, catalizza-
ta dai metalli, può dar luogo ad un danno proteico o alla denaturazione. Questa
reazione coinvolge solo metalli redox.
36
2.4 Ioni metallici e Aβ peptide
Sono stati riportati risultati contradditori su: i) speciazione, ii) costanti di sta-
bilità, iii) siti di legame e ambienti di coordinazione, iv) strutture conformazionali
di metalli del gruppo d con le proteine o con loro frammenti peptidici coinvolti
nell’AD.
Si suppone che l’aggregazione del peptide Aβ sia mediata dall’interazione
con i metalli, e tale ipotesi è avvalorata dal fatto che il peptide ha dei siti con alta
affinità per gli ioni Zn2+ e Cu2+ , ed in misura minore per Fe3+ .
E’ stato mostrato che lo Zn2+ ha una capacità di promuovere l’aggregazione
maggiore rispetto a tutti gli altri metalli [41]. Sono stati condotti numerosi stu-
di usando varie tecniche sperimentali. Alcuni di essi hanno mostrato un possi-
bile coinvolgimento dello ione Zn2+ nel processo di polimerizzazione. I model-
li proposti hanno mostrato che lo Zn2+ è alternativamente coinvolto in legami
intermolecolari [42] e intramolecolari [43].
Studi recenti [43] hanno identificato il frammento minimo coinvolto nel lega-
me col metallo, che risulta essere il frammento Aβ1−16 , la cui sequenza è:
37
La malattia di Alzheimer
Figura 2.13: Andamento temporale dello stato cognitivo di 33 pazienti con diagnosi di probabile
AD
38
Capitolo 3
Dinamica molecolare
Anche il più complesso dei sistemi biologici non è altro che un sistema fisi-
co, composto da atomi, e deve pertanto poter essere analizzato con gli strumen-
ti della fisica. In particolare, il sistema deve essere descritto da una Hamilto-
niana, che a rigore dovrebbe essere trattata in modo quantistico che permetta di
determinarne tutte le proprietà chimico-fisiche. Tuttavia la soluzione esatta del-
l’equazione di Schröedinger è impensabile già per sistemi molto semplici come
le piccole molecole e in genere è sempre necessario adottare approssimazioni
empiriche, giustificate a posteriori solo dalla bontà dei risultati raggiunti. Tutto
questo è a maggior ragione vero per un sistema biologico come una proteina, che
è un insieme di decine di migliaia di atomi. Quindi nel caso in cui si abbia a che
fare con sistemi composti da un numero elevato di atomi, l’approccio quantistico
non è più computazionalmente praticabile, e si ricorre a tecniche alternative, come
la dinamica molecolare classica. Questa tecnica permette di esplorare lo spazio
delle fasi di un sistema nell’ensemble microcanonico.
40
3.1 Teoria
3.1 Teoria
Lo stato all’istante tn di un sistema composto da N atomi può essere rappre-
sentata da un vettore dello spazio delle fasi a 6N dimensioni:
Γ(~r1 (tn ),~r2 (tn ), . . . , r~N (tn ); ~p1 (tn ), ~p2 (tn ), . . . , p~N (tn )) (3.1)
e la sua evoluzione sarà governata dalle equazioni di Hamilton per i vettori po-
sizione e momento:
41
Dinamica molecolare
E’ chiaro che non è possibile procedere in questo modo, poiché il numero di ope-
razioni richiesto è eccessivo per qualunque calcolatore. Esistono vari approcci
alla risoluzione di tale problema, il più raffinato è l’algoritmo di Multiple Time
Step [46], che permette di disaccoppiare i moti veloci da quelli lenti, consentendo
di integrare, con passi differenti, gradi di libertà che corrispondono a frequenze
differenti, permettendoci cosı̀ di velocizzare i tempi di calcolo integrando in modo
appropriato le varie componenti del potenziale. Per prima cosa osserviamo che,
formalmente, la soluzione delle equazioni di Hamilton può scriversi nella forma:
e dove:
N
d ∂ ∂V ∂
iL = ∑ ~ri − (3.7)
i=1 dt ∂~ri ∂~ri ∂ ~pi
42
3.1 Teoria
dove abbiamo posto t = p0 ∆t0 . Dunque p0 è il numero totale dei passi che si
devono effettuare per integrare le equazioni del moto, se ciascun passo evolve
la dinamica per il tempo caratteristico dei moti lenti. A questo punto conviene
suddividere l’operatore differenziale 3.7 in una parte lenta ed una veloce:
dove abbiamo posto ∆t0 = p1 ∆t1 . La formula 3.11 consente di effettuare l’ag-
giornamento dei gradi di libertà che corrispondono ai moti veloci tenendo fissi
quelli che invece corrispondono ai moti lenti. Inoltre si vede che per ogni passo
che viene effettuato sui moti lenti, se ne effettuano p1 su quelli veloci. In questo
modo è possibile realizzare concretamente e rigorosamente la separazione tra i
moti lenti e quelli veloci ed integrare meno frequentemente quelli che richedono
un maggior costo computazionale.
43
Dinamica molecolare
dove:
∂ ~p ∂
iL1 = ~F(~q) , iL2 = (3.14)
∂~p m ∂~q
ed in cui abbiamo definito:
~F(~q) = −∇V (~q) (3.15)
44
3.2 Force field
Nella 3.2 abbiamo assunto che il potenziale non dipenda dagli impulsi, cioè
che non ci sia attrito tra i vari componenti del sistema. L’approccio della dinamica
molecolare consiste nello studiare sistemi di atomi immaginati come punti mate-
riali soggetti ad un potenziale efficace di interazione da determinare fenomeno-
logicamente. A partire da questo si costruisce il campo di forze che determina
il moto degli atomi (force field (ff)). Per rendere questo processo più semplice
si considerano come gradi di libertà quei parametri fisici che sappiamo essere
particolarmente importanti nei legami molecolari, e la cui variazione influisce in
maniera consistente sull’energia del sistema: come per esempio le lunghezze di
legame, gli angoli di legame e gli angoli di torsione. Il potenziale effettivo espres-
so in queste variabili ha una forma semplice e i vari termini che lo compongono
possono essere facilmente confrontati con gli esperimenti. Limitarsi ai gradi di
libertà detti e considerare la struttura interna delle molecole rigida (legami), per-
mette di descrivere comunque una larga parte dello spazio delle configurazioni, e
poiché i termini coinvolgono solo pochi atomi, e non strutture più complesse, pos-
sono essere utilizzati per tutti i tipi di molecole che coinvolgano lo stesso tipo di
45
Dinamica molecolare
46
3.2 Force field
47
Dinamica molecolare
si può pensare che anche l’angolo di legame resterà vicino al suo valore di
equilibrio, ed è pertanto possibile sviluppare anche questo termine in serie.
Fermandosi nuovamente all’ordine più basso si ha un contributo:
1 2
Vθ (θi jk ) = kθ ,i jk θi jk − θ0,i jk (3.25)
2
• Infine si tiene conto anche dell’angolo torsionale (Fig.3.3) ovvero della tor-
48
3.2 Force field
49
Dinamica molecolare
50
3.2 Force field
51
Dinamica molecolare
52
3.3 Gli ensemble statistici
(3.32) per t → ∞. Perchè ciò avvenga è necessario che siano verificate due ipotesi
fondamentali:
2. La dinamica del sistema è ergodica, ovvero c’è una probabilità non nulla di
raggiungere qualunque configurazione accessibile dello spazio delle confi-
gurazioni in un tempo finito. Questa condizione va verificata caso per caso
ed è in genere molto difficile da dimostrare.
Una volta soddisfatte queste due condizioni possiamo simulare un sistema tramite
una opportuna dinamica e misurare le proprietà medie lungo la simulazione. Ve-
diamo ora quali dinamiche sono state sviluppate per riprodurre ensemble di parti-
colare interesse.
L’ensemble microcanonico
53
Dinamica molecolare
L’ensemble canonico
I termostati
54
3.3 Gli ensemble statistici
dove:
1
δt T0 2
λ = 1+ −1 (3.39)
τT T
La costante τT è la costante di accoppiamento temporale che determina la scala
temporale entro la quale la temperatura desiderata è raggiunta. E’ semplice mo-
strare che il termostato proporzionale ha una distribuzione maxwelliana.
Un’altra categoria di termostati sono quelli stocastici, in cui tutti o un sottoin-
sieme di gradi di libertà del sistema sono soggetti a collisioni con particelle vir-
tuali. Questo metodo può essere simulato per mezzo di un equazione differenziale
stocastica di Langevin che descrive il moto di una particella dovuto all’agitazione
55
Dinamica molecolare
∂~pi ∂U
=− − γ~pi + ~F + (3.40)
∂t ∂~qi
56
3.3 Gli ensemble statistici
è:
N ~πi2 πs2
H∗ = ∑ 2
+U(~
ρ ) + + gkB T ln s (3.43)
i=1 2mi s 2Ms
∂ ~ρi ~πi
= 2 (3.44)
∂τ s
∂ ~πi ∂U(~ρ )
= (3.45)
∂τ ∂ ~ρi
∂s πs
= (3.46)
∂τ Ms
∂ πs 1 N ~πi2 gkB T
= 3∑ − (3.47)
∂τ s i=1 mi s
∂~qi ~pi
= (3.48)
∂t mi
∂~pi ∂U(~q)
=− − ζ~pi (3.49)
∂t ∂~qi
∂ ln s
=ζ (3.50)
∂t
!
N
∂ζ 1 ~p2i
= ∑ − gkB T (3.51)
∂t Ms i=1 mi
57
Dinamica molecolare
I barostati
V̇ = 3αV (3.53)
58
3.3 Gli ensemble statistici
Poiché L è una quantità dinamica, il momento non è legato in modo semplice alla
derivata temporale delle coordinate ma
∂~qi ∂ ~ρi ∂L
=L + ~ρi (3.59)
∂t ∂t ∂t
L’Hamiltoniana estesa del sistema sarà allora data da:
N ~πi ~πV2
1
H∗ = ∑ +U(V 1/3~
ρ ) + PexV + (3.60)
V 2/3 i=1 2mi 2MV
dove Pex è la pressione esterna imposta e πV ed MV sono un momento ed una
massa associate alle fluttuazioni del volume. Le equazioni del moto derivate dalla
Hamiltoniana (3.60) sono
∂ ~ρi 1 ~πi
= 2/3 (3.61)
∂t V mi
59
Dinamica molecolare
∂~qi ~pi 1 ∂V
= + ~qi (3.65)
∂t mi 3V ∂t
∂~pi ∂U(~q) 1 ∂V
= − ~pi (3.66)
∂t ∂~qi 3V ∂t
∂V pV
= (3.67)
∂t MV
!
N
∂ pV 1 1 ~p2i ∂U(~q)
= 2/3 ∑ −~qi − Pex (3.68)
∂t 3v V i=1 mi ∂~qi
Queste ultime equazioni sono state ricavate utilizzando il metodo di Parrinello-
Rahman [60], che assicura che il sistema esplori un ensemble canonico.
60
3.4 Schematizzazione dell’algoritmo di dinamica
molecolare classica
noti temperatura e pressione iniziale vengono calcolati i fattori di riscalamento
e/o di attrito associati al termostato e al barostato. A questo punto l’algoritmo
calcola le forze su ogni atomo del sistema risolvendo le equazioni del moto ed
aggiorna prima le velocità e poi le coordinate del sistema. Se abbiamo applicato
vincoli geometrici, vengono a questo punto calcolati, e prima le coordinate e poi
le velocità vengono corrette e salvate. L’algoritmo termina dopo che ha compiuto
un certo numero di passi.
61
Capitolo 4
63
Studio quantistico dei sistemi a molti corpi
dove il termine V (~r, ~R) rappresenta tutte le possibili interazioni tra le particelle del
sistema. Le variabili che compaiono nell’equazione (4.2) sono variabili multidi-
mensionali: ~R è una notazione abbreviata per le 3M coordinate nucleari e simil-
mente ~r è una notazione abbreviata per tutte le coordinate spaziali e di spin degli
elettroni. L’equazione di Schrödinger del sistema è:
h i
TN (~R) + Te (~r)V (~r, ~R) Ψ(~r, ~R) = W Ψ(~r, ~R) (4.3)
Gli autovalori W e le autofunzioni Ψ(~r, ~R) della (4.3) del sistema combinato
elettroni-nuclei sono detti rispettivamente “energie vibroniche” e “funzioni d’on-
da vibroniche” [61]. Classicamente possiamo visualizzare, nella nostra immagi-
nazione, gli elettroni come minuscole particelle leggere che si muovono rapida-
mente attorno ai pesanti nuclei. Questa visione pittoresca della fisica atomica ci
permette di poter fare intuitivamente un’assunzione molto importante. Sapendo
che la massa nucleare è molto più grande di quella elettronica possiamo trattar-
la come infinita. Cosı̀ facendo, l’operatore cinetico nucleare, TN nella (4.2), può
64
4.1 Approssimazione di Born-Oppenheimer
qui sia le funzioni φn (~r; ~R) che gli autovalori En (~R) dipendono parametricamente
dalle posizioni nucleari. Al variare di ~R gli autovalori En (~R) definiscono la cosı̀
chiamata superficie potenziale adiabatica ed i set φn (~r; ~R) le funzioni d’onda adi-
abatiche elettroniche. Con la presente approssimazione, che va sotto il nome di
approssimazione di Born-Oppenheimer [61], siamo riusciti a separare il problema
elettronico da quello nucleare, cioè, una volta “congelati” i gradi di libertà nucle-
ari possiamo risolvere il problema elettronico come se i nuclei fossero fermi. In
un secondo momento si riconsidera il problema nucleare e si vede l’evoluzione
del sistema in base alle nuove energie elettroniche. Il moto nucleare, a seconda
dei casi, può essere trattato sia classicamente che quantisticamente. Supponendo
una superficie adiabatica En (~R) non degenere (ben separata in energia da tutte le
altre), il sistema esplora il fondamentale foglio adiabatico E0 (~R) e l’equazione
classica del moto nucleare può essere scritta come:
d 2~RI ∂VN ∂ E0
MI 2
=− − (4.5)
dt ∂ ~RI ∂ ~RI
65
Studio quantistico dei sistemi a molti corpi
Nella (4.5) la forza agente su ciascun nucleo risulta costituita da due termini, uno
proveniente dal potenziale inter-ionico, l’altro legato alla variazione del contri-
buto adiabatico elettronico con le posizioni nucleari, cioè se si vuole, la forza di
“trascinamento” attraverso la quale agisce la “colla elettronica”. La derivazione
rigorosa della (4.5) nel limite classico costituisce il contenuto del teorema della
forza o di Hellmann-Feynmann [62]. Nel caso in cui si è obbligati a fare una
trattazione quantistica del moto nucleare otteniamo:
" #
2
h̄ ∂2
− + En (~R) χ(~R) = W χ(~R) (4.6)
2M ∂ ~R2
66
4.2 Problema elettronico
HΦ = EΦ (4.11)
dove
E = εi + ε j + . . . + εk (4.12)
67
Studio quantistico dei sistemi a molti corpi
68
4.3 L’approssimazione di Hartree-Fock
elettroni corrisponde a scambiare due righe del determinante di Slater, ossia, cam-
biare il segno del determinante. Ovviamente, un determinante di Slater è una
soluzione esatta dell’equazione di Schrödinger per un sistema di elettroni non in-
teragenti. Ciò nonostante, data l’ortonormalità degli spin orbitali, i determinanti
di Slater formano una base completa anche per lo spazio di Hilbert per un sistema
di elettroni interagenti, per cui l’esatta funzione d’onda può essere descritta da
una loro combinazione lineare. Il formalismo fin qui descritto è relativo ad una
formulazione di tipo one-electron per il sistema elettronico.
Adesso descriveremo quali sono le approssimazioni e i metodi usati per la
determinazione dell’operatore h(i), con il quale sia possibile scrivere l’equazione
agli autovalori (4.9).
N
1 N e2
H = ∑ hi (~ri ) + (4.16)
i=1 2 i6∑
= j ri j
69
Studio quantistico dei sistemi a molti corpi
hχi |χ j i = δi j (4.18)
1 N
A= ∑ (−1) pn Pn
N! n=1
(4.21)
Scegliamo quest’operatore in modo tale che esso sia hermitiano per cui vale A = A2
(operatore proiettivo), ed inoltre che commuti con l’Hamiltoniana.
Calcolando adesso E[Φ] utilizzando le (4.17), (4.19), (4.20), (4.21) e sfruttan-
do le proprietà dell’operatore A, otteniamo:
1
E[Φ] = ∑ Ii +
2∑ ∑[ℑi j − Ki j ] (4.22)
i i j
70
4.3 L’approssimazione di Hartree-Fock
δ E − ∑ ∑ εi j δ hχi |χ j i = 0 (4.23)
i j
χi = ∑ Ui j χ j (4.24)
j
71
Studio quantistico dei sistemi a molti corpi
e2
Z
V jd (q) = χ ∗j (q0 ) χ j (q0 )dq0 (4.28)
|r − r0 |2
Definiamo inoltre l’operatore di scambio in modo tale che quando agisce sugli
orbitali elettronici, si ottenga:
"Z #
e2
V jex (q)χi (q) = χ ∗j (q0 ) χ j (q0 )dq0 χ j (q) (4.29)
|r − r0 |2
Ze2
V (q) = − + ∑ V jd (q) + ∑ V jex (q) (4.30)
r j j
72
4.3 L’approssimazione di Hartree-Fock
Le (4.31) non sono semplici equazioni agli autovalori, una per ogni orbitale elet-
tronico, poiché lo stesso potenziale di HF dipende dagli orbitali elettronici.
Per risolvere il sistema di equazioni integro-differenziali di HF, si procede per
iterazione. Si parte da orbitali elettronici approssimati (scelta del “basis set” -
vedi Appendice D) e si calcola il potenziale di HF corrispondente. Si risolvono
numericamente le equazioni di HF ottenendo nuovi orbitali che a loro volta deter-
minano un nuovo potenziale di HF. La procedura viene ripetuta finché gli orbitali
elettronici producono un potenziale di HF che è identico (nell’approssimazione
desiderata) al potenziale ottenuto nel passo precedente. Il potenziale di HF ot-
tenuto in questo modo è noto come il campo self-consistent dell’atomo (dello
ione) [50].
Osservando la (4.31) verrebbe spontaneo interpretare Ei come l’autovalore
dell’energia del singolo elettrone, però prendendo il prodotto scalare della (4.28)
con χi , usando le definizioni di Ii , del termine di scambio e del termine diretto,
sommando su i e considerando la (4.22) otteniamo:
1 e2
E[Φ] = ∑ Ei − hΦ| |Φi (4.32)
i 2 i6∑
= j ri j
Si vede dunque che l’energia totale dell’atomo non è la somma delle singole ener-
gie. Ciò è dovuto al fatto che, sommando le energie dei singoli elettroni, ciascuna
energia cinetica e ciascuna energia di interazione con il nucleo viene contata una
volta, mentre l’energia d’interazione mutua, che ha come valor medio il secondo
membro della parte destra della (4.32), è contata due volte. Supponiamo adesso
di rimuovere l’i-esimo elettrone dal sistema, ed assumiamo che gli orbitali del si-
stema con N − 1 elettroni rimangano invariati dopo la rimozione, vediamo che la
differenza tra l’energia totale dei due sistemi è:
EN − EN−1 = Ei (4.33)
73
Studio quantistico dei sistemi a molti corpi
74
4.4 Teoria del funzionale densità
della DFT arriva nel 1964 con i teoremi di Hohenberg e Kohn [64] i quali posero
le basi per una rigorosa formulazione della DFT come oggi è conosciuta.
Nel nostro caso l’Hamiltoniana (4.34) coincide con la (4.4), vista nel §4.1 ed il
termine di potenziale esterno non è altro che il potenziale di interazione elettroni-
nuclei, questi ultimi considerati “congelati” in una configurazione assegnata. Il
teorema di Hohenberg e Kohn dà accesso ad alcuni risultati esatti che riguardano
lo stato fondamentale di un sistema di molti elettroni interagenti.
Teorema.
75
Studio quantistico dei sistemi a molti corpi
Dimostrazione.
C :W →Ψ (4.38)
che, per costruzione, è suriettiva1 , in quanto non esiste alcun elemento di Ψ che
non sia associato ad un elemento di W .
Se adesso consideriamo le densità elettroniche corrispondenti allo stato fon-
damentale:
n(~r) = hφs f |φs f i = φs∗f (~r)φs f (~r) (4.39)
questa definisce una nuova relazione, per definizione ancora suriettiva, tra Φ e
l’insieme di tutte le densità elettroniche N:
C :W →Ψ (4.40)
con il codominio, ovvero quando ogni elemento del codominio è immagine di almeno un punto
del dominio.
2 Un’applicazione si dice iniettiva se associa ad un elemento del codominio un solo elemento
del dominio.
76
4.4 Teoria del funzionale densità
(T +W +V )|φs f i = Es f |φs f i
(4.43)
0
(T +W +V )|φs0 f i = Es0 f |φs0 f i
Adesso immaginiamo che, per assurdo, le 4.43 siano soddisfatte dalla stessa fun-
zione d’onda, cioè che:
|φs f i = |φs0 f i (4.44)
I due potenziali esterni differirebbero quindi soltanto per una costante additiva,
e, di conseguenza, sarebbero equivalenti contro l’ipotesi che avevamo assunto.
Abbiamo quindi dimostrato che C associa funzioni d’onda differenti a potenziali
esterni differenti ed è pertanto iniettiva, oltre che suriettiva, e quindi biunivoca.
77
Studio quantistico dei sistemi a molti corpi
H = H 0 +W −W 0 (4.47)
dove:
Z
hφs f |(W −W 0 )|φs f i = d~rφs∗f [w(~r) − w0 (~r)]φs f (~r) =
Z (4.49)
= d~r[w(~r − w0 (~r))]n(~r) = hφs0 f |(W −W 0 )|φs0 f i
quindi i secondi addendi delle equazioni (4.48) sono identici se le due funzioni
d’onda conducono ad uguali densità elettroniche. Allora, sottraendo membro a
membro le due equazioni (4.48) otteniamo l’uguaglianza:
Sfruttando il principio di Ritz (vedi Appendice E), il quale afferma che il valore
d’aspettazione dell’Hamiltoniana sullo stato fondamentale è minore di quello su
qualunque altro stato appartenente allo stesso spazio di Hilbert, ed avendo suppo-
sto che le funzioni d’onda di stato fondamentale sono differenti, dovremmo avere
contemporaneamente:
e
hφs f |H 0 |φs f i − hφs0 f |H 0 |φs0 f i > 0 (4.52)
78
4.4 Teoria del funzionale densità
Figura 4.1: Illustrazione grafica di alcuni aspetti del teorema di Hohenberg e Kohn. Si dimostra
che le relazioni C e D, definite in maniera da garantirne la suriettività (in alto) sono
anche iniettive e quindi biettive e pienamente invertibili.
79
Studio quantistico dei sistemi a molti corpi
prima parte del teorema, un funzionale unico della densità elettronica. D’altra
parte il principio di Ritz garantisce che:
dunque, per la (4.53), il funzionale E[n] assumerà il suo valore minimo in cor-
rispondenza della densità di carica ns f (~r) e tale valore sarà appunto l’energia dello
stato fondamentale. Riscritto in formule otteniamo:
80
4.4 Teoria del funzionale densità
81
Studio quantistico dei sistemi a molti corpi
82
4.4 Teoria del funzionale densità
1
Z
EXC [n] = FHK [n] − d~r0 n(~r)v(~r,~r0 )n(~r0 ) − T 0 [n] (4.68)
2
83
Studio quantistico dei sistemi a molti corpi
δ EXC [n]
vXC ([n],~r) = (4.71)
δ (~r)
Per la prima uguaglianza nella (4.72) abbiamo fatto uso del secondo teorema
di Green e nella sceonda abbiamo sostituito l’equazione (4.61) e la sua com-
plessa coniugata. Osserviamo che, in conseguenza della normalizzazione delle
funzioni d’onda di particella singola, il primo termine dell’ultima uguaglianza è
nullo. Cosı̀, utilizzando la (4.60), troviamo che la variazione del termine cinetico
si riduce a:
N Z
δ T 0 [n] = − ∑ d~rw0 (~r)δ n(~r) (4.73)
i=1
Inserendo quest’ultima relazione nella (4.69) o nella (4.70) e tenendo conto del
fatto che le variazioni δ n(~r) sono arbitrarie, otteniamo finalmente l’equazione che
ci fornisce il potenziale esterno cercato per il sistema non interagente:
Z
0
w (~r) = w(~r) + d~r v(~r,~r0 )n(~r0 ) + vXC ([n],~r) (4.74)
84
4.4 Teoria del funzionale densità
in termini del potenziale esterno e del potenziale colombiano per il sistema reale
oltre che del potenziale di scambio e correlazione. Dalla (4.74) risulta evidente
che il potenziale w0 è esso stesso un funzionale della densità elettronica.
Lo schema di Kohn e Sham consiste nel risolvere l’equazione (4.61) per una
particella singola con la prescrizione (4.74) per w0 . Questa quantità può quindi es-
sere interpretata come un potenziale efficace a particella singola, in quanto ci con-
sente di ottenere le proprietà di un sistema interagente attraverso un’equazione di
Schrödinger a particella singola. In altre parole w0 contiene, risommata, tutta l’in-
formazione relativa alle complicate interazioni a molti corpi presenti nel sistema
reale. Il risultato di questa risomma è che il potenziale w relativo all’interazione
di un elettrone con i nuclei viene indebolito, schermato, dalla presenza della ca-
rica negativa di tutti gli altri elettroni del sistema. Per questa ragione w0 , che
pure dipende localmente da ~r, deve possedere anche una dipendenza funzionale
dal campo della densità, in modo da poter incorporare l’effetto delle interazioni
con tutti gli altri fermioni del sistema. Inoltre l’energia di scambio e correlazione
(4.68) contiene la differenza tra l’energia cinetica del sistema vero e quella del
sistema non interagente (presumibilmente piccola) [66].
Dalla minimizzazione del funzionale dell’energia, tenendo conto del vincolo
sulla densità (4.60), e imponendo che la densità sia ottenuta attraverso un sistema
associato di elettroni non interagenti, otteniamo un’equazione one-electron per gli
orbitali a singola particella detta equazione di Kohn e Sham [68]:
" #
1 2
− ∇i + w0 (~r) φi (~r) = εi φi (~r) (4.75)
2
85
Studio quantistico dei sistemi a molti corpi
Bisogna evidenziare un aspetto molto importante. Gli stati di Kohn e Sham che
si ottengono risolvendo l’equazione (4.75) sono quegli stati che serviranno per la
costruzione del determinante di Slater relativo alla descrizione dello stato fonda-
mentale del sistema di riferimento. Quest’ultimo condivide con il sistema vero
soltanto la densità di carica dello stato fondamentale e nient’altro. Tutto ciò
ci impedisce di dare un significato fisico ai singoli stati di Kohn e Sham che
possono essere visti esclusivamente come uno strumento matematico di ausilio
per la risoluzione del problema originario. Ciò premesso a volte essi vengono
interpretati come stati a singola particella in un campo medio.
Come abbiamo visto nella precedente sezione lo schema di Kohn e Sham non
introduce alcuna approssimazione in quanto, l’errore commesso nella valutazione
del funzionale dell’energia cinetica, viene trasferito e risommato nella grandezza
EXC . Cosı̀ facendo, la teoria del funzionale densità continua ad essere una teoria
esatta.
86
4.4 Teoria del funzionale densità
Adesso è possibile distinguere due casi: un primo caso in cui viene fatta
l’ipotesi che il funzionale di scambio e correlazione sia esatto e corrispondente
alla situazione reale in cui si dispone soltanto di un’approssimazione di tale fun-
zione. Il risultato ottenuto con il funzionale esatto, EXC , coinciderebbe con la
corretta soluzione per lo stato fondamentale dell’equazione di Schrödinger per il
sistema di elettroni interagenti. L’esatta forma del funzionale di scambio e cor-
relazione, però, non è nota ma, anche se lo fosse, probabilmente sarebbe molto
complessa e soprattutto avremmo una dipendenza non-locale dalla densità. La
qualità dei calcoli DFT dipende esclusivamente dalla qualità dell’approssimazione
del funzionale di scambio e correlazione e la risultante energia di scambio e
correlazione.
Sfortunatamente non esistono delle strategie da seguire per ottenere e miglio-
rare la forma del funzionale EXC , comunque, ci sono alcuni vincoli fisici che un
funzionale della densità deve rispettare. Queste regole possono essere usate come
linee guida per lo sviluppo del funzionale di scambio e correlazione. Tra queste ri-
cordiamo, ad esempio, le regole di somma per le buche di scambio e correlazione
(una buca di scambio e correlazione descrive lo svuotamento rinormalizzato del-
la densità elettronica nello spazio quando, piuttosto che elettroni indipendenti, si
considerano elettroni interagenti).
Una delle approssimazioni più usate per il potenziale di scambio e correlazione
è la Local Density Approximation (LDA) [68]. Qui il funzionale EXC viene calco-
lato assumendo che, l’energia di scambio e correlazione di un punto dello spazio
con una data densità elettronica è uguale all’energia di scambio e correlazione,
εXC [n(~r)], di un gas omogeneo di elettroni alla stessa densità. Integrando nello
spazio si ha
Z
EXC = εXC [n(~r)]n(~r)d~r =
Z Z (4.76)
= εX [n(~r)]n(~r)d~r + εc [n(~r)]n(~r)d~r
87
Studio quantistico dei sistemi a molti corpi
La parte di scambio εX [n(~r)] viene espressa sulla base del risultato teorico di Bloch
e Dirac:
r
3 3 3n(~r)
εX [n(~r)] = (4.77)
4 π
88
4.5 Dinamica molecolare ab initio
B3LY P
EXC = (1 − a)ExLDA + aEXexact + b∆EXB88 + (1 − c)ECLDA + cECLY P (4.80)
89
Studio quantistico dei sistemi a molti corpi
ed elettroni e fondato tanto sulla teoria del funzionale della densità (DFT), quanto
sulla dinamica molecolare classica (MD).
Abbiamo già visto che la MD richiede un modello per il potenziale relativo al-
l’interazione efficace classica tra gli atomi. Il metodo di Car e Parrinello (CP) è, al
contrario, in grado di determinare ab initio, ovvero a partire da una descrizione mi-
croscopica in termini di interazioni fondamentali, le interazioni atomiche efficaci.
Nella CP vengono assunte la validità dell’approssimazione di Born-Oppenheimer
e la possibilità di descrivere in maniera classica i gradi di libertà relativi ai nu-
clei. Di conseguenza si ottengono equazioni classiche per il moto dei nuclei che
contengono l’accoppiamento con gli elettroni attraverso un termine di forza di
Hellmann-Feynmann. Gli elettroni, invece, vengono descritti nello schema di
Kohn e Sham in termini di funzioni d’onda a particella singola, ψi , con la densità
di carica fornita dalla 4.60:
N
n(~r) = ∑ |ψi (~r)|2 (4.81)
i=1
con {RI } sono indicate le coordinate nucleari e {αv } sono tutti i possibili vincoli
esterni imposti sul sistema, come il volume, etc. Il funzionale U contiene la re-
pulsione internucleare coulombiana e l’energia potenziale elettronica efficace, che
include i nuclei esterni e i contributi di scambio e correlazione.
90
4.5 Dinamica molecolare ab initio
La lagrangiana (4.82) genera una dinamica per {ψi }, {RI } e {αv } le cui equazioni
del moto sono
δE
µ ψ̈(~r,t) = − + Λik ψk (~r,t) (4.85)
δ ψi∗ (~r,t) ∑
k
1 1 1
Z
K = Kel + Kion = µ ∑ |ψ̇i |2 d 3~r + ∑ MI Ṙ2I + ∑ µv α̇v2 (4.88)
2 i 2 I 2 v
91
Studio quantistico dei sistemi a molti corpi
La DSA è una delle possibili realizzazioni del cosidetto metodo del Simuleated
Annealing che serve per risolvere problemi complessi di minimizzazione, in cui
una funzione oggetto O({β }) è minimizzata rispetto al parametro {β }, generando
una successione di configurazioni {β } con una distribuzione di probabilità con
peso alla Boltzmann attraverso una procedura Monte Carlo. Per T → 0 lo stato
con il più basso O({β }) viene raggiunto senza che il sistema resti intrappolato
all’interno di uno stato metastabile.
92
4.6 Metodi ibridi QM/MM
Spesso quando il sistema che vogliamo studiare è una molecola biologica co-
me un peptide è necessario tener conto anche dell’interazione con il solvente o
con altre molecole biologiche. In questo caso il sistema è troppo grande perchè
si possa trattare con un approccio quantomeccanico. Spesso si tenta di simulare
il sistema, adoperando un modello ridotto dello stesso, per esempio simulando
esclusivamente il sito catalitico o di legame ad un metallo, ed escludendo il resto.
Ciò però non tiene in considerazione le interazioni con il resto del sistema che
possono determinare cambiamenti strutturali importanti, per esempio quando c’è
una interazione con il substrato, o tra proteine e solvente. Potremmo pensare di
utilizzare la dinamica molecolare classica per studiare tali sistemi nella loro in-
terezza, ma, come abbiamo visto, questi metodi si basano sulla conoscenza e sulla
giusta parametrizzazione del ff. Inoltre con i metodi classici non siamo in grado di
descrivere neppure in maniera approssimata la rottura e la formazione di un lega-
me chimico. Per questo motivo negli ultimi anni sono stati sviluppati metodi ibridi
detti QM/MM (Quantum Mechanics/Molecular Mechanics) che suddividono il si-
stema in due parti trattando una parte mediante metodi QM e una parte mediante
una metodologia classica MM.
Lo sviluppo dei metodi ibridi QM/MM, che ha inizio con gli studi di Warshel e
Levitt [73] e di Singh e Kollman [74], segue come linea generale l’idea che un si-
93
Studio quantistico dei sistemi a molti corpi
94
4.6 Metodi ibridi QM/MM
non sono connesse da legami covalenti; mentre nel caso in cui le due regioni siano
connesse da legami covalenti la situazione si complica notevolmente.
Il livello più basso d’interazione è detto mechanical embedding. In questo
caso, solo le energie steriche e di bond delle due regioni sono incluse nei termini
d’interazione, cioè, gli atomi QM sentono forze aggiuntive generate dalla struttura
MM, e vice versa, ma non c’è nessuna interazione tra le parti elettroniche delle
due regioni. Ad ogni atomo QM sono associati dei parametri di Van der Waals che
sono utilizzati in un’espressione simile a quella già vista nel caso della dinamica
molecolare classica, (modellizzata attraverso un potenziale di tipo LJ), cioè:
" !12 !6 #
N M σi j σi j
HQM/MM = ∑ ∑ 4εi j − (4.92)
i j |~Ri −~r j | |~Ri −~r j |
dove l’indice i è associato ai nuclei della regione QM, j agli atomi della regione
MM, ~Ri è il vettore posizione associato all’i-esimo nucleo QM ed ~r j è il vettore
posizione associato al j-esimo atomo MM. Se le due regioni sono unite da legami
chimici ci sono dei termini aggiuntivi corrispondenti alle interazioni di stretching
e bending. Il modello di accoppiamento meccanico è un utile livello di approssi-
mazione quando le funzioni d’onda della regione QM non sono modificate dai
cambiamenti nella regione MM.
Il livello successivo di approssimazione è detto electronic embedding, dove gli
atomi nella regione MM possono polarizzare la regione QM. Le cariche parziali
sugli atomi MM possono essere inserite nell’Hamiltoniana QM e gli atomi QM
sentono il potenziale elettrico dovuto a tutti gli atomi MM, cioè:
n M qj N M Zi q j
HQM/MM = − ∑ ∑ +∑∑ (4.93)
l j |~xl −~r j |
~
i j |Ri −~r j |
dove q j è la carica parziale associata agli atomi MM, ~xl è il vettore posizione as-
sociato all’l-esimo elettrone QM e Zi è la carica associata all’i-esimo nucleo QM.
Questo livello di approssimazione tiene conto di quei cambiamenti geometrici nel-
95
Studio quantistico dei sistemi a molti corpi
la regione MM che hanno effetto sulla regione QM, cioè accoppiano la funzione
d’onda QM alla geometria MM.
Un ulteriore termine può essere incluso assumendo che gli atomi QM polar-
izzano anche la regione MM, cioè il campo elettrico generato dalla regione QM
influenzi il momento elettrico MM. Questo termine prende il nome di polarizable
embedding, e richiede naturalmente uno sforzo computazionale aggiuntivo, dato
che necessita di una procedura doppiamente iterativa dovuta al fatto che tanto il
campo elettrico della regione QM quanto quello MM devono essere determinati
in maniera auto consistente.
La maggior parte dei metodi QM/MM implementati a tutt’oggi utilizzano l’ap-
prossimazione electronic embedding; l’unica eccezione è il metodo dei frammenti
effettivi, dove sia il momento di quadrupolo che le polarizzabilità sono incluse per
gli atomi MM [51].
96
4.6 Metodi ibridi QM/MM
97
Studio quantistico dei sistemi a molti corpi
L’atomo di idrogeno è il più utilizzato a questo scopo che viene inserito lungo
il legame QM/MM che è stato rotto. Possono essere utilizzati anche altri atomi
o gruppi funzionali come gli alogeni o i gruppi metile. I link atoms sono trattati
in maniera esatta nei calcoli QM ma non è permesso loro di interagire con gli
atomi MM (è come se risultassero invisibili agli atomi MM), l’energia e le forze
risultano dai termini energetici di legame. I termini elettrostatici e di Van der
Waals che si considerano di non-legame sono trattati separatamente:
98
4.6 Metodi ibridi QM/MM
99
Studio quantistico dei sistemi a molti corpi
Su questa linea sono stati proposti vari metodi il primo dei quali combina
una tecnica basata su un approccio orbitale (quantistico) con una tecnica MM
(IMOMM) [80]. In un secondo metodo sono combinate due tecniche differen-
ti basate su un approccio quantistico a vari livelli (IMOMO) [81]. La tecnica
ONIOM (Our OwN n-layered Integrated MO and MM) [82] è una combinazione
di queste tecniche in cui il sistema è suddiviso in due o più strati trattati a vari
livelli di approssimazione (Fig.4.6). L’idea base [83] di questa tecnica può essere
spiegata in modo semplice considerando il caso in cui il sistema viene diviso in
due e tre regioni. La Fig.4.7 fornisce uno schema del metodo di estrapolazione
dell’energia. Lo scopo è quello di descrivere il sistema reale al livello più alto
Figura 4.7: Schema di estrapolazione per un sistema partizionato in due e tre strati
100
4.6 Metodi ibridi QM/MM
della teoria (quantistico), cioè il punto 4 nel caso che il sistema sia suddiviso in
due strati oppure il punto 9 se il sistema è suddiviso in tre strati. Nel caso di due
strati, l’energia totale estrapolata EONIOM è definita come:
dove E3 è l’energia dell’intero sistema (real) calcolata al livello più basso (MM)
ed E1 e E2 sono le energie del sistema modello determinate al livello basso (MM)
e al livello alto (QM), rispettivamente. Per un sistema suddiviso in tre strati,
l’espressione per l’energia EONIOM è data da:
link atoms. Questi ultimi sono presenti solo nel sistema modello e le loro coor-
dinate sono descritte dal vettore ~R2 . Nel sistema reale sono sostituiti dagli atomi
descritti dal vettore ~R3 . Gli atomi che appartengono allo strato esterno e non sono
sostituiti dai link atoms sono contenuti nel set 4 e sono descritti dal vettore ~R4 .
La geometria del sistema reale è quindi descritta dai vettori ~R1 , ~R3 , ~R4 che sono
variabili indipendenti per l’energia ONIOM:
Per generare il sistema modello, descritto da ~R1 e dai link atoms ~R2 , si definisce
~R2 in funzione di ~R1 ed ~R3 :
~R2 = f (~R1 , ~R3 ) (4.97)
101
Studio quantistico dei sistemi a molti corpi
La forma funzionale esplicita della dipendenza di ~R2 può essere scelta in maniera
del tutto arbitraria. Tuttavia, considerando il fatto che i link atoms vengono in-
trodotti per mimare i corrispondenti legami covalenti del sistema reale, dobbiamo
in qualche modo simulare i movimenti degli atomi che vengono sostituiti. Lo
schema d’accoppiamento è il seguente: sia A un atomo appartenente al set 1 e B
un atomo appartenente al set 3, ed H un atomo appartenente al set 2 (Fig.4.8) posto
lungo l’asse di legame A-B. In termini delle coordinate interne scegliamo le stesse
coordinate di bond e gli stessi angoli e diedri per gli atomi del set 2 equivalenti a
quelli del set 3. Quindi, nel sistema modello i link atoms sono sempre allineati al
vettore di bond del sistema reale. Per l’esatta posizione~r2 di un singolo atomo H
lungo un bond A-B (~r3 −~r1 ), si introduce un fattore di scala fisso g, ottenendo:
In questo modo se la distanza di legame A-B |~r3 −~r1 | cambia durante la simula-
zione, anche la distanza di legame A-H |~r2 −~r1 | cambierà. La scelta del parametro
g dipende in maniera cruciale dalla natura del legame spezzato, e può dipendere
anche dal livello di teoria usato per descrivere i due strati che sono connessi dal
legame tagliato.
102
Capitolo 5
Come abbiamo visto nel Cap.2 le placche, presenti nelle cellule dei malati
di AD, sembrano essere un “serbatoio metallico”, viste le alte concentrazioni di
Cu2+ , Zn2+ , e Fe3+ . Nel caso dell’Aβ gli ioni Zn2+ e Cu2+ sembrano avere un
ruolo nel processo di misfolding e nella formazione delle fibrille.
Focalizzeremo la nostra attenzione esclusivamente sull’interazione tra Zn2+
e Aβ , in quanto i dati sperimentali attualmente disponibili non forniscono una
chiara interpretazione del suo ruolo. Alcuni sembrano indicare un possibile ruolo
dello ione Zn2+ nella formazione degli aggregati iniziali [42, 43], mentre altri nel
processo del misfolding [44].
L’obiettivo è quello di analizzare la coordinazione dello ione Zn2+ , in partico-
lare, l’interazione con il frammento peptidico Aβ1−16 , identificato come il minimo
frammento coinvolto nel legame con il metallo [43], la cui sequenza è:
NH3+ -DAEFRHDSGYEVHHQL-COO−
104
5.1 Operazioni preliminari
Il nostro primo passo è stato quello di simulare il sistema Zn2+ -Aβ1−16 sen-
za solvatare il sistema utilizzando l’algoritmo di dinamica molecolare classica
implementato nel programma open source Gromacs 3.3.3 [84]. Per la modelliz-
zazione del sistema è stata usata la struttura NMR [44], disponibile sul Protein
Data Bank (PDB) (ID: 1ZE9). Inoltre, sulla base di dati precedentemente ottenuti
attraverso un ottimizzazione geometrica DFT, abbiamo modificato il sito di coor-
dinazione metallico, servendoci del programma GaussView, generando la sfera di
Atomi ri (Å)
Nδ (His-6) 2.11
O1(Glu-11) 2.10
O2(Glu-11) 2.10
Nε(His-13) 2.15
Nδ (His-14) 2.29
Tabella 5.1: Distanza dallo Zn2+ degli atomi che si trovano inizialmente nella sfera di
coordinazione del metallo
105
Simulazioni di dinamica molecolare classica
• Nel primo set di simulazioni non abbiamo inserito i legami chimici dello
ione metallico con i liganti quindi l’unica interazione tra lo Zn2+ e il peptide
è data dalle interazioni nonbonded (Van der Walls ed elettrostatiche).
• Nel secondo set abbiamo inserito invece i legami chimici (nel file di topolo-
gia (.top)). I parametri di bond (Tabella 5.2) ed angles (Tabella 5.3) per
gli atomi che noi riteniamo “legati” allo ione metallico, suggeriti da calcoli
quantomeccanici DFT, effettuati in precedenza dalla dott.ssa Marino presso
l’Università della Calabria, e dalla referenza [88].
lungo la sequenza di Aβ .
3 σ si ottiene attraverso una media aritmetica dei σ dei singoli atomi, ed ε si ottiene
ij ii ij
106
5.1 Operazioni preliminari
Tabella 5.2: Valori della distanza d’equilibrio e costante elastica per gli atomi legati allo Zn2+
Tabella 5.3: Valori dell’angolo d’equilibrio e costante di forza angolare per gli atomi legati allo
Zn2+
107
Simulazioni di dinamica molecolare classica
Tabella 5.5: Condizioni simulative utilizzate per la dinamica molecolare in assenza di solvente
Il secondo passo è stato quello di simulare, sulla base dei dati raccolti in prece-
denza, il sistema Zn2+ -Aβ1−16 in presenza di acqua, utilizzando il force field
CHARMM27. Inoltre abbiamo modificato il force field aggiungendo i potenziali
di bond (Tabella 5.2), angle (Tabella 5.3) ed i parametri di LJ (Tabella 5.4). Il
sistema è stato solvatato in un box cubico con 2517 molecole d’acqua TIP3P [91].
4I box utilizzati sono di forma cubica.
108
5.1 Operazioni preliminari
Tabella 5.6: Condizioni simulative utilizzate per la dinamica molecolare in presenza di solvente
109
Simulazioni di dinamica molecolare classica
5.2 Risultati
Come già visto abbiamo provato a simulare il sistema, oltre che con vari force
field, con tecniche differenti.
Prova 1: Nella prima serie di prove abbiamo escluso legami chimici con lo
ione metallico considerando solo interazioni elettrostatiche e LJ fra lo ione metal-
lico e atomi del peptide. Abbiamo testato i tre force field: Gromos96, OPLS e
CHARMM27.
Gromos96: Abbiamo monitorato gli atomi presenti nella sfera di coordina-
zione del metallo (Fig.5.1), durante tutta la fase di termalizzazione e durante la
dinamica a 300K. Si vede che lo ione metallico rimane nella coordinazione ini-
ziale (Tabella 5.1), per quasi tutta la durata della simulazione a 200K. Verso la
Figura 5.1: Grafico delle distanze di vari Figura 5.2: Grafico dell’energia totale in
atomi dallo Zn2+ in funzione funzione del tempo (Prova 1 -
del tempo (Prova 1 - Gro- Gromos96)
mos96)
110
5.2 Risultati
Figura 5.3: Grafico delle distanze di vari Figura 5.4: Grafico dell’energia totale in
atomi dallo Zn2+ in funzione funzione del tempo (Prova 1 -
del tempo (Prova 1 - OPLS) OPLS)
111
Simulazioni di dinamica molecolare classica
essendo partito da una coordinazione penta. Si noti che i due atomi Nδ apparte-
nenti agli amminoacidi His-6 e His-14 si allontanano dallo ione metallico, e in
certi momenti sono (Fig.5.3) al di là della sfera di coordinazione.
CHARMM27: Prendendo spunto dal lavoro di Sakharov et al. [87] abbiamo
provato a simulare classicamente il sistema Zn2+ -Aβ1−16 usando il programma
Gromacs 3.3.3 abbinato al force field CHARMM27. CHARMM27 non è pre-
sente nel pacchetto del programma Gromacs 3.3.3. La versione da noi utilizza-
ta è stata scaricata dal sito di uno degli autori. Per prima cosa abbiamo dovu-
to modificare CHARMM27 in modo da renderlo consistente con il program-
ma Gromacs 3.3.3. Utilizzando uno script scritto in PERL, abbiamo genera-
to, a partire da un file di parametri di CHARMM (par all27 lipid.prm), i due
file ffcharmmbon.itp e ffcharmmnb.itp che verrano utilizzati dal programma Gro-
macs come nuovi force field. Una volta generato il file.top attraverso il comando
pdb2gmx, e dopo aver scelto come force field CHARMM27, bisogna utilizzare
un programma che modifica il file.top ottenuto in precedenza, aggiungendo alcuni
parametri. Ottenuto un nuovo file.top, siamo pronti a lanciare la nostra simula-
Figura 5.5: Grafico delle distanze di Figura 5.6: Grafico dell’energia totale in
vari atomi dallo Zn2+ in funzione del tempo (Prova 1 -
funzione del tempo (Prova 1 - CHARMM27)
CHARMM27)
112
5.2 Risultati
zione. Notiamo (Fig.5.5) come in questo caso lo ione metallico tenda a rimanere
penta-coordinato e l’energia totale oscilli attorno ad un valore costante alle varie
temperature (Fig.5.6).
Si fa notare che in tutti e tre i casi, nonostante non siano imposti legami chimici
per lo ione metallico, la coordinazione del metallo non tende mai a diminuire,
anzi in due casi aumenta. E’ possibile che questo sia dovuto ad una sovrastima
della carica netta “classica” (+2) dello ione. Inoltre in tutti e tre i force field la
parametrizzazione dello ione metallico (nonbonded) viene effettuata su modelli
Zn2+ -H2 O dove lo ione metallico è sempre esacoordinato.
La Fig.5.7 mostra le deviazioni quadratiche medie (RMSD) [92] del com-
plesso Zn2+ -Aβ1−16 rispetto alla struttura di partenza. Notiamo nella Fig.5.7b
due brusche variazioni fra i 1000ps e i 1200ps, in coincidenza dell’avvicinamen-
to della coda carbossilica (Lys-16). Questo repentino avvicinamento provoca un
profondo mutamento di tutta la struttura del complesso, messo in evidenza dalle
RMSD. Anche la Fig.5.7a mostra un cambiamento, meno brusco del precedente,
in coincidenza dell’entrata all’interno della sfera di coordinazione dell’atomo di
O (His-6), che avviene a ∼ 1400ps. Notiamo che CHARMM27 ha una RMSD
Figura 5.7: RMSD del sistema Zn2+ -Aβ1−16 rispetto alla struttura di partenza, dove con il simbo-
lo h i indichiamo la media e con σ la sua deviazione standard. Prova 1: (a) Gromos96;
(b) OPLS; (c) CHARMM27
113
Simulazioni di dinamica molecolare classica
più bassa rispetto agli altri force field. Ciò potrebbe indicarci che la parametriz-
zazione di CHARMM27 sia più idonea rispetto agli altri, almeno nel caso preso
in esame (metallo-proteina).
Infine, nella Tabella 5.7, riportiamo le distanze medie hri lungo la traiettoria
e le deviazioni standard σ degli atomi che si trovano all’interno della sfera di
coordinazione dello Zn2+ al termine della simulazione.
Atomi r(Å) σ (Å)
Gromos96
Nδ (His-6) 2.09 0.06
O1(Glu-11) 2.10 0.05
O2(Glu-11) 2.10 0.05
Nε(His-13) 2.10 0.06
Nδ (His-14) 2.08 0.06
O(His-6) 2.08 0.08
OPLS
Nδ (His-6) 2.27 0.10
O1(Glu-11) 1.97 0.05
O2(Glu-11) 1.97 0.04
Nδ (His-14) 2.34 0.15
O1(Lys-16) 1.96 0.04
O2(Lys-16) 1.98 0.05
CHARMM27
Nδ (His-6) 2.17 0.05
O1(Glu-11) 1.97 0.04
O2(Glu-11) 1.99 0.04
Nε(His-13) 2.17 0.06
Nδ (His-14) 2.15 0.05
Tabella 5.7: Distanze medie e deviazioni standard degli atomi che si trovano all’interno della sfera
di coordinazione dello Zn2+ al termine della Prova 1
114
5.2 Risultati
Figura 5.8: Grafico delle distanze di Figura 5.9: Grafico dell’energia totale in
vari atomi dallo Zn2+ in funzione del tempo (Prova 2 -
funzione del tempo (Prova 2 - Gromos96)
Gromos96)
115
Simulazioni di dinamica molecolare classica
Figura 5.10: Grafico delle distanze di vari Figura 5.11: Grafico dell’energia totale in
atomi dallo Zn2+ in funzione funzione del tempo (Prova 2 -
del tempo (Prova 2 - OPLS) OPLS)
OPLS: Si vede (Fig.5.10) che nei primi ps di dinamica gli atomi O1 e O2,
che fanno parte della coda carbossilica (Lys-16), si portano all’interno della sfera
di coordinazione. Contemporaneamente i quattro amminoacidi legati inizialmente
allo ione metallico si allontanano, in particolare gli atomi Nδ appartenenti all’His-
6 e His-14. Inoltre l’atomo O2 (Glu-11) esce durante la fase di minimizzazione
dalla sfera di coordinazione, contrariamente a quanto visto nella Prova 1. Il sito di
legame metallico raggiunge la coordinazione esa già a 200K e la mantiene anche
a 300K. L’energia totale ha un regolare andamento (Fig.5.11).
CHARMM27: Già durante la fase di minimizzazione l’atomo O appartenente
al backbone dell’amminoacido Gly-9 entra dentro la sfera di coordinazione dello
Zn2+ . Inoltre, le oscillazioni di quest’atomo attorno alla distanza media sono più
grandi degli altri atomi, aumentando in maniera progressiva all’aumentare della
temperatura (Tabella 5.8). Anche in questo caso l’atomo O2 (Glu-11) presente
nella configurazione iniziale all’interno della sfera di coordinazione dello ione
metallico, si allontana durante la minimizzazione. Lo ione metallico resta stabil-
mente penta-coordinato (Fig.5.12) durante tutta la simulazione. L’energia totale
ha un regolare andamento (Fig.5.13).
116
5.2 Risultati
Figura 5.12: Grafico delle distanze di vari Figura 5.13: Grafico dell’energia totale in
atomi dallo Zn2+ in fun- funzione del tempo (Prova 2 -
zione del tempo (Prova 2 - CHARMM27)
CHARMM27)
La Fig.5.14 mostra l’RMSD per tutti i force field utilizzati. Notiamo che i
valori medi sono più alti rispetto alla prima prova. La Fig.5.14a mostra l’RMSD
relativo a Gromos96. Questo ha un andamento particolarmente oscillante ed inol-
117
Simulazioni di dinamica molecolare classica
Figura 5.14: RMSD del sistema Zn2+ -Aβ1−16 rispetto alla struttura di partenza, dove con il sim-
bolo h i indichiamo la media e con σ la sua deviazione standard. Prova 2: (a)
Gromos96; (b) OPLS; (c) CHARMM27
118
5.2 Risultati
Tabella 5.8: Distanze medie e deviazioni standard degli atomi che si trovano all’interno della sfera
di coordinazione dello Zn2+ al termine della Prova 2
119
Simulazioni di dinamica molecolare classica
Figura 5.15: Grafico delle distanze di vari Figura 5.16: Grafico dell’energia totale in
atomi dallo Zn2+ in fun- funzione del tempo (Prova 3 -
zione del tempo (Prova 3 - Gromos96)
Gromos96)
OPLS: Osserviamo che fin dalle prime fasi della simulazione (Fig.5.17) tutti
gli atomi che nella configurazione iniziale si trovavano all’interno della sfera di
coordinazione si allontanano sensibilmente. In particolare i due atomi Nδ delle
(His-6 e -14) oscillano attorno ad una distanza media maggiore di 2.5Å. All’au-
mentare della temperatura le oscillazioni attorno al valore medio di questi due
120
5.2 Risultati
Figura 5.17: Grafico delle distanze di vari Figura 5.18: Grafico dell’energia totale in
atomi dallo Zn2+ in funzione funzione del tempo (Prova 3 -
del tempo (Prova 3 - OPLS) OPLS)
121
Simulazioni di dinamica molecolare classica
Figura 5.19: Grafico delle distanze di vari Figura 5.20: Grafico dell’energia totale in
atomi dallo Zn2+ in fun- funzione del tempo (Prova 3 -
zione del tempo (Prova 3 - CHARMM27)
CHARMM27)
teriore segnale che la direzione seguita, cioè la modifica dei parametri nonbonded,
possa essere quella giusta.
Figura 5.21: RMSD del sistema Zn2+ -Aβ1−16 rispetto alla struttura di partenza, dove con il sim-
bolo h i indichiamo la media e con σ la sua deviazione standard. Prova 3: (a)
Gromos96; (b) OPLS; (c) CHARMM27
122
5.2 Risultati
possiamo notare che le distanze medie sono molto grandi e gli atomi Nδ degli
amminoacidi His-6 e His-14 sono fuori dalla sfera di coordinazione. Questo ci
Atomi r(Å) σ (Å)
Gromos96
Nδ (His-6) 2.31 0.06
O1(Glu-11) 2.26 0.05
Nε(His-13) 2.28 0.06
Nδ (His-14) 2.33 0.06
O1(Glu-3) 1.69 0.04
O2(Glu-3) 1.69 0.04
OPLS
Nδ (His-6) 2.55 0.07
O1(Glu-11) 2.32 0.06
Nε(His-14) 2.33 0.07
Nδ (His-14) 2.54 0.07
CHARMM27
Nδ (His-6) 2.24 0.05
O1(Glu-11) 2.15 0.05
Nε(His-13) 2.14 0.05
Nδ (His-14) 2.20 0.05
O(Gly-9) 1.74 0.08
Tabella 5.9: Distanze medie e deviazioni standard degli atomi che si trovano all’interno della sfera
di coordinazione dello Zn2+ al termine della Prova 3
induce a pensare che i nuovi parametri di LJ da noi inseriti siano troppo piccoli,
rispetto a quelli presenti nel force field OPLS. Questo da un lato impedisce ad altri
atomi di entrare nella sfera di coordinazione dello Zn2+ , ma probabilmente tiene
lontani quelli presenti nella configurazione di partenza. Possiamo ipotizzare che
OPLS sia costruito in maniera fortemente consistente ed un cambiamento in uno
dei parametri altera in modo evidente la dinamica generata con questo force field.
123
Simulazioni di dinamica molecolare classica
124
5.2 Risultati
totale del sistema a 100K lentamente tenda a stabilizzarsi. Per questo motivo la
dinamica è stata leggermente allungata.
Figura 5.22: Grafico delle distanze durante Figura 5.23: Grafico delle distanze di vari
i primi 20ps di vari atomi dal- atomi dallo Zn2+ in funzione
lo Zn2+ in funzione del tempo del tempo a 100K
a 100K
125
Simulazioni di dinamica molecolare classica
Figura 5.25: Grafico delle distanze di vari Figura 5.26: Grafico dell’energia totale in
atomi dallo Zn2+ in funzione funzione del tempo a 200K
del tempo a 200K
La Fig.5.30 mostra l’RMSD del sistema Zn2+ -Aβ1−16 durante tutta la dina-
mica di 9ns a 300K. Il valore molto alto evidenzia i notevoli cambiamenti che la
Figura 5.27: Grafico delle distanze fra Figura 5.28: Grafico delle distanze di vari
4120 e 4140ps di vari ato- atomi dallo Zn2+ in funzione
mi dallo Zn2+ in funzione del del tempo a 300K
tempo a 300K
126
5.2 Risultati
Figura 5.29: Grafico dell’energia totale in Figura 5.30: RMSD del sistema Zn2+ -
funzione del tempo a 300K Aβ1−16 rispetto alla strut-
tura di partenza durante la
dinamica a 300K
127
Simulazioni di dinamica molecolare classica
Atomi r(Å) σ (Å)
Tabella 5.10: Distanze medie e deviazioni standard degli atomi che si trovano all’interno della
sfera di coordinazione dello Zn2+ al termine della simulazione
5.3 Conclusioni
Le simulazioni di dinamica molecolare classica descritte in questo capitolo
sono state necessarie alla costruzione di modelli da sottoporre all’ottimizzazione
geometrica ONIOM e, in un lavoro successivo a quello presente in questa tesi di
laurea, a simulazioni CP.
Abbiamo confrontato tre force field ma non ci siamo preoccupati di fornire
dei parametri più realistici per le cariche parziali presenti sullo ione metallico
e sui suoi leganti. Analizzando i dati classici abbiamo però notato che le inte-
razioni elettrostatiche sono quelle che hanno maggiormente contribuito a deter-
minare le configurazioni finali dei nostri modelli, fatto che determina un bias per
le simulazioni quantistiche successive.
128
Capitolo 6
Simulazioni QM/MM-ONIOM
• macromolecole biologiche come gli enzimi, il cui sito attivo viene studiato
130
6.1 Operazioni preliminari
Come già detto nel Cap.5, prima di effettuare un’ottimizzazione ONIOM ab-
biamo estratto alla temperatura di 300K alcune configurazioni a bassa energia
ottenute dalle simulazioni di dinamica molecolare classica svolte utilizzando il ff
131
Simulazioni QM/MM-ONIOM
Atomi ri (Å)
Nδ (His-6) 2.26
O(Gly-9) 1.77
O1(Glu-11) 2.16
Nε(His-13) 2.16
Nδ (His-14) 2.23
Tabella 6.2: Distanza dallo Zn2+ degli atomi che si trovano nella sfera di coordinazione del
metallo dopo la fase di minimizzazione classica
CHARMM27 con l’inserimento dei parametri di bond ed angle per lo ione metal-
lico e la modifica dei parametri di LJ. Queste configurazioni sono state minimiz-
zate classicamente utilizzando l’algoritmo CG, i cui parametri simulativi sono
riassunti nella Tabella 6.1.
132
6.1 Operazioni preliminari
• lo spostamento spaziale fra uno step e quello successivo deve essere più
piccolo di un valore di soglia di 0.0018.
Come si vede la differenza di energia tra due step successivi non è un esplici-
to criterio di convergenza. La presenza di quattro distinti criteri di convergenza
previene un’identificazione prematura del minimo. Nel caso dell’ottimizzazione
ONIOM oltre ai quattro criteri di convergenza appena citati si possono aggiun-
gere due nuovi criteri, relativi alle forze, legati alla convergenza della regione
MM. I valori di soglia sono identici a quelli già visti in precedenza. Il termine
d’interazione HQM/MM è calcolato al livello mechanical embedding.
Prima di procedere all’ottimizzazione geometrica abbiamo separato le due re-
gioni utilizzando il software GaussView [93] ed inserito dei link atoms nelle re-
gioni di confine in cui sono stati spezzati dei legami covalenti. Come link atoms
133
Simulazioni QM/MM-ONIOM
134
6.1 Operazioni preliminari
Tabella 6.3: Condizioni simulative utilizzate per la minimizzazione in presenza del solvente
Atomi ri (Å)
Nδ (His-6) 2.521
O(Gly-9) 1.82
O1(Glu-11) 2.26
Nε(His-13) 2.40
Nδ (His-14) 2.43
O(Sol-661) 1.78
O(Sol-774) 1.75
Tabella 6.4: Distanza dallo Zn2+ degli atomi che si trovano nella sfera di coordinazione del
metallo dopo la fase di minimizzazione classica
135
Simulazioni QM/MM-ONIOM
si trovano nella sfera di coordinazione del metallo alla fine del calcolo classico.
Abbiamo escluso dal sistema le restanti molecole di solvente per due motivi; (i) la
maggior parte di esse al termine delle simulazioni di dinamica molecolare classi-
ca e della fase di minimizzazione si trovavano molto distanti dal sito metallico e,
quindi, possiamo supporre non influenzino la coordinazione dello ione metallico;
(ii) i tempi di calcolo ONIOM si dilatano notevolmente in presenza di acqua.
6.2 Risultati
Come già detto per prima cosa abbiamo ottimizzato il sistema Zn2+ -Aβ1−16 in
assenza d’acqua. Il calcolo, dopo 125 cicli, ha soddisfatto tutti e sei i criteri di con-
vergenza. L’energia del sistema ottimizzato è uguale a -1461.167a.u.. Abbiamo
136
6.2 Risultati
Figura 6.3: Grafico delle distanze di vari atomi dallo Zn2+ durante la fase d’ottimizzazione
ONIOM
Figura 6.4: Configurazione finale del siste- Figura 6.5: Ingrandimento del sito metalli-
ma Zn2+ -Aβ1−16 co
137
Simulazioni QM/MM-ONIOM
Atomi ϑ (deg)
La Fig.6.6 mostra l’RMSD del complesso Zn2+ -Aβ1−16 rispetto alla strut-
tura di partenza. La Fig.6.6a rappresenta l’RMSD dell’intero sistema, mentre le
Figg.6.5b-c rappresentano l’RMSD delle regioni QM ed MM rispettivamente. Si
osserva che l’andamento dell’RMSD nelle due regioni è sostanzialmente identico.
Figura 6.6: RMSD del sistema Zn2+ -Aβ1−16 rispetto alla struttura di partenza. (a) Sistema
ONIOM; (b) Regione QM; (c) Regione MM
Abbiamo effettuato un analisi delle cariche parziali usando la tecnica del Natu-
ral Bond Orbital (NBO) (Appendice I) per vedere come queste si sono distribuite
138
6.2 Risultati
Atomi q(e)
Zn2+ 1.306
Nδ (His-6) -0.644
O1(Glu-11) -0.745
O2(Glu-11) -0.822
Nε(His-13) -0.630
Nδ (His-14) -0.655
Tabella 6.6: Cariche parziali dello ione metallico e dei suoi liganti
tra i vari atomi presenti nel sito metallico. Come si vede nella Tabella 6.6 lo Zn,
che inizialmente aveva una carica pari a +2, ha ceduto parte della sua carica agli
atomi circostanti.
In conclusione possiamo osservare che lo ione metallico, inizialmente penta-
coordinato, dopo l’ottimizzazione ONIOM ha una struttura di tetraedro distor-
to. Nella Tabella 6.7 riportiamo le distanze finali dallo Zn2+ degli atomi che si
trovano all’interno della sfera di coordinazione.
Atomi r(Å)
Nδ (His-6) 2.09
O2(Glu-11) 1.97
Nε(His-13) 2.15
Nδ (His-14) 2.09
Tabella 6.7: Distanze degli atomi che si trovano nella sfera di coordinazione dello Zn2+ al termine
dell’ottimizzazione ONIOM
139
Simulazioni QM/MM-ONIOM
al metallo. Solo due criteri di convergenza (relativi alla parte MM) dei sei sono
stati soddisfatti. Nonostante ciò possiamo supporre che il sistema, dopo 287 cicli,
sia ottimizzato a pieno. Infatti osservando le variazioni dell’energia totale noti-
amo che, negli ultimi dieci cicli, sono molto piccole e modificano l’energia alla
quarta cifra decimale. Inoltre l’algoritmo implementato in Gaussian03 predice
quali saranno le variazioni in energia fra un ciclo e l’altro, e questo mostra che da
un certo ciclo in poi le variazioni sono molto piccole. Per questo motivo possia-
mo pensare che il sistema si trova in un minimo globale sulla superficie d’ener-
gia potenziale (PES). L’energia del sistema ottimizzato è uguale a -1614.141a.u..
Come nel caso precedente abbiamo monitorato le distanze degli atomi dallo ione
metallico. Osservando la Fig.6.7 notiamo che l’amminoacido His-6 si allontana
immediatamente dallo ione metallico e si porta a distanze superiori ai 5Å. Ricor-
diamo (Tabella 6.4) che l’atomo Nδ (His-6) già dopo la fase di minimizzazione
classica oscillava attorno a distanze lievemente superiori a quella della sfera di co-
ordinazione dello Zn2+ . Vediamo, inoltre, che dopo ∼140 cicli una delle molecole
di solvente (Sol-661) esce fuori dalla sfera di coordinazione dello ione metalli-
co portandosi progressivamente ad una distanza di ∼4Å. Stesso destino per l’O
(Gly-9) che dopo ∼150 cicli si porta leggermente fuori dalla sfera di coordinazio-
Figura 6.7: Grafico delle distanze di vari atomi dallo Zn2+ durante la fase d’ottimizzazione
ONIOM
140
6.2 Risultati
Atomi ϑ (deg)
Figura 6.8: Struttura finale del sistema Figura 6.9: Ingrandimento del sito metalli-
Zn2+ -Aβ1−16 + H2 O co
141
Simulazioni QM/MM-ONIOM
Figura 6.10: RMSD del sistema Zn2+ -Aβ1−16 + H2 O rispetto alla struttura di partenza. (a)
Sistema ONIOM; (b) Regione QM; (c) Regione MM
142
6.3 Conclusioni
effettuati su modelli piccoli, che l’introduzione del solvente non modifichi notevol-
mente le cariche degli atomi coinvolti nel legame con lo ione metallico [ref].
Conludendo osserviamo che ottimizzando con ONIOM il nostro sistema pas-
sa da una struttura esa-coordinata ad una tetra-coordinata. Nella Tabella 6.9 ri-
portiamo le distanze finali degli atomi che si trovano all’interno della sfera di
coordinazione del metallo.
Atomi r(Å)
O2(Glu-11) 1.96
Nε(His-13) 2.05
Nδ (His-14) 2.03
O(Sol-774) 1.98
Tabella 6.9: Distanze degli atomi che si trovano nella sfera di coordinazione dello Zn2+ al termine
dell’ottimizzazione ONIOM
6.3 Conclusioni
Le simulazioni ONIOM descritte in questo capitolo mostrano delle notevoli
differenze nella coordinazione dello ione metallico in presenza o in assenza di
acqua. Si osserva infatti che in assenza di solvente tutte e tre le istidine presenti
nel peptide Aβ1−16 restano nella sfera di coordinazione del metallo mentre in pre-
senza del solvente una di loro viene sostituita da una molecola d’acqua. Quest’ul-
tima struttura è particolarmente interessante in quanto lo ione metallico potrebbe
legare un secondo peptide perdendo la molecola d’acqua. Ciò supporterebbe alcu-
ni risultati sperimentali discussi nei capitoli precedenti [96] nei quali veniva evi-
denziato come lo Zn fosse maggiormente disponibile di altri ioni metallici (Cu)
a legami “inter-molecolari” che potrebbero rappresentare l’inizio di processi di
aggregazione.
143
Simulazioni QM/MM-ONIOM
Ciò che ci proponiamo di fare nel futuro sono delle simulazioni sia classiche
sia quantistiche di un sistema composto da due peptidi entrambi legati al metal-
lo attraverso le due istidine (His-13 e His-14) per verificare la stabilità di tale
modello.
144
Conclusioni
146
6.3 Conclusioni
ordinato ed una delle tre istidine entra ed esce dalla sfera di coordinazione dello
ione metallico. Entrambe queste strutture sono state sottoposte ad ottimizzazione
strutturale ONIOM.
147
Simulazioni QM/MM-ONIOM
I nostri risultati, oltre a mostrare l’importanza dell’acqua nello studio dei sis-
temi biologici, costituiscono una base per intraprendere nuovi studi di dinamica
molecolare classica e calcoli quantomeccanici su sistemi in cui uno ione Zn2+
sia a ponte tra due peptidi. Per lo studio di questo modello ciò che riteniamo
interessante è confrontare i risultati che si ottengono con tecniche computazio-
nali quali la dinamica molecolare Car-Parrinello, per la quale è però necessario
ridurre il numero di atomi in studio, con tecniche ibride, quali ONIOM o una tec-
nica QM/MM in cui la parte QM sia dinamica (Car-Parrinello) invece che statica
(ottimizzazione ONIOM), con le quali è possibile studiare il sistema completo
immerso nel solvente.
148
Appendice A
• Orientamento temporale;
• Orientamento spaziale;
150
Il MMSE è una scala largamente diffusa in ambito clinico, è utilizzabile da
medici o da altro personale dopo breve addestramento ed è di rapido impiego; la
sua somministrazione richiede un tempo variabile da 5 a 15 minuti e tale brevità
lo rende meno impegnativo per le risorse attentive del soggetto, rispetto ad una
batteria completa di test neuropsicologici, cosı̀ da poter essere utilizzato anche
nelle fasi avanzate della demenza.
I punteggi ai test comunque, non permettono da soli di stabilire una diagnosi
di demenza nè di determinarne l’eziologia; pertanto il MMSE dovrebbe essere uti-
lizzato come strumento in grado di suggerire, in caso di punteggi bassi, il ricorso
a ulteriori approfondimenti. Inoltre, non consentendo una valutazione completa
delle funzioni cognitive, non è sufficientemente sensibile alle fasi iniziali della
demenza (specificità del 96% e sensibilità del 63%) ovvero tende a sottostimare
tali casi.
Il punteggio totale, dato dalla somma delle risposte corrette che il soggetto ha
ottenuto in ciascun item, può andare da un minimo di 0 (massimo deficit cognitivo)
ad un massimo di 30 (assenza di deficit cognitivo).
Il punteggio soglia ai fini della diagnosi di disturbi dell’efficienza intellettiva è
23 e la maggior parte delle persone anziane non dementi ottiene punteggi superiori
a tale soglia. In un ampio studio di revisione del MMSE sono stati proposti tre
cut-score:
151
Mini Mental State Examination
152
153
Appendice B
Su ognuna delle facce del cubo verranno a trovarsi in media N 2/3 atomi e in totale
la frazione di particelle che si trovano sulla superficie del cubo è pari a 6N 1/3 ,
ovvero nel caso N = 103 il 60% si trova sulla superficie, nel caso di N = 106 il
6%. Questo determina un non trascurabile effetto di superficie che viene curato
utilizzando una topologia toroidale, ovvero circondando il box di simulazione con
delle copie del box stesso. Questo artificio impone però di prestare attenzione
nella trattazione delle interazioni a lungo raggio quali quelle di Coulomb e di
Van der Waals, mediante le quali le particelle interagiranno con le proprie copie
contenute nei box confinanti. Per questo motivo, le interazioni di non legame
vengono calcolate solo per gli atomi contenuti entro un certo raggio, detto di cut-
off, il cui valore non deve superare la metà della dimensione minima del box. Al di
155
Condizioni periodiche al contorno
fuori del raggio di cut-off, il potenziale è calcolato effettuando uno “shift” che lo
porta a zero entro un secondo raggio (detto appunto di “shift”), oppure con metodi
di somme su reticolo, o non viene calcolato affatto.
156
Appendice C
Vincoli geometrici
definiamo la matrice
∂ σh
Bhl = (C.3)
∂ rl
se i vincoli sono olonomi:
dσ dr
= B̂ = 0 (C.4)
dt dt
d2σ d 2 r d B̂ dr
= B̂ + =0 (C.5)
dt 2 dt 2 dt dt
158
il sistema può essere scritto in maniera compatta:
d 2~r
− M̂ 2 − B̂T λ + ~f = 0 (C.6)
dt
moltiplicando a sinistra per B̂M̂ −1 si ottiene:
d B̂ d~r
+ B̂M̂ −1 B̂T λ + B̂M̂ −1 ~f = 0 (C.7)
dt dt
da cui si ricava:
d B̂ d~r
B̂T λ = −B̂T (B̂M̂ −1 B̂T )−1 B̂M̂ −1 ~f − B̂T (B̂M̂ −1 B̂T )−1 (C.8)
dt dt
che sostituito nel sistema lineare (C.6):
d 2~r d B̂ d~r
2
= (1 − T̂ B̂)M̂ −1 ~f − T̂ (C.9)
dt dt dt
T̂ = M̂ −1 B̂T (B̂M̂ −1 B̂T )−1 B̂ (C.10)
L’algoritmo noto come LINCS [100], applica le formule (C.9) e (C.10), ma per
evitare di dover calcolare la matrice inversa tra parentesi, usa un piccolo trucco.
Viene costruita la matrice diagonale Ŝ tale che:
s !−1
1 1
S i j = δi j + (C.11)
mi m j
in questo modo
= S(SBM −1 BT S)−1 S
= Ŝ(1̂ − Â)−1 Ŝ
Nel caso in cui i constraint siano solo sulle distanze di legame, la matrice  è
simmetrica, e ha tutti gli autovalori minori di 1, quindi è possibile espandere:
∞
(1̂ − Â)−1 = ∑ Ai (C.13)
i=0
Nelle implementazioni pratiche vengono usati i primi 4-8 termini della somma.
159
Appendice D
Basis set
161
Basis set
esatti per l’atomo d’idrogeno. Tuttavia, le STO non posseggono tutti i nodi radiali;
i nodi nella parte radiale sono introdotti attraverso una combinazione lineare delle
STO stesse. La dipendenza esponenziale assicura una rapida convergenza all’au-
mentare del numero di funzioni, tuttavia i calcoli degli integrali a due elettroni
non possono essere svolti analiticamente.
Gli orbitali di tipo Gaussian [102] possono essere scritti sia in termini di
coordinate polari che in termini di coordinate cartesiane, cioè:
2
χζ ,n,l,m (r, θ , φ ) = NYl,m (θ , φ )r2n−2−l e−ζ r (D.2)
2
χζ ,n,l,m (x, y, z) = Nxlx yly zlz e−ζ r (D.3)
162
dal nucleo rispetto alle STO, e la “coda” della funzione d’onda è rappresentata in
modo scarso. Sia le GTO che le STO possono essere scelte per formare un set
completo di basi, ma per quello che abbiamo appena detto le GTO necessitano
di più funzioni di base per ottenere la stessa accuratezza ottenuta dalle STO con
un numero inferiore di basi. Delle linee guida generali indicano che ci vogliono
circa il triplo di basi GTO per ottenere il medesimo livello di accuratezza ottenuto
con le STO. L’incremento nel numero delle funzioni di base GTO, tuttavia, è più
che compensato dalla facilità con la quale si riesce a calcolare gli integrali. In
termini di efficienza computazionale, le GTO sono quindi preferite e sono usate
quasi universalmente come funzioni di base nei calcoli di struttura elettroniche.
In realtà come componenti del basis set non si prendono direttamente le fun-
zioni (D.2-D.3), bensı̀ delle combinazioni lineari di esse, dette “contrazioni”:
K
χcon = ∑ ai χi (D.4)
i=1
Esistono vari schemi di contrazione. Nella contrazione segmentata ogni primitiva
χi viene utilizzata nell’espansione di una sola funzione di base χcon ; mentre nella
contrazione generale non vi sono limiti all’utilizzo delle primitive. Questo schema
permette di modulare l’accuratezza con cui i singoli orbitali vengono costruiti.
Gli elettroni di core, in genere hanno un maggior peso energetico, per questo
vengono usate delle contrazioni maggiori per riprodurne gli orbitali molecolari,
usare lo stesso numero per gli elettroni più esterni appesantirebbe il calcolo non
apportando migliorie apprezzabili. Tuttavia gli elettroni di valenza sono quelli che
determinano la chimica della molecola in esame, per questo in alcuni casi possono
essere aggiunte delle funzioni diffuse (con esponente ζ piccolo), o dei termini di
momento angolare più alto detti termini di polarizzazione.
Piuttosto che partire con funzioni di base mirate a modellare gli orbitali ato-
mici (STO oppure GTO), ed usare una combinazione lineare di queste per de-
scrivere gli orbitali per l’intero sistema, si possono usare funzioni che mirano
163
Basis set
~
χk (~r) = eik·~r (D.8)
164
solo dalla grandezza della cella periodica, non dal sistema che si trova realmente
all’interno della cella. Mentre nel caso delle Gaussiane centrate sul nucleo queste
crescono linearmente con le dimensioni del sistema, cioè le onde piane sono più
efficienti per sistemi di grandi dimensioni.
Inizialmente le onde piane furono usate per descrivere esclusivamente siste-
mi periodici, ma negli ultimi anni sono state usate per descrivere specie moleco-
lari senza alcuna periodicità, utilizzando l’approccio a supercella [103], dove la
molecola è posta in una cella unitaria sufficientemente grande perchè non inte-
ragisca con la sua immagine. Le onde piane sono ideali per descrivere densità
elettroniche delocalizzate lentamente variabili, come le bande di valenza
165
Appendice E
Il principio di Ritz
H | φn i = En | φn i (E.1)
Dimostrazione
| ψi = ∑ cn | φ i (E.3)
n
167
Il principio di Ritz
e quindi:
2
hψ | H|ψi − hφ0 | H|φ0 i = ∑ |cn| En − E0 = ∑ |cn |2 (En − E0 ) (E.6)
n n
Nell’ultimo passaggio si è fatto uso della (E.4). Nella (E.6) è chiaro che tutti
i termini della sommatoria con n 6= 0 sono positivi (infatti En > E0 ), mentre il
termine con n = 0 è evidentemente nullo. Di conseguenza l’espressione (E.6) avrà
il suo minimo assoluto quando tutti i cn con n 6= 0 sono nulli, in tal caso c0 = 1,
ovvero |φ i = |φ0 i e l’espressione (E.6) vale zero. Pertanto risulta dimostrato che
la (E.6) è maggiore di zero per tutti gli |φ i 6= |φ0 i
168
Appendice F
Algoritmi di ottimizzazione
L’ottimizzazione è una branca della matematica che studia teoria e metodi per
la ricerca dei punti stazionari di un modello che traduce in termini matematici
un dato problema (non occupandosi quindi direttamente di come tale modello sia
stato costruito). L’ambito di ricerca privilegiato dell’ottimizzazione è quello dei
modelli esprimibili in termini di funzioni a più variabili, nei quali i punti stazionari
vengono ricercati ponendo anche vincoli qualitativi espressi in termini di derivate
successive.
Un qualsiasi problema di ottimizzazione può essere espresso nella seguente
forma:
min f (x)
x
(F.1)
~x ∈ X ⊆ ℜ
170
Generalmente la ricerca è indirizzata sui minimi o sui punti di sella, sia globali
che locali. Potremmo pensare che per calcolare il minimo di una funzione a più
variabili, basti derivare in successione rispetto ad ognuna delle variabili lasciando
le altre fisse, e cosı̀ facendo raggiungere il minimo dopo un certo numero di passi.
Quest’approccio però suppone che le variabili siano indipendenti fra di loro, ed
inoltre diviene impraticabile per funzioni con un numero di variabili relativamente
piccole, da cinque a dieci.
Il problema può essere schematizzato in questo modo. Supponiamo di avere
una funzione f (~x), dove x è una variabile sconosciuta, e di voler trovare un mini-
mo di f . Naturalmente la prima cosa da fare sarà scegliere il punto di partenza ~x
da cui iniziare la ricerca, che ipotizziamo essere un certo x~0 in cui la f è nota. Na-
turalmente la scelta di tale punto può influenzare la ricerca del minimo, per questo
motivo si cerca sempre questo punto sulla base della conoscenza dell’andamento
di f (~x), in maniera tale da trovare un valore che non sia troppo distante dal mi-
nimo ricercato. A questo punto dobbiamo decidere due cose: (i) la direzione da
seguire nella ricerca del minimo, cioè dobbiamo scegliere la direzione da seguire
nel cammino dal punto x~0 a quello successivo; (ii) quanto grande deve essere
lo spostamento lungo quella direzione. Avremmo la seguente rappresentazione
iterativa:
dove îk rappresenta la direzione e |λk îk | la grandezza dello spostamento (step).
La differenza tra i vari metodi di ottimizzazione è determinata dalla scelta della
direzione e dello step. Possiamo classificare i metodi approssimativamente in
questo modo:
1. Metodi che usano le derivate prime della funzione nel punto (metodi a
gradiente).
171
Algoritmi di ottimizzazione
172
dove T sta per trasposto. Ponendo quest’espressione uguale a zero, vediamo quali
direzioni scegliere via via e notiamo che sono ortogonali tra loro. Il prossimo pas-
so è allora spostarci nel punto che si trova nella direzione negativa del gradiente e
dopo n step otterremo un cammino a zig-zag (Fig.F.1). Quest’iterazione continua
fino ad un estremo che viene determinato in base all’accuratezza da noi scelta.
Questo è un problema di minimizzazione lungo una linea, dove la linea è data
dalla (F.4) per differenti valori di λk . Quindi la ricerca del minimo si riduce ad
una sequenza di ricerche lineari.
Figura F.1: Illustrazione bidimensionale del metodo SD. Si vede proprio come l’approccio al
minimo avviene attraverso una traiettoria a zig-zag.
173
Algoritmi di ottimizzazione
d~iT · Ŵ · d~ j = 0 (F.6)
174
L’ idea è quella di prendere la direzione ~di dipendente da tutte le altre direzioni
possibili per localizzare il minimo di f attraverso l’equazione (F.6). Un insieme
di tali direzioni di ricerca è detto Ŵ -ortogonale o coniugato. Il modo migliore di
visualizzare il metodo CG è il seguente: supponiamo di avere lo spazio in cui sti-
amo lavorando e uno spazio “allungato” (Fig.F.3). La Fig.F.3a mostra la forma del
contorno di una funzione quadratica nello spazio reale, che è ellittica per ~b 6= 0.
Figura F.3: Schema d’ottimizzazione del CG (a). Tutte le direzioni presenti sono ortogonali. (b)
Lo stesso problema visto nello spazio “allungato”, dove le linee sono Ŵ -ortogonali.
Qualsiasi coppia di vettori che appare perpendicolare in questo spazio deve essere
ortogonale. La Fig.F.3b mostra lo stesso contorno in uno spazio che è allunga-
to lungo gli assi delimitati dagli autovettori della matrice Ŵ in modo tale che il
contorno ellittico diventi circolare. Una qualsiasi coppia di vettori che è perpen-
dicolare in questo spazio è, infatti, Ŵ -ortogonale. La ricerca del minimo nel caso
di una funzione quadratica parte dal punto x~0 (Fig.F.3a) ed effettua uno step lungo
la direzione d~0 e si ferma nel punto x~1 . Questo è un punto di minimo lungo quella
direzione, determinato con lo stesso metodo dello SD, cioè il minimo che si trova
nella direzione dove la derivata direzionale è nulla (F.5). La differenza essenziale
175
Algoritmi di ottimizzazione
tra il CG e lo SD è legata alla scelta della direzione da seguire dopo il primo pas-
so. Mentre il metodo SD dovrebbe seguire la direzione delimitata dal vettore ~r1
(Fig.F.3a), il metodo CG dovrebbe scegliere come direzione quella delimitata dal
vettore d~1 . Come riesce il CG a trovare la giusta direzione che porta al minimo ~x?
La risposta si trova nella Fig.F.3b; nello spazio allungato, la direzione d~0 sembra
essere tangente al contorno circolare nel punto x~1 . Poichè la direzione successiva
d~1 è costretta ad essere Ŵ -ortogonale alla precedente, essa apparirà perpendico-
lare in questo spazio modificato. Quindi, d~1 ci porterà direttamente al minimo
della funzione quadratica f (~x). Per evitare di ricercare lungo direzioni già esplo-
rate, il CG garantisce che la minimizzazione di f lungo una direzione non “rovini”
la minimizzazione lungo le altre; cioè che dopo i step, la f sarà minimizzata lungo
tutte le direzioni esplorate.
La F.6 afferma che una ricerca lungo d~i mostra dove il gradiente di f è ortogo-
nale a d~i e come muoverci adesso lungo la nuova direzione d~ j . Il gradiente subisce
una variazione:
δ (∇ f ) = Ŵ · d~ j (F.7)
Per non interferire con la minimizzazione lungo d~i , si richiede che il gradiente
resti perpendicolare a d~i ; cioè che il cambiamento nel gradiente sia esso stesso
ortogonale a d~i , come si vede dalla F.6. Per una funzione quadratica, come nella
F.3, la procedura è la seguente. Lo step iniziale è nella direzione dello SD:
176
La lunghezza dello step lungo ogni direzione è data da:
d~kT ·~gk
λk = (F.11)
d~T · (Ŵ · d~k )
k
Il metodo NR differisce dai metodi appena visti nel fatto che il minimo viene
localizzato utilizzando l’informazione relativa alle derivate seconde della funzione
177
Algoritmi di ottimizzazione
f (x). Questo porta ad una più rapida convergenza, ma non ad un tempo compu-
tazionale necessariamente inferiore. Infatti, il calcolo delle derivate seconde e
la gestione della matrice Ĥ può richiedere uno sforzo computazionale enorme,
specialmente per sistemi molto grandi.
L’idea del metodo NR è approssimare la f (~x) in ogni iterazione attraverso una
funzione quadratica (F.3) e muovere questa funzione verso il minimo. La funzione
quadratica su un punto ~x in un adatto intorno del punto ~xk è data da una serie di
Taylor troncata:
1
f (~x) ≈ f (~xk ) + (~x −~xk )T ·~gk + (~x −~xk )T · Ĥk · (~x −~xk ) (F.13)
2
dove sia il gradiente g~k che la matrice Ĥk sono valutate in ~xk . Le derivate della
F.13 sono:
1 1
∇ f (~x) = ~gk + Ĥk · (~x −~xk ) + ĤkT · (~x −~xk ) (F.14)
2 2
La matrice Hessiana Ĥ è sempre simmetrica se la funzione ~xk è continua, dif-
ferenziabile due volte e con derivate continue su ogni punto. Quindi quest’ultima
equazione si riduce a:
∇ f (~x) = ~gk + Ĥk · (~x −~xk ) (F.15)
che è un sistema lineare di equazioni. Il metodo NR usa il punto ~x∗ come punto
successivo presente nella formula iterativa:
178
Benchè la convergenza possa sembrare molto rapida, osservando il numero di
step, ogni iterazione include il calcolo delle derivate seconde e la gestione della
matrice Ĥ. Le performance del metodo sono quindi dipendenti da alcune qualità di
questa matrice. Una di queste qualità è che la matrice deve essere definita positiva.
Infatti, affinchè il metodo converga verso il minimo, la direzione di Newton deve
essere una direzione in discesa. Perciò, si richiede che:
179
Algoritmi di ottimizzazione
180
Appendice G
Pseudopotenziali
Gli elementi chimici che si trovano nella parte più bassa della tavola periodica
hanno un elevato numero di elettroni di core. Questi elettroni non sono impor-
tanti dal punto di vista dello studio del legame chimico, poichè non partecipano
direttamente alla sua formazione o rottura ed inoltre poiché è necessario usare un
gran numero di funzioni di base per espandere i corrispondenti orbitali si potrebbe
pensare ad una via alternativa per la loro rappresentazione attraverso le funzioni di
base. Bisogna però tenere a mente che gli orbitali di valenza per essere descritti in
maniera adeguata necessitano di una buona rappresentazione degli orbitali di core,
a causa degli effetti di repulsione elettrone-elettrone. Inoltre negli elementi ap-
partenenti alla parte più bassa della tavola periodica gli effetti relativistici giocano
un ruolo essenziale complicando notevolmente i calcoli. Questi problemi pos-
sono essere “risolti” simultaneamente modellando gli elettroni di core attraverso
un’adatta funzione, e trattando solamente gli elettroni di valenza esplicitamente.
182
si (cristalli), e questa differenza è rispecchiata anche nei corrispondenti PP. Quan-
do usiamo funzioni Gaussiane per descrivere gli orbitali di valenza, è ovvio usare
le stesse funzioni per descrivere lo PP. Poichè le funzioni Gaussiane sono conti-
nue, non c’è una distanza fissa che caratterizza l’estensione del potenziale di core
e la qualità dell’ECP è determinata dal numero di elettroni che attraverso esse si
voglione rappresentare. Per i metalli di transizione, è chiaro che gli orbitali più
esterni, cioè (n + 1)s, (n + 1)p ed (n)d, costituiscono lo spazio di valenza. Se
scegliamo di modellizzare il resto degli orbitali attraverso uno PP, otterremo dalle
geometrie ragionevoli, ma troveremo che in alcuni casi le energie del sistema non
sono soddisfacenti. Il risultato migliore si ottiene includendo nel modello anche
gli orbitali che sono prossimi a quelli di valenza. In questo modo aumentano i
costi computazionali, ma si è certi che le energie derivate sono soddisfacenti. Per
esempio, nel caso dell’argento che ha numero atomico 47, possiamo considera-
re due differenti scelte per gli elettroni di core, oltre quella di cosiderarli tutti
esplicitamente, da simulare con un PP:
Il guadagno nell’usare l’ECP è maggiore per atomi che appartengono alla parte
bassa della tavola periodica, specialmente per quelli dove gli effetti relativistici
sono importanti.
Il numero di onde piane da usare per rappresentare i nostri orbitali dipende
direttamente dal vettore d’onda, associato all’energia più grande, che è inversa-
183
Pseudopotenziali
mente proporzionale alla più piccola variazione della funzione d’onda che può
essere descritta. La singolarità del potenziale nucleare (Vne ) e i corrispettivi elet-
troni di core fortemente localizzati sono essenzialmente impossibili da descrivere
attraverso un numero ragionevole di onde piane. Gli PP sono quindi utilizzati per
“spalmare” la carica nucleare e modellare gli elettroni di core. Questi potenziali
sono caratterizzati tipicamente da un “raggio di core” rc . Cioè gli PP usati in
connessione con le onde piane hanno un estensione fisica finita. Il potenziale per
distanze più piccole di rc è descritto da un’adatta funzione analitica, tipicamente
un polinomio oppure una funzione sferica di Bessel. E’ chiaro che uno PP “duro”
(piccolo rc ) richiederà un numero maggiore di onde piane per descrivere la re-
gione al di sotto di rc rispetto ad uno “soffice” (grande rc ), ma quando rc diviene
troppo grande la qualità dei risultati calcolati si deteriora ed anche la trasferibilità
dello PP si perde.
Mentre non genera nessun errore l’uso delle onde piane come funzioni di base
per espandere gli orbitali, l’uso dello PP per descrivere la regione di core comporta
un limite fondamentale all’accuratezza con cui possiamo svolgere i nostri calcoli.
Per sistemi composti da elementi delle prime due righe della tavola periodica, l’er-
rore nei calcoli DFT è approssimativamente equivalente a quello che si otterrebbe
usando basi Gaussiane per descrivere tutti gli elettroni [109]. Una limitazione im-
184
plicita al metodo degli PP è l’incapacità di descrivere le proprietà molecolari che
dipendono direttamente dagli elettroni di core (come nella X-ray photoelectron
spectroscopy (XPS)) o dalla densità elettronica vicino ai nuclei (come nell’NMR)
[51].
185
Appendice H
Algoritmo di Berny
187
Algoritmo di Berny
Gram-Schmidt.
a−1
r̃ia = (xia − xi0 ) − ∑ rib ∑(xia − xi0)rbj (H.1)
b=1 j
a−1
∑ (gai − g0i )rib − ∑ kcb ∑(xia − xi0)ric
i c=1 i
knab = a 0 a
6= knba
∑ (xi − xi )ri
i
vamene in termini del decremento dei contributi non quadratici. Quindi, per b > a,
188
ci aspettiamo che knab abbia una componente non quadratica più piccola rispetto a
knba .
Se l’ampiezza del vettore gradiente è molto grande, possiamo migliorare la
situazione evitando di calcolare la matrice delle derivate seconde fino a che i gra-
dienti non assumono un valore al di sotto di una certo cutoff. Questo comporta
che negli ultimi step dell’ottimizzazione non spenderemo tempo nel correggere le
derivate seconde ottenute nei primi step. L’effetto netto di questo approccio è uti-
lizzare nei primi step un algoritmo SD, fino a che i gradienti non scendono sotto
un valore ragionevole.
Punto b. Bisogna ricercare un minimo per E nella direzione ~xa −~xb attraver-
L’inverso della matrice delle derivate seconde è calcolato attraverso una diago-
nalizzazione della matrice. Gli autovalori possono allora essere testati in maniera
tale da assicurare che la matrice sia definita positiva e perciò che lo step d’ottimiz-
zazione proceda realmente verso un minimo. Se troviamo un autovalore negativo,
invertiamo il suo segno. L’effetto netto è quello di forzare una procedura di SD
lungo la direzione degli autovettori associati agli autovalori negativi. La matrice
189
Algoritmo di Berny
inversa può allora essere costruita attraverso gli autovalori e gli autovettori. Inoltre
si richiede che la variazione totale su un qualsiasi xi deve essere al dı̀ sotto di una
certa soglia. Se il cambiamento in ~c è troppo grande, allora la variazione prevista
è scalata in modo tale che i limiti non siano eccedenti.
Lo step successivo è simile a quelli già visti, in cui si utilizzano le informazioni
che vengono dagli n step precedenti e si eliminano le informazioni derivanti dagli
step più “vecchi”. Ad ogni step il gradiente e i vettori spostamento sono testati
per vedere se l’ottimizzazione è giunta a convergenza. Lo scarto quadratico medio
(RMS) del gradiente e del valore assoluto della più grande componente del gra-
diente devono essere al di sotto del loro rispettivo valore di soglia. Lo spostamento
spaziale stimato deve superare un test simile sull’RMS e sul valore assoluto della
componente più grande. Se queste quattro condizioni sono soddisfatte, l’ottimiz-
zazione è completata. Se le condizioni non sono soddisfatte, il processo continua
a meno che non superi il numero massimo di step imposti [110].
190
Appendice I
192
legami molecolari. L’idea del Natural Atomic Orbital (NAO) e dell’analisi NBO
[111] è quello di usare la matrice densità one-electron per definire la forma degli
orbitali atomici nell’ambiente molecolare e derivare i legami molecolari tra gli
atomi attraverso la densità elettronica.
Assumiamo che le funzioni di base siano distribuite in maniera tale che tutti
gli orbitali localizzati sul centro A vengono prima di quelli localizzati sul centro
B, che sono prima di quelli localizzati sul centro C, ecc.
193
Natural Bond Orbital (NBO)
NAO dipende dalla grandezza del set di base atomico, in particolare, il numero
di NAO-Rydberg incrementa all’aumentare del set di basi. Quindi è desiderabile
che la procedura d’ortogonalizzazione preservi la forma degli orbitali fortemente
occupati quanto più possibile. La procedura è la seguente (Fig.I.1):
4. Ogni NAO debolmente occupato è reso ortogonale a tutti gli altri NAO
debolmente occupati associati a centri diversi attraverso una procedura di
diagonalizzazione pesata attraverso il numero d’occupazione.
194
Il set finale di orbitali ortogonali sono semplicemente denotati NAO, e gli
elementi diagonali della matrice densità in queste basi sono le popolazioni degli
orbitali. Sommando tutti i contributi degli orbitali appartenenti ad uno specifico
centro si ottiene la carica atomica. Generalmente si trova che i “natural mini-
mal” NAO contribuiscono a più del 99% della densità elettronica, e formano una
rappresentazione molto compatta della funzione d’onda in termini degli orbitali
atomici. L’ulteriore vantaggio dei NAO è che sono definiti dalla matrice densità,
garantendo in questo modo che il numero d’occupazione associato agli elettroni
sia compreso tra 0 e 2, e che convergono ad un ben definito valore all’aumentare
del basis set. Inoltre, l’analisi può essere effettuata anche per funzioni d’onda
correlate. Lo svantaggio delle NAO è che la loro estensione può superare di gran
lunga le dimensioni dell’atomo a cui sono associate e, quindi, possono descrivere
densità elettroniche che si trovano vicine ad un nucleo differente da quello per cui
sono state create e su cui sono centrate.
Quindi la matrice densità è stata trasformata nelle basi NAO ed i legami tra
gli atomi possono essere identificati attraverso i blocchi che si trovano fuori dalla
diagonale. La procedura coinvolge i seguenti step:
1. I NAO per un blocco atomico nella matrice densità che hanno numero d’oc-
cupazione molto vicino a 2 (diciamo >1.999) sono identificati come orbitali
di core. Il loro contributo dalla matrice densità è rimosso.
2. I NAO per un blocco atomico nella matrice densità che hanno un grande
numero d’occupazione (diciamo >1.90) sono identificati come orbitali con
un doppietto elettronico. Il loro contributo dalla matrice densità è rimosso.
3. Ogni coppia di atomi (AB, AC, BC, . . .) viene adesso considerata, ed i sot-
toblocchi a coppie di due della matrice densità sono diagonalizzati (senza i
contributi di core e dei doppietti elettronici). Gli NBO sono identificati co-
195
Natural Bond Orbital (NBO)
Una volta ricavati gli NBO questi possono essere scritti come una combinazione
lineare dei NAO, generando un modello localizzato degli “orbitali” atomici che
sono coinvolti nei legami.
196
Ringraziamenti
Il mio primo pensiero va alla mia famiglia. Non esiste un modo per ripagare
appieno i sacrifici che avete fatto per me. A mia madre che non si è mai stancata
e con forza mi ha sostenuto nei momenti bui. A mio padre che mi fa sentire or-
goglioso di essere suo figlio quando vedo quello che ha fatto e come lo ha fatto. A
mio fratello perché è il mio “miracolo”, un mio punto di riferimento e mi è sempre
stato accanto anche da lontano. A mia sorella che riesce a sorprendermi sempre
e che per troppo tempo non ha avuto un “buon” fratello ma che col tempo, spero,
riuscirà ad apprezzare. Alla mia famiglia allargata Rossana e Fabrizio. Infine ai
miei due nipotini che riescono a farmi piangere come un bambino. Grazie.
Un grazie a Cristina per avermi sopportato, per aver curato tutte le mie ferite
198
e per avermi fatto ridere nei momenti peggiori.
Vorrei ringraziare con sincero affetto Filippo per tutto quello che ha fatto, per
l’aiuto nella stesura della tesi e per avermi dato forza in un momento di sconforto.
Un grazie a tutti i miei amici. Non me ne vogliano se non li cito uno per uno, ma
“sta tesi adda fernı̀”.
Infine il più importante dei ringraziamenti lo riservo per Dio che mi ha da-
to la forza di superare dei momenti molto bui, mi ha protetto e non mi ha mai
abbandonato lungo questo cammino. Grazie Signore.
199
Bibliografia
[2] G.A. Petsko, D. Ringe. Protein Structure and Function. New Science Press
Ltd, London, 2004.
[3] R.B. Martin. Peptide bond characteristics. Met. Ions Biol. Syst., 38: 1-23,
2001.
[4] D. Voet, J.G. Voet. Biochemistry, 2nd edition. Wiley, New York, 1995.
[6] L. Pauling, R.B. Corey, and H.R. Branson. The structure of proteins: two
hydrogen-bonded helical configurations of the polypeptide chain. PNAS,
37: 235-240, 1951.
[7] S.H. Gellman. Minimal model systems for beta sheet secondary structure
in proteins. Curr. Opin. Chem. Biol., 37: 729-740, 1998.
200
BIBLIOGRAFIA
[11] C.B. Anfinsen. Principles that govern the folding of protein chains. Science,
181: 223230, 1973.
[13] C.M. Dobson. Protein misfolding, evolution and disease. Trends Biochem.
Sci., 24: 329-332, 1999.
[14] D.J. Selkoe. Folding protein in fatal ways. Nature, 426: 900-904, 2005.
[16] V.N. Uversky , A. Talapatra, J.R. Gillespie and A.L. Fink. Protein de-
posits as the molecular basis of amyloidosis. Part I. Systemic amyloidoses.
Medical Science Monitor, 5: 1001-1012, 1999.
[17] A.I. Bush. Metals and Neuroscience. Curr. Opin. Chem. Biol., 4:184-191,
2000.
201
BIBLIOGRAFIA
202
BIBLIOGRAFIA
[29] M.R. Meyer, J.T. Tschanz, M.C. Norton, K.A. Welsh-Bohmer, D.C.
Steffens, B.W. Wyse, J.C. Breitner. APOE genotype predicts when-not
whether-one is predisposed to develop Alzheimer disease. Nature Genetics,
19: 321-322, 1998.
[31] D.M. Mann, P.O. Yates, B. Marcyniuk, C.R. Ravindra. The topography
of plaques and tangles in Down’s syndrome patients of different ages.
Neuropathology & Applied Neurobiology, 12: 447-457, 1986.
203
BIBLIOGRAFIA
[36] J. Kyte, R.F. Doolittle. A simple method for displaying the hydrophobic
character of a protein. J. Mol. Biol., 157: 105-132, 1982.
[39] V.N. Uversky, A.L. Fink. Conformational constraints for amyloid fibrilla-
tion: the importance or being unfolded. Biochimica et Biophysica Acta,
1698: 131-153, 2004.
[40] A.I. Bush, W.H. Pettingell, G. Multhaup, M. Paradis, J.P. Vonsattel, J.F.
Gusella, K. Beyreuther, C.L. Masters, R.E. Tanzi. Rapid induction of
Alzheimer Aβ -amyloid formation by zinc. Science, 265: 1464-1467, 2003.
204
BIBLIOGRAFIA
[49] W.F. van Gunsteren, P.K. Weiner, A.P. Wilkinson. Computer simulation of
biomolecular systems, theoretical and experimental applications, Vol. 3.
Escom Science Publishers, Leiden, 1993.
[50] G.C. Rossi. Appunti di Teoria dei Sistemi a Molti Corpi. Università degli
studi di Roma “Tor Vergata”, Roma, 2002.
205
BIBLIOGRAFIA
[54] K. Huang. Statistical mechanics, 2nd edition. Wiley, New York, 1987.
[55] L.V. Woodcock. Isothermal molecular dynamics calculation for liquid salt.
Chemical Physical Letters, 10: 257-261, 1971.
[56] H.J.C. Berendsen, J.P. Postma, W.F. van Gunsteren, A. Dinola, J.R. Haak.
Molecular dynamics with coupling to an external bath. Journal of Chemical
Physics, 81: 3684-3690, 1984.
206
BIBLIOGRAFIA
[61] J.M. Ziman. Principles of the Theory of Solids, 2nd edition. Cambridge
University Press, Cambridge, 1979.
[62] R.P. Feynmann, A.R. Hibbs. Quantum Mechanics and Path Integrals.
McGraw Hill, New Jersey, 1965.
[66] E. Bruno. Appunti del corso di fisica dello stato solido. Università degli
studi di Messina, Messina, 1996.
207
BIBLIOGRAFIA
[75] A. Warshel, R.M. Weiss. An empirical valence bond approach for com-
paring reactions in solutions and in enzymes. Journal of the American
Chemical Society, 102: 6218-6226, 1980.
[76] S.P. Webb, M.S. Gordon. Solvation of the Menshutkin Reaction: A Rigor-
ous Test of the Effective Fragment Method. Journal of Physical Chemistry
A, 103: 1265-1273, 1999.
208
BIBLIOGRAFIA
[84] D. Van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A.E. Mark, H.J.C.
Berendsen. GROMACS: Fast, Flexible and Free. Journal of Computational
Chemistry, 26: 1701-1718, 2005.
209
BIBLIOGRAFIA
[87] N. Foloppe; A.D. MacKerell Jr. All-Atom Empirical Force Field for Nu-
cleic Acids: I. Parameter Optimization Based on Small Molecule and
Condensed Phase Macromolecular Target Data. Journal of Computational
Chemistry, 21: 86-104, 2000.
[91] W.L. Jorgensen; J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein.
Comparison of simple potential functions for simulating liquid water.
Journal of Chemical Physics, 79: 926-935, 1983.
[93] M.J. Frisch, G.W. Trucks et al. Gaussian 03, Revision C.02. Gaussian, Inc.,
Wallingford, 2004.
210
BIBLIOGRAFIA
[95] A.K. Rappe, C.J. Casewit, K.S. Colwell, W.A. Goddard III, W.M. Skiff.
UFF, a Full Periodic Table Force Field for Molecular Mechanics and
Molecular Dynamics Simulations. Journal of the American Chemical
Society, 114: 10024-10035, 1992.
[97] C.D. Syme, R.C. Nadal, S.E.J. Rigby, J.H. Viles. Copper Binding to the
Amyloid-β (Aβ ) Peptide Associated with Alzheimer’s Disease: Fold-
ing, Coordination Geometry, pH Dependece, Stoichiometry, and Affini-
ty of Aβ -(1–28): Insights from a range of complementary spectroscopic
Techniques. J. Biol. Chem., 279: 18169-18177, 2004.
[98] M.F. Folstein, S.E. Folstein, P.R. McHugh. “Mini-mental state”. A practical
method for grading the cognitive state of patients for the clinician. Journal
of psychiatric research, 12(3): 189-198, 1975.
211
BIBLIOGRAFIA
[101] J.C. Slater. Atomic Shielding Constants. Physical Review, 36: 57-64, 1930.
[102] S.F. Boys. Electronic Wave Functions. II. A Calculation for the Ground
State of the Beryllium Atom. Proc. R. Soc. Lond. A, 201: 125-137, 1950.
[105] T.R. Cundari, M.T. Benson, M.L. Lutz, O.S. Sommerer. Effective Core Po-
tential Approaches to the Chemistry of the Heavier Elements. Reviews in
computational chemistry, 8: 145-202, 1996.
[106] M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias, J.D. Joannopoulos. Itera-
tive minimization techniques for ab initio total-energy calculations: molec-
ular dynamics and conjugate gradients. Rev. Mod. Phys., 64: 1045-1097,
1992.
212
BIBLIOGRAFIA
213