Universita degli Studi di Bologna

Facolta di Scienze Matematiche, Fisiche e Naturali

Corso di laurea in Fisica

MODELLI DINAMICI PER ILRIPIEGAMENTO DELLA

PROTEINA

Tesi di Laurea di: Lorenzo BERNACCHIONI

Relatore: Prof. Armando BAZZANI

Co-Relatore: Prof. Sandro RAMBALDI

Parole Chiaveprotein folding, dinamica molecolare, catene di corpi rigidi

algoritmo di integrazione, quaternioni

I Sessione

Anno Accademico 2003/2004

Indice

Introduzione 7

1 Proteine e ripiegamento 11

1.1 Il problema del Protein Folding . . . . . . . . . . . . . 11

1.2 La proteina . . . . . . . . . . . . . . . . . . . . . . . . 13

1.2.1 Amminoacidi . . . . . . . . . . . . . . . . . . . 13

1.2.2 Struttura secondaria . . . . . . . . . . . . . . . 17

1.3 Interazioni fisiche . . . . . . . . . . . . . . . . . . . . . 17

1.4 Un primo sguardo al modello . . . . . . . . . . . . . . 18

2 Quaternioni 21

2.1 Algebra dei quaternioni . . . . . . . . . . . . . . . . . . 23

2.1.1 Osservazioni . . . . . . . . . . . . . . . . . . . . 26

2.2 Applicazioni dei quaternioni . . . . . . . . . . . . . . . 28

2.3 Dinamica Molecolare e Quaternioni . . . . . . . . . . . 29

2.3.1 Riferimento per un corpo rigido . . . . . . . . . 30

2.4 I quaternioni nel modello per la proteina . . . . . . . . 31

3 Dinamica Molecolare e algoritmi di integrazione 33

3.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . 33

3.2 Generalita sugli integratori . . . . . . . . . . . . . . . . 33

3.3 Metodologie di integrazione . . . . . . . . . . . . . . . 34

3.3.1 Approccio atomico . . . . . . . . . . . . . . . . 34

3.3.2 Approccio molecolare . . . . . . . . . . . . . . . 35

3.4 Dinamica e integratori per una catena . . . . . . . . . 38

3

4 Integratore per un pendolo sferico 39

4.1 Dinamica . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.1.1 Pendolo sferico . . . . . . . . . . . . . . . . . . 39

4.1.2 Corpo rigido . . . . . . . . . . . . . . . . . . . . 41

4.2 Algoritmo per il pendolo sferico . . . . . . . . . . . . . 42

4.2.1 Integratore al 2o Ordine . . . . . . . . . . . . . 43

4.2.2 Integratore al 3o Ordine . . . . . . . . . . . . . 44

4.3 Catena di pendoli sferici . . . . . . . . . . . . . . . . . 45

4.4 Test al calcolatore . . . . . . . . . . . . . . . . . . . . . 47

4.4.1 Risultati . . . . . . . . . . . . . . . . . . . . . . 47

5 Modello per la proteina 51

5.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . 51

5.2 Potenziali . . . . . . . . . . . . . . . . . . . . . . . . . 52

5.3 Bagno termico . . . . . . . . . . . . . . . . . . . . . . . 55

5.3.1 Dinamica di Langevin . . . . . . . . . . . . . . 56

5.3.2 Dinamica di Langevin all’interno dell’algoritmo 58

5.3.3 Test termodinamico . . . . . . . . . . . . . . . . 58

5.4 Similarita χ . . . . . . . . . . . . . . . . . . . . . . . . 60

6 Simulazioni e risultati 63

6.1 Ipotesi teoriche sul protein folding . . . . . . . . . . . . 63

6.2 Simulazione . . . . . . . . . . . . . . . . . . . . . . . . 64

6.2.1 Sequenza . . . . . . . . . . . . . . . . . . . . . . 65

6.2.2 Dimensionalita . . . . . . . . . . . . . . . . . . 68

6.2.3 Modalita della simulazione . . . . . . . . . . . . 68

6.3 Analisi . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6.3.1 Classificazione dei minimi . . . . . . . . . . . . 69

6.3.2 Grafici Energia-Popolazione . . . . . . . . . . . 71

6.3.3 Grafici Energia-Similarita-Popolazione . . . . . 71

6.3.4 Sequenza idrofoba . . . . . . . . . . . . . . . . 74

6.4 Risultati . . . . . . . . . . . . . . . . . . . . . . . . . . 75

Conclusioni 77

4

Appendice A 79

Appendice B 83

Bibliografia 89

5

6

Introduzione

In questo lavoro di tesi si affronta il problema della costruzione di un

modello tridimensionale dinamico per una catena di elementi connessi

ed interagenti in grado di emulare il processo di ripiegamento di una

proteina noto come protein folding.

Tale processo, scoperto e investigato dalla biofisica sperimentale fin

dagli anni ’70, rappresenta ancora oggi una sorta di paradosso insolu-

to dal punto di vista teorico.

Stabilire il nesso tra la sequenza di amminoacidi che compongono una

data proteina e la configurazione spaziale (responsabile della specifica

attivita biologica) assunta al termine di ogni ripiegamento, e l’obiet-

tivo principale di fisici e biologi.

Cio rappresenterebbe uno strumento diretto per elaboare e sintetizzare

proteine con specifiche funzionalita biologiche a partire dalla codifica

(ricavabile dal DNA) delle sequenze amminoacide.

Nel Capitolo 1 si introduce il problema del protein folding dal punto

di vista fenomenologico: si descrivono gli aspetti salienti della fisica

della proteina, le sue strutture, i suoi costituenti, si delineano infine le

possibili metodologie di ricerca fra le quali, appunto, la modellistica

computazionale.

La prima questione da risolvere per la costruzione di un modello di-

namico e la creazione di un integratore che risolva il moto di quella

che possiamo per il momento definire una catena di aste connesse.

Questo e l’obiettivo principale di questa tesi.

Secondariamente l’integratore verra messo alla prova su di una cate-

na che, dotata degli opportuni potenziali di interazione, simulera un

filamento di amminoacidi in grado di ripiegare correttamente.

7

Per far sı che la distanza tra elementi contigui sia preservata im-

plicitamente, la catena viene fatta evolvere sotto l’azione di rotazioni

nello spazio.

Tali rotazioni saranno compiute tramite operatori matematici noti

come quaternioni.

Questi oggetti, la cui algebra fu investigata per la prima volta da Sir.

Hamilton nel secolo scorso, identificano una estensione quadridimen-

sionale dei numeri complessi; la loro composizione da una rappresen-

tazione non singolare delle rotazioni in R3 particolarmente efficiente

dal punto di vista computazionale.

Il Capitolo 2 e quindi una parentesi matematica che descrive i quater-

nioni e la loro applicazione per la manipolazione di oggetti nello spazio.

La particolare versatilita dei quaternioni nella programmazione e nel

calcolo gli ha resi uno degli strumenti maggiormente utilizzati per le

simulazioni di Dinamica Molecolare; a riguardo, nel Capitolo 3 verra

data una panoramica dei principali algoritmi di integrazione utilizzati

in questo campo e, genericamente, per la soluzione del moto di un cor-

po rigido. Maggiori dettagli ed esempi sono inclusi nelle Appendici.

Nel Capitolo 4 viene presentato l’integratore che verra poi utilizzato

nel corso delle simulazioni.

Per prima cosa si affronta il problema del moto di una singola asta;

ad ogni passo di integrazione le forze esterne modificano la velocita

angolare la quale, a sua volta, induce una rotazione planare della po-

sizione dell’asta compiuta attraverso la composizione di quaternioni.

L’algoritmo e strutturato con uno schema molto simile a quello de-

nominato leap-frog per la risoluzione di moti puramente traslazionali;

il passo temporale e suddiviso in tappe intermedie la cui composizione

fa raggiungere all’integratore una precisione di ordine O(h3 ), dove h

e il passo di integrazione.

Per estendere il problema da una ad N aste e necessario trovare il

contributo da sommarsi alle forze esterne dato dalle reazioni vincolari

che si propagano lungo la catena e che la mantengono connessa.

8

La stabilita dell’algoritmo e stata testata nel caso del comune pendolo

sferico e della catena di pendoli sferici nel campo uniforme di gravita.

Nel Capitolo 5 si opera la modellizzazione della proteina. Il varie-

gato e complesso mondo fenomenologico della chimica organica viene

drasticamente semplificato con la costruzione di potenziali “fittizi” che

riflettono gli aspetti piu rilevanti dell’interazione tra amminoacidi.

Lo stesso codice del filamento, composto in realta dai 20 tipi di ammi-

noacidi presenti in natura, viene ridotto alla composizione di 3 uniche

classi definite in base alla idrofobicita.

La presenza del solvente acquoso dentro al quale si svolgono gli esperi-

menti di ripiegamento della proteina viene simulata con una dinamica

browniana (o di Langevin) che riproduce gli urti casuali delle molecole

d’acqua e la dispersione energetica sotto l’effetto della frizione.

Alle equazioni differenziali deterministiche (o Lagrangiane) viene cosı

aggiunta una componente stocastica che unisce all’aspetto meccanico

del sistema il contributo termodinamico.

L’introduzione del cosiddetto bagno termico e stata suffragata da

prove di simulazione di moto browniano per le quali si riscontra una

buona attinenza fra i valori di oscillazione cinetica media e della dis-

tribuzione maxwelliana delle velocita registrati e le aspettative teoriche.

Ora abbiamo tutti gli ingredienti per procedere con le simulazioni di

ripiegamento.

Nell’ultimo Capitolo il problema del protein folding e ripreso ed anal-

izzato da un punto di vista teorico.

La trattazione cinematico-termodinamica di un sistema a molti gradi

di liberta capace di trovare autonomamente e, soprattutto, in un lasso

di tempo breve rispetto alla scala temporale, una e un’unica config-

urazione di equilibrio detta stato nativo presenta molteplici difficolta

concettuali.

Le simulazioni di protein folding compiute in questo lavoro sono state

eseguite su di un filamento relativamente breve (22 amminoacidi) la

cui configurazione nativa simile a una forcina e nota come β-sheet e

che rappresenta il prototipo di quello che possiamo definire un buon

9

ripiegatore.

In ogni simulazione tale proteina e immersa in un bagno termico a

temperatura costante.

Si e focalizzata l’attenzione sulla distribuzione dei punti di equilibrio

raggiunti dal sistema al termine dell’intervallo di tempo necessario al

ripiegamento e alla termalizzazione con l’ambiente.

Queste configurazioni vengono classificate in base alla loro distanza en-

ergetica (minimi del potenziale) e alla loro distanza configurazionale

dallo stato nativo; quest’ultima grandezza (χ) rappresenta una stima

per la similarita fra due elementi dello spazio 3N-dimensionale delle

configurazioni fornendo per essa un unico valore reale.

Dall’analisi statistica compiuta sui risultati di differenti serie di sim-

ulazioni indipendenti a temperatura costante emerge per la suddetta

sequenza la presenza di una configurazione di equilibrio di gran lunga

piu frequente e coincidente, come ci si puo aspettare, con lo stato na-

tivo β-sheet. Il relativo valore del potenziale coincide con il minimo

meccanico dell’energia del sistema.

Al crescere della temperatura nello spettro di popolazione compaiono

minimi situati in una regione configurazionalmente ed energeticamente

prossima allo stato nativo, che definiamo regione metastabile.

Il sistema all’equilibrio, sotto l’azione delle sollecitazioni termiche,

evolve visitando i minimi all’interno di detta regione come ci si puo

aspettare dalle ipotesi teoriche sull’attraversamento delle barriere di

potenziale che li separano.

Tutti i codici dei programmi di simulazione, archiviazione dei dati, sta-

tistica, immagini e animazione sono stati scritti in linguaggio C++;

per la parte grafica si e fatto uso delle librerie OpenGL.

10

Capitolo 1

Proteine e ripiegamento

1.1 Il problema del Protein Folding

La proteina e un eteropolimero formato da una catena di amminoacidi

la cui sequenza ordinata ne identifica la natura.

Sotto normali condizioni fisiologiche (solvente acquoso, pH neutro e

temperatura ambiente) la proteina ha la proprieta di ripiegarsi (fold-

ing) fino ad assumere una configurazione di equilibrio che ne determi-

na gran parte delle proprieta biochimiche; tale configurazione e detta

stato nativo.

Intervenendo dall’esterno la proteina puo essere denaturata con

agenti chimici o riscaldamento e, al ripristinarsi delle condizioni otti-

mali dell’ambiente, torna ad assumere la configurazione di stato nativo

(Figura 1.1).

Nel gergo tecnico la sequenza codificata di amminoacidi e denominata

struttura primaria, mentre la configurazione tridimensionale assun-

ta dal filamento e detta struttura terziaria.

Il legame diretto esistente tra struttura primaria e terziaria fu scop-

erto alla fine degli anni ’50 da Anfinsen [1] e da allora restano ancora

aperte due importanti problematiche riassumibili nelle domande:

• In che modo la struttura terziaria e codificata in quella primaria?

• Data una struttuta terziaria, quali sono le strutture primarie che

la ammettono come stato nativo?

11

Figura 1.1: Esempio di rappresentazione grafica della strutturaterziaria. Immagine tratta dal sito web del TCB Group, Universityof Illinois.

note con il termine di Protein Folding Problem e Inverse Folding

Problem.

Comprendere la dinamica soggiacente al processo del ripiegamento e

di grandissimo interesse biomedico per la costruzione (design) di nuove

proteine con funzioni biologiche desiderate a partire dalla conoscenza

della semplice sequenza di amminoacidi.

Dal punto di vista teorico cio significa comprenedere la causalita nello

schema

sequenza primaria −→ struttura terziaria −→ proprieta biologiche

Ottenere la struttura terziaria, ossia visualizzarne l’immagine tridi-

mensionale, infatti, comporta lunghe e costose operazioni di laborato-

rio come la cristallografia a raggi X o la risonanza magnetica nucleare

(NMR); per ricavare la struttura primaria invece ci si avvale di pro-

cessi molto meno dispendiosi attraverso la codifica del relativo gene di

DNA.

Quindi, ammesso di sapere come la struttura terziaria implichi deter-

minate proprieta biologiche, queste ultime, tramite la comprensione

del meccanismo del protein folding, sarebbero direttamente correlabili

a un codice sequenziale per il quale si dispone di un database vastis-

simo e in continua crescita.

12

Uno degli approcci per la risoluzione del problema del protein folding

si fonda sull’analisi diretta di questo database alla ricerca di corre-

lazioni statistiche tra struttura primaria e terziaria.

Per far questo algoritmi di elaborazione a reti neurali vengono previ-

amente “addestrati” mediante campioni gia noti.

Un secondo approccio, che sara quello esaminato in questo lavoro di

tesi, si basa sulla costruzione di un modello in grado di simulare il

comportamento della proteina.

1.2 La proteina

Le proteine sono catene costituite da amminoacidi connessi da legami

peptidici.

In natura sono presenti 20 varieta di amminoacidi ed essendo la pro-

teina tipicamente costituita da almeno un centinaio di amminoacidi,

esiste idealmente un potenziale esorbitante di possibili strutture pri-

marie; basti pensare che una generica catena di 100 elementi puo essere

generata in 20100 ≈ 10130 codici differenti.

La natura sotto la costante pressione evolutiva ha selezionato una pic-

colissima frazione di tali sequenze favorendo quelle a cui e associato

uno stato nativo stabile; in prima istanza chiameremo queste proteine

buoni ripiegatori.

La caratterizzazione di un buono o di un cattivo ripiegatore e una delle

tematiche cruciali per il problema del protein folding e verra discussa

nel Capitolo 6.

1.2.1 Amminoacidi

L’amminoacido e un composto organico formato da un atomo di car-

bonio centrale (Cα) al quale sono connessi un gruppo carbossilico, un

gruppo aminico, un gruppo H e una catena laterale R che varia con

13

il tipo di amminoacido.

Nella tabella 1.3 sono elencati i 20 amminoacidi con le relative catene

laterali che vanno dalla piu semplice della Glicina formata da un solo

atomo di idrogeno alle piu complesse come la Isoleucina o la Tripto-

fanina.

Solitamente la prima classificazione viene fatta in base alla idrofobic-

ita, una delle caratteristiche salienti in presenza di un solvente acqu-

oso. Gli amminoacidi lungo il filamento, infatti, tendono a disporsi

e impacchettarsi in modo da formare un nucleo (core) idrofobico che

minimizzi la superficie esposta al solvente ricoperto da una superficie

idrofilica.

Le conformazioni spaziali assunte dai gruppi attorno al nucleo Cα (det-

to residuo) sono quelle energeticamente favorite (dette rotameri) e

sono catalogate in vaste librerie. In generale ogni amminoacido puo

essere sintetizzato stabilmente in due conformazioni speculari : chi-

ralita L (sinistrorsa) e chiralita D (destrorsa) (Figura 1.2).

Per isotropia dello spazio amminoacidi con differenti chiralita han-

no proprieta chimico-biologiche identiche. Per legarsi chimicamente,

pero, sono necessari amminoacidi con la stessa chiralita. Cio compor-

ta l’ipotetica formazione di distinte proteine L e D; d’altronde anche

per una corretta interazione biologica fra filamenti distintie necessaria

la stessa chiralita.

In natura tutte le proteine sono composte da amminoacidi L e, non

essendovi alcuna ragione apparente, sembra esservi stata una sorta di

scelta evoluzionistica casuale che ha selezionato la chiralita sinistrorsa

per gli organismi di questo pianeta1.

I legami peptidici tra amminoacidi contigui lungo la catena avven-

gono tra l’atomo N del gruppo aminico che cede un ione H+ e l’atomo

C del gruppo carbossilico che cede lo ione OH- con la conseguente

condensazione di una molecola d’acqua.

1La esobiologia si interessa nella ricerca di amminoacidi con chiralita Dall’interno di meteoriti.

14

Figura 1.2: Schema strutturale dell’amminoacido nelle oppostechiralita.

Amminoacidi contigui rimangono cosı imperniati tramite questi lega-

mi sull’asse dei quali sono liberi di ruotare.

Come vedremo il modello standard tridimensionale schematizza ogni

amminoacido come un punto identificabile con l’atomo Cα e con un’

etichetta di appartenenza ad una delle tre classi: idrofobica (B), idro-

filica (L) e neutra (N), passando cosı da un codice a 20 lettere a un

codice a 3 lettere.

Questa drastica semplificazione permette di facilitare la modellizzazione

dei potenziali di interazione con l’introduzione del minor numero di

parametri necessari ad riprodurre le principali caratteristiche fisiche.

La scala di idrofobicita rappresenta grossolanamente la prima carat-

terizzazione fisica degli amminoacidi all’interno della sequenza e in

questo modello sara scelta come unica per la classificazione (Capitolo

5).

In generale le proprieta chimiche degli amminoacidi concorrono solo

nella formazione della struttura terziaria della proteina e non influen-

zano direttamente le sue qualita biologiche. Vale a dire: se due se-

quenze diverse formano la stessa struttura terziaria di fatto formano

la stessa proteina.

15

Figura 1.3: Codice, composizione, struttura e massa espressa in Dalton(g/mol) per i 20 tipi di amminoacidi

Figura 1.4: Struttura secondaria composta da tratti di α-elica (trattoa zig-zag) e β-sheet (freccia) associata ad una sequenza primaria nelcodice a 20 lettere.

16

1.2.2 Struttura secondaria

La sequenza degli amminoacidi e caratterizzata su piccola scala dalla

comparsa di conformazione ricorrenti quali la α− elica e la β− sheet

composti approssimativamente da una decina di elementi. Come raf-

figurato in Figura 1.4 il codice della struttura primaria della proteina

puo essere spezzato e raggruppato in blocchi che costituiscono la strut-

tura secondaria.

La struttura terziaria della proteina puo quindi essere schematizzabile

con la composizione di oggetti ricorrenti orientati (eliche o strisce)

connessi da tratti di filamento (loop) non riconducibili ad esse.

E importante sottolineare come gia in questo primo passo di aggregazione

non vi sia un nesso diretto fra sequenza di amminoacidi e α -eliche o

β -sheet, vale a dire tra sequenza primaria e secondaria.

Sequenze diverse di amminoacidi possono dare luogo localmente alle

stesse strutture secondarie e concorrere globalmente a diverse terziarie.

La struttura secondaria rappresenta percio un ponte tra la primaria

e la terziaria; di notevole aiuto per il riconoscimento e la classificazione

delle proteine ci offre una base di partenza o un test di prova per

modelli e simulazioni con sequenze brevi.

Come si vedra nel Capitolo 6, la sequenza campione utilizzata per le

simulazioni di ripiegamento sara composta da 22 amminoacidi il cui

stato nativo corrisponde alla β -sheet rappresentata in Figura 6.2 (a).

Oltre la struttura terziaria, su scale piu estese, la combinazione di

due o piu filamenti di amminoacidi da luogo a complessi intrecci di

filamenti a cui viene dato il nome di struttura quaternaria.

1.3 Interazioni fisiche

Per dare una interpretazione dinamica del processo di ripiegamento e

necessario individuare le forze che determinano la struttura terziaria.

I legami che connettono i diversi gruppi dell’amminoacido e il legame

peptidico tra amminoacidi contigui lungo la catena sono tutti covalen-

17

ti e rappresentano un contributo costante che si cancella calcolando

differenze energetiche tra diverse configurazioni.

Le interazioni non covalenti restanti sono tre ordini di grandezza piu

deboli e sono i motori effettivi nella dinamica del folding.

Schematizzando:

• Repulsione a corto range: fra coppie di atomi che si avvicinano

al punto da sovrapporre i rispettivi orbitali elettronici; la repul-

sione cresce enormemente e mantiene gli amminoacidi a debita

distanza.

• Forze elettrostatiche : fra coppie di atomi parzialmente ioniz-

zati in accordo con le leggi Coulombiane e con le modulazioni

dielettriche del mezzo circostante.

• Interzioni van der Waals : mutuate dagli effetti di polarizzazione

indotta.

• Ponti a idrogeno : fra coppie di atomi elettronegativi; sono re-

sponsabili della crezione di strutture locali come le α - eliche e

le β - sheet.

• Interazione idrofobica : generata dai dipoli presenti nel solvente;

l’acqua trasmette forze effettive fra atomi non polari e induce la

creazione di un nucleo idrofobico.

1.4 Un primo sguardo al modello

Diamo ora una breve introduzione sugli aspetti legati alla modelliz-

zazione che verra discussa approfonditamente nel Capitolo 5.

La struttura terziaria e individuata dalle coordinate spaziali dei residui

Cα avendo cosı ridotto gli amminoacidi a oggetti puntiformi.

Ogni residuo e vincolato a mantenersi a distanza fissa da quello che lo

precede nella sequenza; la sua posizione relativa e percio individuata

dalle due coordinate angolari sulla sfera di raggio costante centrata

18

nel residuo precedente.

Modelli semplificati studiano la proteina discretizzando le possibili po-

sizioni angolari del generico residuo rispetto al precedente; assumendo

che tali posizioni siano ad esempio 7, una proteina costituita da 100

amminoacidi ha a disposizione 7100 possibili configurazioni e fra tutte

queste, al variare delle condizioni iniziali e dell’ambiente esterno, e

in grado di assumere, tramite il ripiegamento, la configurazione dello

stato nativo.

Lo stato nativo e presumibilmente il minimo assoluto in un panorama

energetico che dipende da un elevato numero di gradi di liberta.

La simulazione computazionale ha lo scopo di gettare luce sulla con-

formazione di questo panorama (presenza di minimi locali, loro dispo-

sizione e stabilita, presenza di traiettorie preferenziali che facilitino il

ripiegamento) in relazione alla sequenza assegnata.

Il problema e strutturalmente inaffrontabile dal basso, ossia a partire

dai costituenti atomici dei singoli amminoacidi; la complessita del sis-

tema ne fa perdere il controllo, inoltre il divario di scala temporale che

contrappone processi locali e globali rende necessario un adattamento

dei passi di integrazione a tempi ridottissimi (≃ 10−15s) aggravando

ulteriormente il carico computazionale.

Gli amminoacidi, in quella che si puo definire una modellistica stan-

dard sono percio i costituenti primi della proteina; masse puntiformi

vincolate e soggette a potenziali configurazionali creati per emulare le

reali interazioni esistenti.

Ammettendo di poter disporre dei potenziali da introdurre nel mod-

ello, il primo problema e quello di costruire un integratore dinamico

efficiente in grado di descrivere l’evoluzione temporale di un sistema

di punti vincolati.

Il problema del vincolo puo essere aggirato sostituendo le ipotetiche

aste rigide che connettono i residui con oscillatori armonici ad alta fre-

quenza; in questo caso il problema di una catena di N residui e ridotto

al problema di N punti nello spazio soggetti alla dinamica newtoniana;

ci si trova pero di fronte al problema di dover innalzare fortemente la

19

rigidita elastica delle molle e conseguentemente ridurre il passo di in-

tegrazione per assecondarne la dinamica veloce a spese dell’efficienza

di calcolo.

D’altro canto considerando le aste rigide e necessario avvalersi di un

integratore che preservi il vincolo delle distanze tra residui contigui.

La costruzione di un integratore con queste caratteristiche e uno dei

principali scopi di questo lavoro.

Un’ altra classe di modelli detta a reticolo (o lattice models) elude

questo problema schematizzando le posizioni dei residui all’interno di

una griglia tridimensionale di caselle; le posizioni angolari di residui

contigui sono cosı discretizzate e la dinamica assume un’ impronta

probabilistica con enormi agevolazioni computazionali.

Da questi modelli si possono ottenere importanti risultati pur se la

marcata semplificazione pone limiti al parallelismo con la meccanica

della proteina reale.

Estendendo la dinamica ad uno spazio continuo (fuori-reticolo, off-

lattice) il modello di riferimento per molti lavori ([20], [21], [19]) svolti

negli ultimi anni e quello di Thirumalai - Honeycutt [18] che utilizzer-

emo ed esamineremo dettagliatamente nel Capitolo 5.

Per la conservazione delle distanze l’algoritmo di integrazione fara uso

di rotazioni nello spazio. Tali rotazioni, come e stato detto nell’Intro-

duzione, verranno compiute attraverso operatori vettoriali noti come

quaternioni.

Nel prossimo Capitolo verra esposta tale algebra e vaerra chiarito

il motivo di un cosı largo sviluppo nel campo del calcolo e della

elaborazione grafica tridimensionale.

20

Capitolo 2

Quaternioni

Nel 1843 Sir William Hamilton presento alla Royal Irish Academy

una ricerca personale in cui veniva proposto un tentativo di generaliz-

zazione dei numeri complessi come strumento di manipolazione della

geometria spaziale.

L’idea di Hamilton trae spunto dall’isomorfismo esistente tra il grup-

po (C, ·) , campo complesso e moltiplicazione, e le rotazioni sul piano.

Un numero c ∈ C, formato da parte reale e immaginaria (x, ıy) iden-

tifica un vettore bidimensionale che giace sul piano complesso; alla

moltiplicazione tra elementi e associata la rotazione su tale piano dei

rispettivi vettori.

Ricordiamo brevemente che tramite la notazione di Eulero ogni

numero complesso e esprimibile come

c = a+ ıb = R cos(θ) + ıR sin(θ) = Reıθ (2.1)

e il prodotto di una coppia di complessi come

c · c′ = R · R′eı(θ+θ′) (2.2)

Se R′ = 1, ossia c′ e unitario, il prodotto corrisponde ad una rotazione

del vettore c (Figura 2.1).

Hamilton si propose di trovare un gruppo (Q, ) che estendesse

questo isomorfismo alle rotazioni nello spazio tridimensionale.

c ∈ C ⇐⇒ Rotazioni 2D

21

C’ φ’

φ

φ’C*C’

Re

Im

C

Figura 2.1: Al prodotto tra numeri complessi corrisponde unarotazione sul piano.

q ∈ Q ⇐⇒ Rotazioni 3D

Che dimensionalita dare agli elementi di Q ?

Se alla rotazione 2D e associato un binomio (x, ıy), penso Sir Hamil-

ton, a quelle 3D deve essere associato un trinomio. Bisognava dotare

l’elemento q di una terza parte oltre a quella reale e immaginaria e

costruire quindi un corredo di regole algebriche consistenti per definire

il gruppo (Q, ). Sir Hamilton non riuscı a venirne a capo e concluse,

oltre al fatto che la topologia gioca brutti scherzi, che un trinomio non

era sufficiente a suoi scopi.

Associato alle rotazioni 3D vi e uno spazio Q popolato da quadrinomi

composti da una parte reale e tre distinte parti immaginarie ı, , k.

A tale quadrinomio Hamilton diede il nome di quaternione. 1

1In una lettera ad un collega Sir Hamilton racconta di come, dopo una decinad’anni di laboriose riflessioni, fu folgorato dall’idea risolutiva mentre camminavacon la moglie lungo il Royal Canal di Dublino. Preso dall’eccitazione, incise conun temperino la nota formula ı2 = 2 = k2 = ı k = −1 su una pietra del ponteche successivamente i suoi alunni ribattezzarono Quaternion Bridge.

22

2.1 Algebra dei quaternioni

Lo spazio Q e generato da una unita reale (1, oppure omesso) e

tre unita immaginarie ı k. Un elemento di tale spazio e percio

univocamente individuato da quattro variabili indipendenti.

q = q0 + q1ı+ q2+ q3 k = (q0, q1, q2, q3)

La somma e il prodotto esterno per uno scalare che compaiono

nella suddetta espressione sono definite in analogia con l’algebra or-

dinaria; (Q,+, ·) risulta percio uno spazio vettoriale isomorfo a R4.

Introduciamo ora il prodotto interno tra elementi di Q costruendolo a

partire da regole algebriche fondamentali con cui si combinano le tre

unita immaginarie.

ı2 = 2 = k2 = −1 (2.3a)

ı · = − · ı = k (2.3b)

· k = −k · = ı (2.3c)

k · ı = −ı · k = (2.3d)

Con riferimento a tali regole il prodotto tra due elementi

q = (q0, q1, q2, q3), p = (p0, p1, p2, p3)

risulta essere una combinazione lineare delle rispettive coordinate

q p = (q0p0 − q1p1 − q2p2 − q3p3)

+ i(q0p1 + q1p0 + q2p3 − q3p2)

+ j(q0p2 − q1p3 + q2p0 + q3p1)

+ k(q0p3 + q1p2 − q2p1 + q3p0)

(2.4)

esprimibile nella forma matriciale

q p =

q0 −q1 −q2 −q3q1 q0 −q3 q2q2 q3 q0 −q1q3 −q2 q1 q0

p0

p1

p2

p3

(2.5)

23

=

p0 −p1 −p2 −p3

p1 p0 p3 −p2

p2 −p3 p0 p1

p3 p2 −p1 p0

q0q1q2q3

(2.6)

Si dimostra facilmente che l’ operazione binaria interna in Q del

prodotto risulta associativa e distributiva rispetto alla somma, ma,

come si osserva nelle eq. 2.5 e 2.6, incontrovertibilmente non com-

mutativa (ne anticommutativa). A ragione di cio si ricordi che la

commutativita di cui godono le rotazioni 2D (Equazione 2.2) e una

proprieta che viene intrinsecamente persa nelle rotazioni 3D dove il

risultato di una combinazione di rotazioni dipende dall’ ordine con

cui sono concatenate.

Sebbene le quattro variabili qi i = 0 . . . 3 siano indipendenti e di

conseguenza l’ugulianza tra due elementi di Q implichi la simultanea

ugualianza di quattro quantita reali, risulta utile identificare il gener-

ico quaternione q come l’unione di una parte reale (q0) e di una parte

vettoriale rappresentata dal trinomio immaginario

iq1 + jq2 + kq3 = ~q

Si noti come la curiosa coincidenza di simbolismo tra le unita immag-

inarie ı, , k e i consueti versori euclidei i, j, k getti le fondamenta del

dualismo tra un’ algebra estesa dei complessi e lo spazio tridimension-

ale.

Dall’ algebra di C possiamo importare l’operatore di coniugazione

estendendolo alle tre dimensioni immaginarie

q = (q0,−q1,−q2,−q3) = (q0, −~q) (2.7)

e introdurre quindi la norma di un quaternione definita come

‖q‖2 = q q = q20 + q2

1 + q22 + q2

3 (2.8)

Se ‖q‖2 = 1 il quaternione e detto unitario.

Il campo (Q,+, ·) dotato della metrica d(q, p) =√

‖q − p‖ rappresenta

un’algebra di Banach.

24

Mettiamo ora in corrispondenza lo spazio R3 con lo spazio Q de-

finendo per il generico vettore ~x ∈ R3 il corrispettivo quaternione

x = (0, ~x) e dato l’omeomorfismo esistente tra R3 e il sottospazio

vettoriale

QR3 = q ∈ Q | q0 = 0

ci permettiamo di identificare concettualmente e notazionalmente il

quaternione x ∈ QR3 con il vettore ~x.

Ora abbiamo tutti gli ingredienti per illustrare come l’algebra quater-

nionale racchiusa nelle leggi di combinazione 2.3 possa essere una rap-

presentazione non singolare del gruppo delle rotazioni SO(3) 2.

Enunciamo il seguente

Teorema 1 Dato un vettore ~x ∈ R3 e un quaternione unitario q

q ~x q e ancora un elemento di R3

e la trasformazione

~x→ q ~x q (2.9)

rappresenta un rotazione.

Seguendo la traccia di Hamilton, si puo verificare che ogni quaternione

unitario puo essere associato ad un versore v e ad un angolo α nella

forma

qv,αdef= (cos

α

2, v · sin α

2) (2.10)

Come corollario al Teorema 1 vale quindi

Corollario 1 Dato ~x ∈ R3 e qv,α

~x → qv,α ~x qv,α (2.11)

rappresenta una rotazione antioraria del vettore ~x compiuta rispetto

al versore v di un angolo α.

La composizione di due rotazioni si ottiene, per proprieta associa-

tiva, dal prodotto di due quaternioni

q2(q1~x q1)q2 = (q2q1)~x (q1 q2)

2Si puo generalizzare l’algebra dei quaternioni assumendo che esista una base1, e1, e2, e3 che soddisfi la condizione ei · ej = −δij + ǫijkek.

25

j

k

i

v

x

x’

α

q = cos(α/2) + v sin(α/2)

x’ = q x q

Figura 2.2: Al prodotto tra quaternioni e associata una rotazione nellospazio ~x→ ~x′ di un angolo α attorno all’asse individuato da v.

2.1.1 Osservazioni

D’ ora in avanti, dove non vi siano ambiguita, verra omesso l’operatore

.Per prima cosa riscriviamo le Eq. 2.5 e 2.6 in notazione vettoriale.

Dati q = (q0, ~q) e p = (p0, ~p) il loro prodotto si esprime

r = (r0, ~r) dove

r0 = q0 · p0 − ~q · ~p

~r = q0 · ~p+ p0 · ~q + ~q × ~p

(2.12)

dove sono utilizzati il prodotto scalare e vettoriale ordinari.

• Si puo notare come q p = p q ⇔ ~q ‖ ~p. A corollario di cio si

ottiene che due rotazioni commutano solo se posseggono un asse

comune.

qv,α qv,β~x qv,β qv,α = qv,β qv,α~x qv,α qv,β (2.13)

• Cambiando il versore v con il suo opposto e l’angolo α con il suo

complemento a π si ottiene il quaternione q−v,π−α = −qv,α. Di

conseguenza ogni coppia q,−q da luogo alla stessa rotazione,

proprieta comune all’algebra delle matrici di Pauli e in generale

caratteristica di quelle che vengono dette skew algebra

26

v

x x’

v xv x

α

Figura 2.3: Rotazione planare (v ⊥ ~x) di un angolo α. ~x′ = ~x cosα +(v × ~x) sinα.

• Un caso particolare di rotazione si ha quando essa e planare ossia

quando vi e ortogonalita tra vettore oggetto e asse di rotazione;

vale il seguente

Teorema 2 1 Dati un vettore ~x ∈ R3 e un quaternione unitario qv,α

con v ⊥ ~x si ha

qv,α~x qv,α = qv,2α ~x (2.14)

Battezziamo la 2.14 espressione ridotta per la rotazione.

Dimostrazione Sviluppiamo brevemente l’algebra dei quaternioni

per una rotazione planare .

qv,2α ~x = (cosα, v · sinα) · (0, ~x)= (−~xv sinα, ~x cosα+ (v × ~x) sinα)

= (0, ~x cosα + (v × ~x) sinα)

(2.15)

che, come mostrato in Figura 2.3, rappresenta appunto la rotazione

antioraria di ~x rispetto all’asse v.

• Osserviamo infine che applicando l’espressione ridotta nel caso

in cui v 6⊥ ~x continua a valere la conservazione della norma

‖qv,2α~x‖ = ‖~x‖ (2.16)

27

ma poiche

(qv,2α~x)0 6= 0 ⇒ qv,2α~x /∈ R3 (2.17)

essa cessa di rappresentare una rotazione in R3.

2.2 Applicazioni dei quaternioni

Dal punto di vista teorico la quadridimensionalita di questi oggetti

trova vaste applicazioni in fisica moderna: in relativita generale per le

assonanze con la geometria di Minkowsky, in meccanica quantistica si

sfrutta l’isomorfismo tra quaternioni e matrici di Pauli per lo svilup-

po di teorie unificatrici di campo forte e elettrodebole; le equazioni

di Maxwell ad esempio si scrivono come un’unica equazione d’onda

quaternionale non omogenea e l’ equazione di Klein Gordon come un

semplice oscillatore armonico quaternionale; in matematica tramite

successioni iterative di quaternioni si possono generare insiemi frattali

tridimensionali (il set di Julia consueto e una figura 2D e si genera con

una successione di numeri complessi)

Dal punto di vista tecnico e tecnologico, i quaternioni sono un ve-

loce strumento di calcolo per la manipolazione geometrica dello spazio

tridimensionale.

Elenchiamo dispositivi di controllo nella robotica, percezione e ri-

conoscimento di immagini per l’intelligenza artificiale e la neuroscien-

za, dinamica dei vortici nelle simulazioni meteorologiche, altimetri,

giroscopi e altra strumentazione a bordo di aereoplani, elicotteri e,

con mio grande rammarico, missili teleguidati; triangolazioni per il

sistema satellitare GPS, reti radar e ingegneria aereospaziale.

Infine c’e il grande supporto dato alla computer grafica: dalle simu-

lazioni di dinamica molecolare agli effetti speciali di Hollywood fino

alle evoluzioni turbinanti della visione soggettiva in videogiochi come

Tomb Raider.

Questo sodalizio tra quaternioni e calcolo computazionale e dovuto

principalmente a due fattori.

1. La rappresentazione non singolare del gruppo delle rotazioni

28

rende possibile la descrizione della geometria sferica tramite una

sola mappatura, evitando cosı la delicata operazione del cambio

di carta necessaria in presenza di singolarita.

2. La rotazione, espressa in forma quaternionale, e una combi-

nazione di somme e moltiplicazioni le quali, confrontate con il

carico computazionale necessario per gli sviluppi trigonometrici,

rendono piu agevole il calcolo.

2.3 Dinamica Molecolare e Quaternioni

La dinamica molecolare (MD) descrive un gas di molecole consideran-

dole alla stregua di corpi rigidi.

Di conseguenza all’ integrazione del moto traslazionale va sovrapposta

quella del moto rotazionale per ogni singola molecola.

Come sappiamo quest’ ultimo e intrinsecamente piu complicato del

primo a causa della dipendenza posizionale delle velocita angolari e

dalla non-separabilita della relativa Hamiltoniana.

Percio se si vuole ottenere e simulare al calcolatore l’evoluzione di un

gran numero (∼ 106) di molecole e naturale pensare che l’algoritmo

integratore del moto debba essere veloce.

I principali metodi di integrazione per il moto del corpo rigido in

dinamica molecolare sono parametrizzati in termini di

• Angoli di Eulero

• Quaternioni

• Matrice degli assi principali

• Vincoli atomici

Quest’ultimo tratta la dinamica del corpo rigido come la composizione

di moti puramente traslazionali dei singoli atomi (puntiformi) che si

muovono sotto l’azione delle forze potenziali e delle le forze vincolari

introdotte per mantenere costanti le distanze interatomiche. Tale ap-

proccio e riconosciuto essere il migliore in termini di stabilita ed e

29

per questo considerato come il modello di confronto per lo sviluppo di

nuovi integratori. Lo svantaggio risiede nel carico di calcolo indotto

dai complicati sistemi non lineari da risolvere ad ogni passo temporale

per ogni molecola e dal conseguente calo di efficienza al crescere del

numero di atomi per molecola.

L’integrazione numerica con angoli di Eulero e poco efficiente a causa

delle singolarita presenti nelle equazioni che impongono cambi di car-

ta durante il processo di calcolo. Da questo inconveniente sono esenti

l’approccio matriciale e quello quaternionale che sono entrambi uti-

lizzati in dinamica molecolare come alternativa al metodo dei vincoli

interatomici. Il moto viene parametrizzato tramite rispettivamente 6

e 4 variabili indipendenti , ma l’integrazione deve affrontare il prob-

lema della conservazione del vincolo di ortonormalita per le matrici e

di norma unitaria per i quaternioni.

Questo vincolo puo essere imposto forzatamente ad ogni passo di inte-

grazione tramite riscalature o proiezioni, oppure essere implicitamente

garantito utilizzando trasformazioni che lo conservano.

Una discussione approfondita dei diversi integratori verra fatta nel

prossimo Capitolo e in Appendice B.

2.3.1 Riferimento per un corpo rigido

L’orientazione di un corpo rigido e univocamente determinata dall’ori-

entazione del sistema ad esso solidale esprimibile con una matrice 3×3

ortonormale A che ha per colonne i versori degli assi principali

Generalizzando la rotazione quaternionale di un vettore

~x→ q ~x q

a quella di una matrice

A → qA q

si puo individuare l’orientazione spaziale di un riferimento mobile che

abbia come versori di base la terna (e1, e2, e3) rispetto ad un sistema

fisso di base (i, j, k) tramite il quaternione q.

30

A =

| | |e1 e2 e3| | |

= q

| | |i j k| | |

q =

| | |q i q qj q qk q| | |

= q I q

(2.18)

Istante per istante sia la matrice A che il quaternione q individuano

l’orientazione del corpo rotante, ma mentre per A sono necessarie 6

variabili indipendenti, utilizzando q il moto puo essere parametrizzato

dalle sue 4 coordinate integrando il sistema di equazioni differenziali

q =1

2q · ω (2.19)

accoppiato con il sistema di equazioni di Eulero per le componenti di

ω

ωx =τ bx

Ixx+

(

Iyy − Izz

Ixx

)

ωx ωz x→ y → z (2.20)

La soluzione a questo sistema (non lineare) e vincolata a rimanere

sulla sfera unitaria S3, ossia si deve verificare

‖q(t)‖ = 1 ∀t (2.21)

Il mantenimento di questo vincolo costituisce una delle questioni spinose

per un algoritmo di integrazione quaternionale e verra affrontata e

discussa in Appendice B.

2.4 I quaternioni nel modello per la pro-

teina

Come e stato accennato nell’Introduzione la catena polimerica di am-

minoacidi che costituisce la proteina e modellizzata come una succes-

sione di oggetti puntiformi connessi da legami rigidi. L’intento sara

quello di riuscire a trattare singolarmente ognuno di questi oggetti.

Potremmo definire tale sistema corpo semirigido, ossia un set di corpi

rigidi vincolati a mantenere un punto in comune con i contigui (le es-

tremita dove risiedono gli amminoacidi); nella fattispecie questi corpi

rigidi sono puntiformi, ma non adimensionali in quanto il vincolo li

31

rende idealmente associabili ad una asta monodimensionale (sottile)

la cui massa e concentrata ad un’ estremita.

Mettiamo in luce un tratto saliente della questione: la trattazione

dinamica illustrata in Sezione 2.3.1 non si presta al caso di un asta

monodimensionale. Tale oggetto rappresenta infatti un corpo rigido

degenere e la sua patologia, racchiusa nella sua monodimensionalita,

si manifesta nelle seguenti constatazioni che si implicano mutuamente

• L’asta non individua un sistema di riferimento, poiche ad essa e

associato un vettore, non una terna.

• La matrice di inerzia della singola asta e singolare.

• Sono sufficienti due coordinate per individuare l’orientazione di

una asta. 3

• Il vettore asta e invariante rispetto ad una rotazione (o rivoluzione)

attorno al suo asse.

• La posizione dell’asta e (supposta unitaria) rispetto ad un ver-

sore di riferimento (ad esempio k) non e univocamente associata

ad un solo quaternione, ossia esistono infiniti q tali che

q k q = e

dove e e il vettore asta.

Di conseguenza per la dinamica di questa catena di asticelle non fare-

mo uso di equazioni del tipo 2.19; la configurazione del sistema verra

descritta da un set di vettori e ∈ R3.

Ai quaternioni verra dato il compito operatoriale di aggiornare tale

configurazione ad ogni passo di integrazione, senza che essi assumano

il ruolo di variabili dinamiche.

Per risolvere la dinamica di una singola asta procederemo nel Capitolo

4 analizzando per prima cosa il caso del pendolo sferico e, successiva-

mente, di una catena di pendoli connessi.

3Utilizzando i consueti angoli di Eulero di precessione, nutazione e rivoluzione

(ϕ, θ, ψ), per parametrizzare la posizione dell’asta sono sufficienti due angoli(ϕ, θ, ·)

32

Capitolo 3

Dinamica Molecolare ealgoritmi di integrazione

3.1 Introduzione

Gran parte dei problemi di fisica chimica, biochimica e biologia pre-

suppongono la risoluzione di esperimenti al calcolatore che simulino

processi che possiamo generalmente definire di dinamica molecolare

(MD); il gas o il liquido in esame e normalmente studiato come una

collezione di corpi rigidi in interazione. Anche se la trattazione di tale

problema ha origini tutt’altro che recenti, lo sviluppo di un algoritmo

di integrazione stabile ed efficiente richiama ancora oggi gli sforzi di

fisici e matematici.

3.2 Generalita sugli integratori

Le svariate metodologie di integrazione si qualificano principalmente

in base all’ ordine della precisione e al peso del calcolo.

Altri punti di forza sono la stabilita, ossia la gestione di eventuali

situazioni critiche che eviti l’incorrere di singolarita numeriche e, in

visione di un vasto e diversificato utilizzo, la versatilita e la adatta-

bilita dell’algoritmo alla conduzione di diversi esperimenti (energia,

temperatura, pressione costanti, ecc...).

Con particolare riguardo ai sistemi vincolati, un algoritmo puo garan-

33

tire il mantenimento del vincolo in modo implicito o esplicito. Nel

primo caso il ciclo del passo di integrazione e strutturato in modo

da preservare intrinsecamente il vincolo; nel secondo caso al termine

di ogni ciclo viene compiuta una (piccola) modifica delle variabili di-

namiche che vengono forzatamente proiettate sulla varieta descritta

dalle equazioni di vincolo.

Come esempio particolare, per un sistema hamiltoniano e possibile

trovare un integratore simplettico 1 che ha la qualita di conservare

implicitamente tutti gli invarianti geometrici del moto (invarianti di

Cartan). Di conseguenza in un sistema che ammette integrali primi la

loro conservazione e garantita a meno di fluttuazioni numeriche pro-

porzionali al passo di integrazione.

Precisione, velocita di calcolo, stabilita, simpletticita e mantenimen-

to implicito dei vincoli sono qualita non necessariamente esclusive;

trovare un algoritmo con determinati requisiti, rappresenta la sfida

matematica alla base della progettazione di un integratore.

3.3 Metodologie di integrazione

In questa sezione viene presentata una panoramica qualitativa delle

piu comuni strategie algoritmiche in campo molecolare.

Per un approfondimento matematico specifico si rimanda all’ Appen-

dice A e B.

3.3.1 Approccio atomico

Il metodo tutt’ora largamente piu utilizzato in MD e quello dei vin-

coli atomici; ogni molecola del sistema e un aggregato di atomi (pun-

tiformi) la cui traiettoria nello spazio delle fasi e integrata secondo

le equazioni di Newton; la conservazione della struttura molecolare e

1la matrice di trasformazione S ad ogni passo temporale verifica la condizionedi simpletticita SJS

T = J.

34

garantita dall’introduzione di forze vincolari che mantengono costanti

le distanze interatomiche ad ogni passo di integrazione.

Nell’approccio atomico la maggior parte degli algoritmi (il piu uti-

lizzato dei quali e il metodo SHAKE [8]) si avvale dell’aggiunta di

termini nell’Hamiltoniana fattorizzati da moltiplicatori di Lagrange.

A questo punto lo schema di integrazione puo essere ricondotto al co-

mune metodo Verlet (con le sue varianti velocity-Verlet e leap-frog) che

rappresenta lo standard per l’evoluzione di moti puramente traslazion-

ali (vedi Appendice A).

Naturalmente l’aggiunta delle condizioni vincolari accresce il numero

di equazioni, generalmente non lineari, e il peso di calcolo computazionale

rendendo inefficiente questa strategie per strutture molecolari molto

complesse.

Per aggirare questa difficolta Ahlrichs e Brode [10] sviluppano un

metodo (detto degli assi principali) con il quale la generica struttura

atomica viene sostituita da 4 pseudo-particelle (3 per molecole planari,

2 per molecole lineari) rigidamente ancorate agli assi del corpo; fatto

cio il moto e integrato come sopra.

3.3.2 Approccio molecolare

Le limitazioni imposte dalla strutturalita molecolare nel caso dei vin-

coli atomici sono assenti in una trattazione dinamica della molecola

come corpo rigido; l’integrazione delle relative equazioni garantisce

implicitamente la conservazione delle distanze interatomiche e il moto

si traduce nella composizione di una traslazione e di una rotazione

rispetto al baricentro.

Il moto traslazionale del centro di massa e generalmente integrato con

metodo Verlet o affini.

L’applicazione di questi schemi presuppone che la derivata posizionale

sia indipendente dalla posizione stessa o similmente che l’accelerazione

sia indipendente dalla velocita.

Questa condizione e intrinsecamente violata nel caso delle variabili

35

di orientazione, siano esse angoli di Eulero, matrici di rotazione o

quaternioni, dove vi e un esplicita dipendenza dell’accelerazione ango-

lare dalla velocita.

Per l’orientazione di un corpo rigido, come si e detto nel Capitolo 1,

sono possibili diverse parametrizzazioni;

• angoli di Eulero

• quaternioni

• matrice di rotazione

Angoli di Eulero La parametrizzazione diretta tramite angoli di

eulero utilizza il numero minimo di 3 coordinate che, al contrario delle

6 coordinate della matrice di rotazione e le 4 del quaternione, soggette

rispettivamente ai vincoli di ortonormalita e di norma unitaria, risul-

tano non vincolate.

Il principale inconveniente che emerge da questo tipo di rappresen-

tazione sono le singolarita presenti nelle zone polari a cui si deve far

fronte mediante un adeguato cambio di mappa (ad esempio una inver-

sione degli assi).

E’ stato proposto recentemente [11] un metodo di splitting completa-

mente simplettico e reversibile in cui l’hamiltoniana viene separata in

quattro termini e tramite il quale il moto risultante del corpo rigido

si scompone nella concatenazione di rotazioni planari.

Quaternioni Nell’ambito dell’approccio molecolare alla MD, la parametriz-

zazione quaternionale proposta da Evans [12] per i moti rotazionali del

corpo rigido e largamente la piu diffusa.

Come si e appena detto non e possibile utilizzare lo schema Verlet

per l’evoluzione delle 4 variabili qi; l’hamiltoniana in rappresentazione

quaternionale risulta inoltre non separabile ed e quindi necessario l’u-

tilizzo di schemi non simplettici come Runge-Kutta o Gear predictor-

corrector, solitamente al quarto ordine.

36

Un metodo di integrazione alternativo e stato introdotto da Fincham

[13] e rappresenta il primo tentativo di trasportare lo schema leap-frog

in campo rotazionale meritandosi l’appellativo di leapfrog-like.

Il metodo Fincham risulta efficiente nel caso di simulazioni a tem-

peratura costante, ma le fluttuazioni energetiche nelle simulazioni a

energia costante sono molto grandi paragonate a quelle implementate

con i vincoli atomici.

Inoltre in questo metodo non vi e conservazione implicita della norma

del quaternione il quale viene sottoposto ad una forzosa procedura di

riscalatura ad ogni ciclo perche vi sia rigidita molecolare.

Recentemente un nuovo metodo a velocita angolari leapfrog-like e sta-

to sviluppato da Omelyan [7]. Lo schema puo essere adattato sia a

quaternioni che a matrici di rotazione come variabili di orientazione e

prevede una conservazione implicita delle norme.

L’algoritmo originale prevede la risoluzione iterativa di un sistema

non lineare ad ogni ciclo che, pur essendo computazionalmente meno

onerosa di quella necessaria per ottenere le forze vincolari nell’approc-

cio atomico, rappresenta un inconveniente per simulazioni a lunga du-

rata.

Successivamente lo stesso Omelyan ha presentato una variante modi-

ficata del suo integratore [14] nella quale il sistema non lineare viene

risolto semi-analiticamente e che mostra eccellente stabilita nelle sim-

ulazioni a energia costante.

Approfondimenti tecnici sui metodi di integrazione Evans, Fincham e

Omelyan citati sono riportati in Appendice B.

Matrice di rotazione Kol presenta infine un algoritmo [15] in cui

la matrice di rotazione non e parametrizzata, ma le 6 coordinate sono

fatte evolvere direttamente.

L’approccio che viene fatto e di tipo Hamiltoniano per cui alla matrice

di rotazione che individua l’orientazione del corpo viene affiancato un

37

momento coniugato corrispondente.

Alla Hamiltoniana che ne deriva viene aggiunto un termine che coin-

volge un moltiplicatore di Lagrange matriciale per assicurare l’ortonor-

malita della matrice e la rigidita della molecola.

Tale metodo in assonanza con SHAKE e stato battezzato R-SHAKE

(Rotational Shake) ed e approfonditamente discusso in Appendice B.

3.4 Dinamica e integratori per una cate-

na

Come e stato discusso nel Capitolo 1 i composti polimerici come la pro-

teina e il DNA possono essere trattati meccanicamente come catene

di elementi incernierati.

Esiste gia una vasta gamma di studi teorici riguardanti la dinamica di

corpi rigidi connessi [16] e applicati sull’assemblaggio e l’articolazione

di bielle, perni e pistoni per scopi ingegneristici.

L’approccio [17] si basa sulla traduzione del concetto di contatto in

un set di forze e momenti vincolari in grado di mantenere tale con-

tatto. Per prima cosa i tipi di contatti (punto comune, asse comune,

superficie comune, punto su asse, punto su superficie ecc...) vengono

trasformati in equazioni lineari che vincolano le variabili del moto.

Le forze e momenti trasmessi dalle strutture rigide dell’oggetto sono

le soluzioni del sistema che ne deriva.

Nel modello standard per la proteina Honeycutt-Thirumalai si e in

presenza di una collezione di aste con mutue connessioni punto-punto

alle estremita occupate dai residui.

Il procedimento preso in esame in questo modello e del tutto simile e

verra illustrato dettagliatamente nel Capitolo 4.

38

Capitolo 4

Integratore per un pendolosferico

Siamo in presenza di una massa puntiforme soggetta a forze esterne

vincolata a muoversi sulla superficie di una sfera.

Come e noto, anche se il problema e integrabile analiticamente, la

soluzione fa ricorso a funzioni ellittiche spinose dal punto di vista

computazionale; e preferibile percio la ricerca di una integrazione nu-

merica piu maneggevole.

Trattiamo un punto vincolato a distanza fissa dall’origine alla stregua

di un corpo rigido unidimensionale (che chiameremo asta) la cui mas-

sa e concentrata su di una estremita. Individuare l’orientazione dell’

asta equivale a individuare la posizione del punto sulla sfera.

4.1 Dinamica

4.1.1 Pendolo sferico

Definiamo

e vettore (versore) che indidua l’asta (supposta di lunghezza unitaria).

~ω velocita angolare

~L momento angolare

Tutti i vettori sono definiti rispetto ad un riferimento esterno.

39

Figura 4.1: Pendolo sferico. n individua il piano istantaneo dirotazione, t la direzione tangente.

Come e stato anticipato nel Capitolo 3 l’asta (un qualsiasi oggetto

unidimensionale) rappresenta un corpo rigido degenere; il suo tensore

d’inerzia ha determinante nullo rispetto a qualsiasi polo e, come illus-

trato in Figura 4.1, i vettori ~L e ~ω sono paralleli e perpendicolari al

piano di rotazione istantanea dell’asta.

Notiamo che non viene tolta generalita al problema supponendo la

massa dell’asta unitaria e concentrata alla sua estremita.

Sotto tale ipotesi, con riferimento alla Figura 4.1, otteniamo.

~v ≡ ~e = ~ω × e = |ω|t (4.1)

Dalla definizione di momento angolare

~L ≡∑

i

mi · ~ri × ~vi (4.2)

dove la sommatoria e estesa alle particelle che compongono il corpo

rigido, si ha nel nostro caso

~L ≡ m · e× ~v =

m · |ω| · e× t =

|ω|n = ~ω

(4.3)

Il momento angolare risulta coincidere con la velocita angolare.

40

Questo risultato valido per oggetti unidimensionali semplifica notevol-

mente l’integrazione del sistema rispetto al caso generico del corpo

rigido.

4.1.2 Corpo rigido

Definiamo con A la matrice di rotazione che identifica la posizione

del corpo (degli assi principali d’inerzia) rispetto agli assi del sistema

esterno (una trattazione analoga puo essere fatta indifferentemente

per il quaternione).

La matrice A permette di trasportare qualsiasi vettore dal riferimento

esterno a quello solidale (accentato)

~ω ′ = A~ω e ~L′ = A~L (4.4)

ricordiamo inoltre che

~L = J ~ω (4.5)

dove J rappresenta la matrice d’inerzia riferita agli assi del laboratorio

legata alla matrice d’inerzia degli assi principali J′ dalla relazione

J = A†J′A (4.6)

Descriveremo la dinamica del corpo rigido dall’evoluzione delle grandezze

A e ~ω per le quali ricaviamo le variazioni temporali.

A = ~ω × A (4.7)

(il prodotto vettoriale e distribuito sulle colonne della matrice) o in

forma matriciale

A = Ω · A dove Ω =

0 ωz −ωy

−ωz 0 ωx

ωy −ωx 0

(4.8)

Dalla meccanica sappiamo che

~L = ~τ (4.9)

41

dove τ rappresenta il momento delle forze totale.

Accoppiando le equazioni (4.8) e (4.9) si ottiene il sistema di equazioni

di Eulero per il moto del corpo rigido.

Nel metodo di impronta leapfrog-like di Fincham [13] l’algoritmo e

strutturato

~L(t) = ~L(t− h2) + h

2~τ (t)

A(t+ h) = A(t) + hA(t+ h2)

(4.10)

Rimandando all’Appendice B per i dettagli, osserviamo nel sistema

(4.10):

• la seconda equazione deve essere risolta iterativamente poiche,

come si osserva in (4.7) e confermemente alle caratteristiche

dei moti rotazionali, la derivata temporale della posizione ha

dipendenza posizionale.

• essa deve inoltre poter garantire la rigidita del corpo conservando

la ortonormalita di A.

4.2 Algoritmo per il pendolo sferico

Il sistema (4.10) si semplifica notevolmente se si passa al caso degenere

dell’asta (o pendolo sferico). La relazione (4.3) premette di tradurre

immediatamente l’equazione per ~L in una equazione per ~ω e, paral-

lelamente ad un aggiornamento di tale vettore, si provvedera ad ogni

passo d’integrazione all’avanzamento del vettore e.

Il sistema da integrare e il seguente

˙e = ~ω × e

~ω = ~τ (~ω, e)

(4.11)

dove ~τ = e× ~f (e, ~w) rappresenta il momento totale delle forze agenti1.

1sono incluse anche le forze vincolari lungo l’asta qualora sia imperniata(pendolo) o connessa ad altre aste (catena).

42

L’integrazione sara portata ad un ordine di errore h3 e garan-

tira implicitamente la conservazione di ‖e(t)‖ tramite un avanzamento

compiuto con rotazioni sul piano ortogonale a ~ω.

4.2.1 Integratore al 2o Ordine

Procediamo con lo sviluppo di Taylor per un incremento h del vettore

e (per semplicita si omette il cappelletto)

e(t+ h) = e(t) + he(t) +h2

2e(t) + o(h3) =

e(t) + h(ω(t) × e(t)) +h2

2[ω × e+ ω × (ω × e)] + o(h3)

(4.12)

Se ω e costante, il moto rotazionale uniforme di e si integra facil-

mente nel modo (vedi Capitolo 2)

e(t) = qω,α(t) e0 qω,α(t) = q(ω,2α(t))e0 (4.13)

Alla seconda identita si giunge sfruttando l’ortogonalita di e e ~ω e

usando la forma ridotta (1.10).

Dato l’incremento h definiamo il quaternione che compie la rotazione

q(ω, h) = cosα + ω sinα con α = |ω|h (4.14)

L’equazione (4.13) si discretizza percio

e(t+ h) = q(ω, h)e(t) (4.15)

Proseguiamo l’analisi intoducendo l’esponenziale ottenuto dallo svilup-

po di seni e coseni e le parentesi di Poisson per l’operatore quater-

nionale che, come sappiamo, non commuta.

Si dimostra che

q(ω,h

2) e q(ω,

h

2) = exp

ωh

2e exp−ωh

2= exp [

ωh

2, · ]e (4.16)

Si ottiene che

e(t+ h) = exp [(h

2)ω(t+

h

2)] e(t) exp [−(

h

2)ω(t+

h

2)] + o(h3) (4.17)

43

e

ωτ

t t+ht+h/2 t t+ht+h/2 t t+ht+h/2

Figura 4.2: Evoluzione della posizione e, velocita angolare ω e mo-mento delle forze τ all’interno del passo di integrazione ordine 2. Lefrecce indicano il calcolo a partire da grandezze note (in grigio).

coincide con lo sviluppo di e al secondo ordine scritto in (4.12).

Per completare l’algoritmo e necessario aggiornare ω al tempo t + h

secondo uno schema leapfrog-like

ω(t+h) = ω(t)+hτ(t+h

2) dove τ(t+

h

2) = e(t+

h

2)×f(t+

h

2) (4.18)

L’algoritmo, utilizzando la forma ridotta 2 si compone percio di un

passo intermedio (midstep)

ω∗(t+ h2) = ω(t) + h

2τ(t)

e∗(t+ h2) = q(ω(t), h

2) e(t)

(4.19)

e di un passo completo (fullstep)

ω(t+ h) = ω(t) + hτ ∗(t+ h2) + o(h3)

e(t+ h) = q(ω∗(t+ h2), h) e(t) + o(h3)

(4.20)

Il passo di integrazione per l’ordine 2 e schematizzato in Figura 4.2.

4.2.2 Integratore al 3o Ordine

Per passare ad ordini superiori si espande lo sviluppo (4.12) in

e(t+ h) = e(t) + he(t) +h2

2e(t) +

h3

6

...e (t) + o(h4) (4.21)

2Con la discretizzazione dell’evoluzione temporale si perde in linea di principiol’ortogonalita tra e(t) e ω(t + h

2) che non sono piu simultanei. L’errore che si

commette risulta comunque largamente inferiore alla precisione raggiunta.

44

e si suddivide il passo h in opportuni avanzamenti parziali per ω e

rotazioni concatenate per e che mantengono la unitarieta della trasfor-

mazione.

Il passo completo e dato da

e(t+ h) = q[ω(t) + h3τ ∗(t+ h

2), h

3] q[ω(t) + h

4τ(t), 2h

3] e(t) + o(h4)

ω(t+ h) = ω(t) + 3h4τ ∗(t+ 2h

3) + h

4τ(t) + o(h4)

(4.22)

dove il termine τ ∗(t + h2) deve essere fornito all’ordine o(h2) tramite

midstep mentre il termine τ ∗(t+ 2h3

) deve essere fornito all’ordine o(h3)

tramite un fullstep relativo al passo 2h3

.

L’evoluzione e(t) −→ e(t+ h) al terzo ordine comporta un totale di 4

cicli.

4.3 Catena di pendoli sferici

Il moto di una catena composta da N pendoli e integrabile come il

moto di N pendoli indipendenti a patto di conoscere le forze vincolari

che si trasmettono lungo la stessa.

Tali forze si estraggono risovendo un sistema che imponga che il mo-

to relativo di due masse contigue sia puramente rotazionale (velocita

tangente all’asta che le connette e accelerazione solo centripeta).

Le forze vincolari possono trasmettersi solo lungo le aste,

~φi = φiei

L’accelerazione della massa mi vale percio

~ai = ~F toti = ~fi + φiei − φi−1~ei−1 (4.23)

dove ~fi sono le forze sterne e ~F toti le forze totali.

Affinche l’accelerazione relativa sia solo centripeta deve valere

ei · (~ai+1 − ~ai) = − (~vi+1 − ~vi)2 = − ~ωi

2 (4.24)

45

Inserendo (4.23) in (4.24) si ricava il sistema lineare tridiagonale

φi−1(eiei−1) − 2φi + φi+1(eiei+1) = − ei(~fi+1 − ~fi) − ~ωi2 (4.25)

dove il termine a destra e noto.

Il sistema 4.25 permette di ottenere i coefficienti delle forze vincolari

φi e conseguentemente le forze totali F toti (esterne + vincolari) agenti

su ogni massa mi.

I momenti delle forze τi da introdurre nei sistemi 4.20 (2o ordine) o

4.22 (3o ordine) sono ricavati

τi = ei × (F toti+1 − F tot

i ) (4.26)

Il calcolo dei momenti τi, essendo essi una funzione di ei e ωi, va

compiuto separatamente all’interno di un ciclo di integrazione ogni

qual volta sia necessario, ad esempio

τ(t+h

2) = τ [e(t+

h

2), ω(t+

h

2)].

Data percio una catena di N aste, si integrano indipendentemente

i moti relativi di ogni massa mi rispetto alla precedente mi−1

ei = ωi × ei

ωi = τi (ω1...N , e1...N)

(4.27)

la posizione assoluta del punto i-esimo al tempo t e data da

~Ri(t) =i−1∑

j=0

ei(t)

e la sua velocita

~vi(t) = ~vi−1(t) + ~ωi(t) × ei(t)

Ricordiamo che l’ordine dell’errore di integrazione e riferito ai vet-

tori locali e e non alle posizioni globali ~R.

Lo stesso vale per le velocita locali (angolari) ~ω e globali ~v.

46

0 20 40 60 80 100tempo

10-9

10-8

10-7

h = 0.0005

h = 0.002

h = 0.004

h = 0.001

EE

∆

Figura 4.3: Deriva energetica per l’integraore al 3o ordine nel casodel pendolo semplice (N=1) con differenti passi di integrazione h. Inordinata in scala logaritmica lo scarto energetico relativo.

4.4 Test al calcolatore

Per la programmazione si e fatto uso del linguaggio C + + che dispone

dell’utilizzo di classi sulle quali strutturare l’algebra voluta (compo-

sizione di quaternioni e vettori).

Sono stati compiuti ripetuti test sull’integratore per verificare e con-

frontare la stabilita al variare del passo h (Figura 4.3), dello schema

di integrazione (Figura 4.5) e del numero di aste (Figura 4.4 e 4.6).

E stato sperimentato in aggiunta un algoritmo che prevede la modu-

lazione del passo di integrazione in relazione al crescere della cinetica

del sistema (come riferimento si prende max(~ωi) ad ogni passo) in mo-

do da mantenere la raffinatezza di integrazione consona all’ incremento

delle velocita delle aste (vedi Figura 4.7).

4.4.1 Risultati

Gli esperimenti di ripiegamento della proteina (Capitolo 6) prevedono

un’ analisi statistica fatta su di un gran numero di simulazioni (meto-

do MonteCarlo).

La scelta dello schema e del passo di integrazione e percio fatta anche

47

50 60 70 80 90 100tempo

0

2×10-5

4×10-5

6×10-5

EE

∆

N = 4

N = 3

N = 1

N = 2

Figura 4.4: Deriva energetica dell’integratore ordine 3 con h = 0.005per una catena di N pendoli. In ordinata lo scarto energetico relativo.

1,0e-041,0e-03

h

1e-12

1e-10

1e-08

1e-06

∆EE

Figura 4.5: Deriva energetica per gli integratori ordine 2 (linea a punti)e 3 (linea tratto-punto) in relazione alla riduzione del passo h. Inordinata l’errore relativo massimo a t = 100.

48

0 2 4 6 8 10N

1e-06

1e-05

1e-04

EE

∆

Figura 4.6: Deriva energetica dell’integratore ordine 3 con h = 0.005 alvariare del numero N di aste. In ordinata lo scarto energetico relativomassimo a t = 100 in scala logaritmica.

0 10 20 30 40tempo

ordine 2

ordine 2

ordine 3∆EE

-5 10-5.

.5 10. -5

-1 10-4.

h modulato

0

Figura 4.7: Deriva energetica per una catena di 50 pendoli con differ-enti schemi di integrazione. In nero lo schema con passo h modulatoal variare di ωmax. Nella configurazione iniziale la catena e orizzon-tale, acquistando velocita si ha (t = 30s) lo shock cinetico del ”colpodi frusta“ sull’ultimo elemento della catena. L’integratore ordine 3reagisce meglio dell’ ordine 2 con h modulato. In ordinata lo scartorelativo.

49

in relazione alla velocita di calcolo.

Il tempo macchina 3 richiesto per l’evoluzione del moto di una cate-

na di 10 pendoli fino a t = 100 per l’integratore ordine 3 con passo

h = 0.005 e di circa 5 secondi.

Vale a dire che il tempo necessario per un singolo passo e di circa

2.5 · 10−4 secondi.

L’integratore al terzo ordine compie quattro cicli interni per passo

mentre al secondo ordine solamente due (full-step e mid-step). Percio

nel caso del pendolo (campo costante) il tempo di calcolo cresce lin-

earmente con il numero di aste coinvolte raddoppiandosi se si passa

dal secondo al terzo ordine.

Nel modello della proteina verranno introdotte interazioni a coppie; di

conseguenza la misura delle distanze reciproche fa sı che il tempo di

calcolo cresca come ∼ N2. Per economizzare i processi, dato il veloce

decadimento a zero del potenziale a lungo range, verra eseguito un

cut-off compatibile con la dinamica.

Nel prossimo Capitolo verra descritta inoltre la procedura di inte-

grazione adottata in presenza di contributi stocastici che simulino la

presenza di un bagno termico.

3processore PentiumIII 866MHz

50

Capitolo 5

Modello per la proteina

5.1 Introduzione

La difficolta nella costruzione di un modello per la proteina sta nel

definire un set di potenziali che colga gli aspetti salienti della dinami-

ca con l’introduzione del minor numero di parametri.

Gran parte dei modelli off-lattice, ossia al continuo, interpretano le

interazioni come forze che intervengono tra coppie di amminoacidi.

La proteina e schematizzata come una catena di punti identificabili

con i residui Cα dove si considera concentrata la massa dell’ammi-

noacido (considerata costante e unitaria).

Le catene laterali R non sono prese in considerazione; quello che ri-

mane e percio la spina dorsale (backbone) del filamento formata dai

legami peptidici.

Una ulteriore semplificazione si attua restringendo la varieta delle clas-

si di apparteneza degli amminoacidi coinvolti rispetto alle 20 presenti

nel mondo organico.

La trattazione piu semplice e probabilmente quella fatta nel modello

H-P dove gli amminoacidi sono divisi in due tipologie in base alla loro

tendenza idrofobica (H) o polare (P).

Stillinger su questo schema ha costruito un modello minimale (toy-

model) dove i residui si muovono in uno spazio bidimensionale [23].

Il modello Thirumalai-Honeycutt [18], al quale fa riferimento questo

lavoro, coinvolge tre classi di amminoacidi: idrofobici (B), idrofilici (L)

51

Figura 5.1: Rappresentazione del filamento di proteasi del Virus HIV.Immagine tratta dal sito web del dipartimento di Chimica, Universitadel Wisconsin.

e neutri (N); le sequenze sono espressa percio in codici di tre lettere.

5.2 Potenziali

Nella Sezione 1.3 sono state elencate le principali forze che agiscono

nella fisica della proteina e che ora vanno tradotte in potenziali per il

modello.

Il legame peptidico presente tra residui contigui e gia rappresentato

dalle forze vincolari del modello a catena che conserva le distanze.

Nei modelli ad approccio atomico tale legame e simulato con poten-

ziali armonici quadratici e quartici ad alta frequenza.

Schematizziamo i restanti potenziali ricordando che chiarificazioni riguar-

do a dimensionalita e unita di misura verranno fornite nel prossimo

Capitolo (Sezione 6.2.2). Per chiarezza riportiamo in anticipo i valori

delle unita riportate sui grafici che seguono; l’energia (i potenziali V )

e espressa in unita energetiche ε ≃ 9600J/mol; le distanze sono in

52

unita metriche a ≃ 4A corrispondente alla distanza tra elementi del

filamento; il tempo e in unita τ ≃ 1ps.

Angolo di legame Tra due segmenti contigui di catena, ossia tra i

residui i - i+2, poniamo una molla armonica ideale che a riposo forma

un angolo θ0 = 105o fra i segmenti.

Vθ =1

2Kθ(d− d0)

2 con d0 =√

2 − 2 cos θ0 (5.1)

La costante Kθ = 20ε/rad2 dove ε rappresenta l’unita energetica

del sistema. Questa interzione e quantitativamente la piu intensa

ed e temporalmente la prima ad intervenire nel processo di folding

inducendo la formazione locale di “strisce a zig-zag”.

Angolo diedrale Tre segmenti contigui individuano due piani fra i

quali definiamo l’ angolo diedrale φ il cui potenziale

Vφ = A(1 + cosφ) +B(1 + cos 3φ) (5.2)

ha un minimo principale (trans) e due minimi transienti (gauche)(vedi

Figura 5.2).

L’angolo φ rappresenta la rotazione attorno all’asse del legame pep-

tidico che connette due amminoacidi.

Le costanti sono fissate A = B = 1.2ε eccetto nel caso in cui siano

presenti due o piu residui neutri nel qual caso A = 0, B = 0.2ε;

questo induce la formazione di tornanti (turn regions) che separano

due β-sheet dando luogo alle tipica struttura a forcina (vedi prossimo

Capitolo).

Potenziale Lennard-Jones Per esprimere le interazioni a lungo

range viene utilizzato un potenziale che schematizza una forte repul-

sione idrofobica a stretta vicinanza e una debole attrazione che decade

rapidamente con la distanza.

Il potenziale Lennard-Jones caratterizza di per se una comune tipolo-

gia di modelli molecolari interpretando in modo soddisfacente le tipiche

53

0 90 180 270 360φ

V

minimi transientiminimo principale

Figura 5.2: Potenziale relativo all’angolo diedrale; la linea tratteggiatacorrisponde al potenziale in presenza di due o piu residui N. (la scalaenergetica non e riportata).

peculiarita delle interazioni a coppie 1.

Sappiamo dalla fisica organica che le semplici interazioni a due corpi

sono limitanti; del resto, al livello descrittivo di questo modello, esse

danno una descrizione fisica ragionevolmente buona; le interazioni a

tre o piu corpi, inoltre, sono raramente utilizzate in Dinamica Moleco-

lare per il loro alto costo in termini di tempo di calcolo.

Il potenziale Lennard-Jones e correntemente definito

VLJ = 4εS1

[

1

r12− S2

1

r6

]

(5.3)

dove le costanti introdotte ponderano l’ omofilia delle tre classi:

• S1 = S2 = 1 per l’interazione B-B

• S1 = 2/3, S2 = −1 per le interazioni L-L e L-B

• S1 = 1 , S2 = 0 per tutte le interazioni che coinvolgono residui

N

1Gli studi (Bobetic Baker Maitland Smith, 1971) effettuati su una grande basedi dati sperimentali (scattering, spettroscopia, proprieta dello stato solido ...) han-no permesso di estrarre la curva del tipico potenziale a coppie tra elementi atomici(Potenziale BBMS studiato per atomi di Argon).

54

0,8 1 1,2 1,4r

0

100

200

300

V

N - X

L - LB - B

Figura 5.3: Potenziale Lennard-Jones. Idrofobi (B), Idrofili (L), Neu-tri (N), generico (X). In ordinata la grandezza e in unita energeticheε, in ascissa in unita metriche a.

L’interazione B-B e l’unica delle tre a possedere una componente at-

trattiva che, come si vede in Figura 5.3, forma un minimo di potenziale

a r ≃ 0.9.

Questo potenziale e stato modificato con un cut-off smussato a r = 5

per agevolare i processi di calcolo.

In Figura 5.4(b) si osservano i valori dei tre potenziali suddetti nel

corso di una simulazione. Il potenziale Lennard-Jones, che come si puo

notare fornisce il contributo energetico inferiore, e in realta decisivo

nella formazione della struttura terziaria del filamento.

Come e stato discusso nel Capitolo 4, l’integratore scelto per le sim-

ulazioni del modello per la proteina e quello ordine 3 con passo h =

0.005.

5.3 Bagno termico

Gli esperimenti sul ripiegamento della proteina sono condotti in pre-

senza di solventi, perlopiu soluzioni acquose a temperatura ambiente

e pH neutro.

Il mezzo polare induce meccanismi di interazione tra coppie di am-

55

0 20 40 60 80 100

0 20 40 60 80 100tempo

0

20

40

60

(a)

(b)

E

V

3 10

6 10

0

-5.

. -5

E∆

Figura 5.4: (a) Deriva energetica dell’integratore ordine 3 per una pro-teina composta da 22 residui (N=21) in assenza di bagno termico (en-ergia costante). In ordinata lo scarto relativo. (b) Potenziali dell’an-golo di legame (nero), angolo diedrale (blu), Lennard-Jones (rosso).Passo h = 0.005. In ordinata l’energia in unita ε. In entrambi igrafici il tempo in unita τ .

minoacidi che sono gia state inclusi nel modello tramite il potenziale

Lennard-Jones.

Resta da considerare il contributo termodinamico del solvente che

schematizzi gli urti tra molecole d’acqua ed elementi del filamento.

In gran parte delle simulazioni di dinamica molecolare il bagno termi-

co viene schematizzato con la cosiddetta dinamica di Langevin per il

limite di bassa frizione o dinamica browniana [25].

Utilizzeremo per la catene di residui l’analoga procedura che si applica

per insiemi di molecole libere (non soggette a vincoli rigidi) dopo aver

apportato alcune modifiche ed averne verificato la validita.

5.3.1 Dinamica di Langevin

Su di una particella immersa in un liquido agiscono una forza di

frizione (damping) e una forza dipendente da una variabile stocas-

tica (rumore).

La rappresentazione piu semplice nel limite di bassa frizione e espressa

56

(caso unidimensionale)

x = v

v = f(x, v, t) − γ v + ξ

(5.4)

dove f e il campo di forze noto, γ e il parametro d’attrito nella formu-

lazione di Stokes, ξ e una variabile stocastica con una distribuzione a

media nulla

〈ξ(t)〉 = 0 (5.5)

e funzione di correlazione proporzionale alla funzione δ.

〈ξ(t)ξ(t′)〉 = q2δ(t− t′) (5.6)

ξ e detta anche forza di Langevin o, date le condizioni 5.5 e 5.6, ru-

more bianco.

Il sistema 5.4 si traduce nella equazione di Fokker-Planck per la

funzione di distribuzione ρ(v, t) che indica la probabilita di trovare

una velocita compresa nell’intervallo (v, v + dv).

∂ρ

∂t= γ

∂(vρ)

∂v+ 2q2∂

2ρ

∂v2(5.7)

La soluzione dell’eq. 5.7 nel caso stazionario ∂ρ∂t

= 0 ci da la dis-

tribuzione di velocita ρ(v) all’equilibrio termico.

La media quadratica della velocita vale

〈v2〉 =q2

2γ

e per il teorema di equipartizione

〈E〉 =1

2m〈v2〉 =

1

2mq2

2γ=

1

2kT (5.8)

Questo risultato, estendibile al caso tridimensionale, mette in relazione

i parametri γ e q con la temperatura T (la costante di Boltzmann k e

assorbita nella grandezza T ).

T =q2

2γ(5.9)

Ancora una volta, dettagli sulle dimensionalita e le unita di misura

della simulazione verrano esposti nel prossimo Capitolo.

57

5.3.2 Dinamica di Langevin all’interno dell’algo-ritmo

La dinamica del sistema e quindi composta da una parte determinis-

tica (Lagrangiana) e da una parte stocastica rappresentata dagli urti

casuali (browniani).

Per integrare quest’ultima si traduce l’impatto casuale delle molecole

del solvente in un impulso (kick) che modifica istantaneamente le vari-

abili dinamiche ~ωi al termine di ogni passo.

Le forze stocastiche ~ξ dell’eq. 5.4 (esteso al caso 3D) sono generate

alla fine del ciclo di integrazione Lagrangiano descritto nella Sezione

4.2.

Per la compatibilita vincolare le forze ~ξ sono trasformate in momen-

ti delle forze con una procedura simile a quella descritta in 4.3 con

la differenza che sono escluse dall’eq. 4.25 le componenti inerziali ~ω2i

poiche esse non svolgono funzione in una interazione istantanea.

Il sistema tridiagonale si traduce in

φi−1(eiei−1) − 2φi + φi+1(eiei+1) = − ei(~fi+1 − ~fi) (5.10)

e da qui, in analogia con 4.25, estraiamo i momenti delle forze ~τR.

Le velocita angolari ωi vengono cosı aggiornate con la sostituzione

~ω −→ ~ω +√h~τR (5.11)

dove i vettori ~τR sono i momenti random la cui interazione, seguendo

l’analogia con la dinamica browniana, e pesata dal coefficiente√h.

5.3.3 Test termodinamico

Per validare l’interpretazione dinamica che si e data al bagno termico

e verificare la consistenza della simulazione dal punto di vista termod-

inamico e statistico eseguiamo una prova di moto browniano.

Dopo aver ”spento” tutti i potenziali (angolare, diedrale, Lennard-

Jones), un filamento di polimero (N=21) e lasciato evolvere sotto

l’azione degli urti e della frizione del solvente.

58

0 0,1 0,2 0,3 0,4 0,5T0

5

10

15

20

<K>

0 0,5 1 1,5

||v||

0

0,5

1

1,5

2

ρ

(a)

(b)

Figura 5.5: (a) Valori dell’energia cinetica media per un filamento diN=21 aste all’equilibrio termodinamico sotto la unica azione del bagnotermico (urti+frizione); in grassetto la retta 〈K(T )〉 = 2N+3

2kb T . K e

T sono espresse in unita energetiche ε. (b) Relativa distribuzione dellevelocita (modulo) di un singolo elemento del filamento; in grassetto la

curva ρ(v) = 4πv2

(3πT )3/2exp

(

− v2

3T

)

. In ascissa la velocita e espressa in

a · τ−1.

59

Dopo un breve intervallo di termalizzazione, il sistema raggiunge l’e-

quilibrio con l’ambiente con una energia che oscilla intorno al valore

medio 〈K(T )〉.Come si vede in Figura 5.5 l’andamento di 〈K(T )〉 e ben descritto

dalla retta

〈K(T )〉 =Ng

2kb T =

2N + 3

2kb T (5.12)

dove Ng rappresenta il numero di gradi di liberta della catena (2 gradi

di rotazione per ogni asta piu 3 gradi di traslazione per il centro di

massa) in accordo con la legge di equipartizione per sistemi in equi-

librio.

Una ulteriore verifica si ha sulla distribuzione maxwelliana delle ve-

locita per un singolo elemento della catena. Come si osserva in Figura

5.5, l’istogramma dei dati raccolti riproduce la curva

ρ(v) =4πv2

(3πT )3/2exp

(

− v2

3T

)

(5.13)

di una particella che, in presenza di vincolo, possiede 2 gradi di liberta.

5.4 Similarita χ

Il problema del ripiegamento della proteina esamina il nesso esistente

tra la sequenza di residui assegnata (struttura primaria) e il raggiung-

imento di una configurazione tridimensionale di equilibrio (struttura

terziaria).

Per questa analisi e utile costruire una funzione che valuti il grado di

similarita tra la configurazione raggiunta e lo stato nativo relativo alla

sequenza.

Affinche la similarita sia invariante per rotazione, le variabili in gioco

sono le distanze inter-residuali rij .

Definiamo la distanza configurazionale [28]:

χ = 1 − 2

2 + n(n− 3)

n−2∑

i=1

n∑

j=i+2

Θ(ǫ− |rij − rN

ij|) (5.14)

dove n e il numero di residui, rN

ij sono le distanze inter-residuali dello

stato nativo, Θ e la funzione a gradino di Heaviside e il parametro ǫ

60

0 20 40 60 80 100tempo

0

0,2

0,4

0,6

0,8

1

χ

Figura 5.6: χ(t) durante un processo di folding per un filamento di21 elementi a T=0.1. A t ≃ 70 viene stabilmente raggiunta la con-figurazione di stato nativo. Il tempo in ascissa e espresso in unitatemporali τ .

(fissato a 0.2) fa si che si trascurino piccole oscillazioni intorno allo

stato nativo sotto l’azione del bagno termico.

Il fattore 22+n(n−3)

normalizza il conteggio delle sommatorie sul numero

di coppie esaminate (sono escluse le coppie di residui contigui), cos-

icche χ = 0 corrisponde ad una struttura perfettamente coincidente

mentre per χ & 0.8 si ha una struttura del tutto dissimile.

In Figura 5.6 si osserva l’evoluzione temporale di χ durante un pro-

cesso di folding.

61

62

Capitolo 6

Simulazioni e risultati

6.1 Ipotesi teoriche sul protein folding

Gli esperimenti compiuti sul ripiegamento della proteina sono condotti

attraverso la denaturazione della proteina (alterazioni dell’ambiente,

tipicamente innalazamento della temperatura) e studiano il ripristino

spontaneo della struttura terziaria e delle relative attivita enzimatiche

[2].

Cio conduce all’ipotesi che il ripiegamento del filamento polimerico sia

un processo termodinamico particolarmente delicato e che per la sua

completa comprensione sia necessaria una visione globale del panora-

ma energetico del sistema.

La cosiddetta ipotesi termodinamica afferma che la struttura tridi-

mensionale a cui si giunge al termine del ripiegamento (stato nativo) e

quella corrispondente al minimo assoluto dell’energia libera di Gibbs.

Si delinea cosı la visione di un profilo energetico a imbuto (funnel land-

scape) in cui ogni cammino dinamico conduce rapidamente al minimo

suddetto senza che il sistema rimanga intrappolato negli stati metasta-

bili dei minimi locali.

L’approccio cinematico delle simulazioni , d’altra parte, suggerisce la

presenza di barriere energetiche che ostacolano il percorso del sistema

nella ricerca dello stato nativo.

A riguardo, Levinthal [29] stabilisce su base teorica che il tempo per

la ricerca casuale cinematica dello stato di minima energia libera da

63

parte di filamenti anche di moderate dimensioni cresce smisuratamente

e incompatibilmente con i tempi biologici (ordine che va dal millisec-

ondo al secondo).

Per risolvere questo paradosso l’intricato e multidimensionale panora-

ma energetico dovrebbe essere caratterizzato da un numero ridottis-

simo di cammini preferenziali (pathways) che incanalano il sistema

direttamente verso la stato nativo evitando una lunga ricerca casuale.

Thirumalai [18] avanza una seconda ipotesi sulla base dell’approccio

cinematico.

L’ ipotesi metastabile presuppone la presenza di una regione lo-

calizzata di minimi energeticamente distinti, ma configurazionalmente

prossimi allo stato nativo.

Tali minimi metastabili sono confinati da tenui barriere di potenziali

che non precludono occasionali transizioni all’interno della regione.

In questo modo il sistema, una volta concluso il ripiegamento, evolve

visitando statisticamente minimi energeticamente differenti, ma che

conservano un forte carattere di similarita con lo stato nativo.

Cio equivale a dire che da uno stato denaturato (non ripiegato, random

coil) esistono molteplici cammini che conducono a una zona metasta-

bile pseudo-nativa.

6.2 Simulazione

Come e stato detto nel Capitolo 1, per un approccio cinematicamente

controllato del protein folding sono state adottate diverse metodologie

di simulazione mirate a gettare luce sui diversi aspetti del problema.

Il modello tridimensionale off-lattice di Honeycutt & Thirumalai, ri-

costruito in questa tesi, e stato utilizzato per svariati tipi di ricerca

sulle correlazioni tra i parametri di ingresso nella simulazione (tipolo-

gia e coefficienti delle interazioni, metodo di codifica e assegnazione

della sequenza di amminoacidi, procedura termodinamica e coeffieci-

enti per il bagno termico, etc..) e valori ricavati all’uscita (configu-

64

razione, frequenza e stabilita dei minimi visitati, tempi di rilassamen-

to, similarita χ con lo stato nativo, panorami energetici, etc..).

Con l’analisi statistica elaborata su di un gran numero di simulazioni

indipendenti (metodo MonteCarlo) e possibile cosı evidenziare le in-

fluenze che hanno sul processo, ad esempio, le caratteristiche del sol-

vente [35], le mutazioni nella sequenza primaria [20] o la parametriz-

zazione dei potenziali [33] [34] e avanzare ipotesi sulla dinamica sog-

giacente: carattere diffusivo del moto [32], contrapposizione (frus-

trazione) tra ordine locale e globale della configurazione [30] [31], dis-

tinzione fra sequenze relative a buoni o cattivi ripegatori.

In questo lavoro si e focalizzata l’attenzione sulla ricerca

dei minimi dell’energia meccanica del sistema per una se-

quenza associata ad un presumibile buon ripiegatore tramite

un metodo MonteCarlo con simulazioni in un bagno termico

a temperatura costante.

Le posizioni di equilibrio raggiunte a seguito del ripiegamento sono

state archiviate in registri di variabili a cui puo avere accesso un pro-

gramma di visualizzazione della struttura (vedi Figura 6.1).

Tramite un’ analisi con ripetute simulazioni indipendenti, ad ogni min-

imo e stata associata una probabilita o popolazione; in una ideale rap-

presentazione del panorama energetico ad un minimo molto popolato

dovrebbe corrispondere un bacino di attrazione ampio a cui converge

un gran numero di cammini e una buca di potenziale profonda che

imprigioni stabilmente il sistema.

6.2.1 Sequenza

Studi condotti sull’accessibilita cinetica dei minimi e sulla loro stabilita

[36] mostrano come queste si rivelino ben correlabili al parametro σ

65

Figura 6.1: Configurazioni di equilibrio raggiunte al termine dellasimulazione

definito

σ =Tθ − TF

Tθ

0 < σ < 1 (6.1)

dove Tθ e TF rappresentano la temperatura del bagno sotto alla quale

statisticamente avvengono rispettivamente il collasso e il ripiegamen-

to.

A piccoli valori di σ corrispondono ripiegatori veloci con un bacino

d’attrazione dominante (Native Basin of Attraction) ossia buoni rip-

iegatori, mentre ad alti valori di σ corrispondono sequenze con molte

trappole cinetiche (minimi locali o Competing Basin of Attraction).

La sequenza utilizzata principalmente in queste simulazioni, denom-

inata Beta21, e costituita da 22 elementi residuali (21 aste che rap-

presentano i legami peptidici tra gli amminoacidi) il cui codice a tre

lettere e riportato in Figura 6.2.

Lo stato nativo di Beta21 e la caratteristica β-sheet che rappresenta il

ricorrente elemento topologico della struttura secondaria (vedi Capito-

lo 1). I due residui neutri al centro della sequenza favoriscono, tramite

il potenziale diedrale (5.2), la formazione del cardine nella struttura

a forcina; l’attrazione idrofobica induce le due meta del filamento a

porsi parallelamente.

Per Beta21 si ha σ = 0.14, indice di un buon ripiegatore [36].

A titolo di confronto e stato esaminato anche il tipico prototipo di

cattivo ripiegatore costituito da una catena completamente idrofobica

(Figura 6.2 c ).

66

Figura 6.2: (a) La configurazione dello stato nativoper Beta21 e una β-sheet. Il codice della sequenza eLBNBBBLBBBNNBBLBLBBBLB = (LB)N(B)3L(B)3(N)2B(BL)2(B)3(LB)dove alle lettere B, L, N corrispondono residui idrofobici (rosso),idrofilici (blu) e neutri (bianco). (b) La configurazione dello stato na-tivo per Beta45 e una β-sheet tripla. (c) La sequenza completamenteidrofobica (B)21 e un esempio di cattivo ripiegatore.

67

6.2.2 Dimensionalita

Nello sviluppo del modello e nel corso della simulazione, tutte le

grandezze sono state espresse in unita di misura ridotte.

Con semplici riscalature queste possono essere ricondotte a tre unita

di base [37]: massa M , lunghezza a ed energia ε .

M rappresenta la massa molecolare media per un amminoacido (vedi

Tabella 1.3) uguale a 110 g/mol; a e la distanza inter-residuale media

uguale a 3.8 A; ε = 2.3 kcal/mol = 9628 J/mol e l’unita di misura

energetica con la quale sono espressi i valori sui grafici.

Consistentemente l’unita di tempo ridotta τ =√

ε/Ma2 corrisponde

approssimativamente ad 1 ps.

La temperatura, assorbendo la costante k di Boltzmann, e espressa in

unita di energia T ∗ = kT/ε .

I tempi di folding ottenuti per le simulazioni isoterme risultano di gran

lunga piu brevi rispetto a quelli reali (& ms). Cio e dovuto in parte ai

valori di viscosita del solvente (vedi prossimo Paragrafo) e per il resto

alle macroscopiche semplificazioni adottate nel modello.

6.2.3 Modalita della simulazione

Sono state effettuate diverse serie di simulazioni indipendenti per dif-

ferenti temperature mantenute costanti.

Le configurazione iniziale e generata associando ad ogni vettore di

partenza ei(t = 0) una direzione casuale nell’ottante (+,+,+) in mo-

do che essa riproduca un filamento steso.

Si osserva che nelle condizioni di temperatura scelte, la proteina Beta21

impiega circa 100τ per termalizzare; l’istantanea della configurazione

di equilibrio raggiunta e presa a t = 200τ insieme al valore del minimo

energetico. Il coefficiente di viscosita del solvente e fissato al valore

γ = 0.05τ−1 che, tradotto in unita reali, da una frizione circa 10 volte

piu intensa di quella dell’acqua a 25oC [22]; cio riduce i tempi di calco-

lo, ma non altera sostanzialmente il risultato dell’analisi del processo

di ripiegamento.

68

Con le considerazioni fatte nel Capitolo 4, per l’integrazione della di-

namica si adotta l’algoritmo del 3o ordine con un passo di integrazione

h = 0.005τ . Il tempo reale di calcolo per una simulazione e dell’ordine

del minuto1.

Sono state processate per la sequenza Beta21 tre serie da 1000 sim-

ulazioni, rispettivamente alla temperatura T = 0.05, T = 0.10 e

T = 0.15; quest’ultima corrisponde in unita reali all’incirca alla tem-

peratura ambiente.

6.3 Analisi

6.3.1 Classificazione dei minimi

Ogni simulazione fornisce una istantanea della configurazione tridi-

mensionale a t = 200τ e il corrispondente valore del potenziale.

Un algoritmo procede alla classificazione dei minimi in base alla loro

distanza energetica e alla similarita χ.

Due configurazioni C1 C2 sono considerate identiche se la similarita

χ (C1, C2) 6 0.05 e se contemporaneamente i potenziali differiscono

per meno di uno scarto energetico ∆V ≃ ∆K ≃ 〈K〉 = 2N+32

kb T ≃21kb T dove T e la temperatura del bagno.

Quest’ultima stima e basata sulla comparabilita numerica esistente fra

le oscillazioni del potenziale V e dell’energia cinetica K per un sistema

in equilibrio oscillante e sulla relazione di equipartizione riportata in

5.12. Le oscillazioni V (t) per un filamento in equilibrio a diverse tem-

perature sono graficate in Figura 6.4.

Naturalmente anche l’ampiezza delle oscillazioni χ(t) variano con T

(Figura 6.3). Il valore costante dello scarto di discriminazione ∆χ =

0.05 e stato scelto in modo non critico considerando che due configu-

raioni all’interno di questo range di similarita hanno la stessa confor-

mazione tridimensionale, nella fattispecie, gli stessi valori di equilibrio

degli angoli diedrali.

1processore PentiumIII 866 MHz

69

0 50 100 150 200 250tempo

0

0,2

0,4

0,6

0,8

1

χ

zona metastabile

Figura 6.3: Fluttuazioni di χ rispetto alla configurazione finale (inquesto caso t = 250τ) durante simulazioni a T = 0.05 (blu), T = 0.15(rosso), T = 0.3 (viola); per T = 0.15 la χ oscilla all’interno dellazona metastabile delimitata dalle linee tratteggiate; per T = 0.3 latemperatura e troppo alta per consentire un corretto ripiegamento, laproteina si trova in uno stato random coil (filamento caotico).

50 100 150 200 250tempo

-20

-15

-10

-5

0

5

V

fluttuazioni di V

Figura 6.4: Fluttuazioni di V nelle stesse condizioni del grafico inFigura 6.3. T = 0.05 (blu), T = 0.15 (rosso), T = 0.3 (viola). Nelriquadro sono riportati i segmenti di scarto ∆V che verranno utilizzaticome discriminanti per la classificazione dei minimi.

70

6.3.2 Grafici Energia-Popolazione

Con tale procedura di classificazione dei minimi configurazionali rac-

colti si ottiene cosı uno spettro delle popolazioni per diverse temper-

ature del bagno.

Il grafico in Figura 6.5 e caratterizzato dalla presenza di un picco in

corrispondenza del minimo assoluto dell’ energia meccanica del sis-

tema. La relativa configurazione tridimensioanle visibile in Figura 6.2

corrisponde, come ci si poteva aspettare, alla struttura nativa a forci-

na della sequenza Beta21.

L’aumento della temperatura da T = 0.05 a T = 0.15 comporta una

sostanziale crescita della popolazione dello stato nativo.

Le oscillazioni dell’energia cinetica all’equilibrio termico per una sim-

ulazione a temperatura maggiore crescono in accordo con l’Eq. 5.12

e facilitano lo scavalcamento delle barriere di potenziale che intrap-

polano il sistema nei minimi locali.

A T = 0.15 notiamo inoltre la comparsa di minimi locali energetica-

mente prossimi allo stato nativo; per un analisi piu approfondita intro-

duciamo nello spettro la distanza configurazionale (similarita) definita

nel Capitolo 5.

6.3.3 Grafici Energia-Similarita-Popolazione

Stabilita l’esistenza di un minimo privilegiato che denominiamo stato

nativo, possiamo estrarre per i minimi locali la loro distanza configu-

razionale o similarita (vedi 5.4) ottenendo cosı una panoramica ”bidi-

mensionale” dello spettro popolativo di Beta21 (Figure 6.7 e 6.6).

Come si puo notare, i minimi enegeticamente prossimi allo stato na-

tivo, evidenziati anche in Figura 6.5, sono allo stesso tempo configu-

razionalmente simili.

Possiami definire minimi metastabili quelli la cui similarita o distanza

configurazionale con lo stato nativo e compresa tra 0.05 (sotto al quale

la similarita e considerata identita) e 0.1 (arbitrario).

71

-20 -15 -10 -5 0

V0

5

10

15

20

25

30

% Stato nativo

minimi metastabili a T=0.15

Figura 6.5: Spettro delle popolazioni per i minimi configurazionali dellasequenza Beta21 per T = 0.05 (blu), T = 0.1 (arancio), T = 0.15(rosso). In ascissa il valore energetico dell’equilibrio, in ordinata lapopolazione del minimo (con il criterio di classificazione riportato in6.3.1) su un totale di 1000 simulazioni. Si evidenzia la comparsa diminimi locali metasatabili (scarto di similarita χ ⊂ [0.05, 0.1], vedi6.3.3) in prossimita dello stato nativo per T = 0.15. Le configurazionidi stato nativo per le tre temperature coincidono (χ ≃ 0); lo scartoenergetico e dovuto al contributo termico medio dato al potenziale.

72

-0,2 0 0,2 0,4 0,6 0,8

χ

-25

-20

-15

-10

-5

0

V

Stato Nativo

T=0.05

zona metastabile

Figura 6.6: Spettro energia-similarita-popolazione per i minimi config-urazionali a T = 0.05. Il raggio dei cerchi e proporzionale alla popo-lazione, il simbolo + al centro del cerchio individua il punto (χ, V )dove χ e la similarita con lo stato nativo, V e l’energia del minimo.Delimitata da linee tratteggiate la regione entro la quale i minimi localisono considerati metastabili. E visibile un minimo metastabile.

0 0.1 0,2 0,4 0,6 0,8

χ-25

-20

-15

-10

-5

0

VStato Nativo

T=0.15

zona metastabile

Figura 6.7: Spettro energia-similarita-popolazione per i minimiconfigurazionali a T = 0.15. Sono visibili 3 minimi metstabili

73

-30 -25 -20 -15 -10 -5 0

0

2

4

6

8

0 0,2 0,4 0,6 0,8-35

-30

-25

-20

-15

-10

-5

0

5

V

%

V

χ

zona metastabile

Figura 6.8: Spettro delle popolazioni e spettro energia-similarita peruna serie di simulazioni a T = 0.15 per la sequenza completamenteidrofoba (B)21. Il minimo energetico assoluto e stato assunto comestato nativo di riferimento anche se non risulta essere coincidente conil massimo della popolazione.

Le tre configurazioni metastabili possiedono locali alterazioni sugli an-

goli diedrali (spesso nel tornante della forcina) che le rendono distin-

guibili nella classificazione, anche se nelle immagini visibili in Figura

6.9 appaiono sostanzialmente uguali.

Come si osserva in Figura 6.3 le fluttuazioni termiche a T = 0.15 fanno

si che χ oscilli prevalentemente all’interno della regione metastabile.

L’istantanea a t = 200τ fotografa il filamento in una delle sue possibili

permutazioni; la frequenza con la quale ricorre ognuna di esse e correla-

bile alla permanenza temporale in quella configurazione in accordo con

la teoria di Kramers [38] sulla probabilita di transizione e in qualche

modo ricalca l’ipotesi di metastabilita avanzata da Thirumalai.

6.3.4 Sequenza idrofoba

Per avere un raffronto qualitativo sulla bonta del ripiegatore Beta21 e

stata eseguita una ulteriore serie di simulazioni a T = 0.15 per la se-

quenza completamente idrofoba (B)21 rappresentata in Figura 6.2(c).

74

Figura 6.9: Configurazioni metastabili a T = 0.15. La differen-za energetica e dell’ordine 20 kb T . La mutua similarita e compresanell’intervallo [0.05, 0.1].

Gli spettri di popolazione ottenuti riportati in Figura 6.8 delineano la

radicale differenza del panorama energetico fra buono e cattivo ripie-

gatore e l’assenza per quest’ultimo di una configurazione probabilisti-

camente dominante.

6.4 Risultati

Una spiegazione teorica riguardo alla dinamica evolutiva del sistema

all’interno del panorama energetico esula dagli obiettivi preliminari di

questo lavoro.

Lo studio presentato in questo ultimo Capitolo inoltre e un’ analisi

meccanica e non termodinamica del processo.

Una distinzione non trascurabile sta nel fatto che la classificazione dei

minimi qui compiuta si basa sull’energia interna e non sull’energia lib-

era di Gibbs a cui fa riferimento l’ipotesi termodinamica discussa nel

Paragrafo 6.1.

Si puo comunque pensare che le visioni dello scenario da queste due

ottiche non siano troppo dissimili o comunque complementari. Il fatto

ad esempio che il minimo meccanico assoluto sia anche il piu frequente

lo rende identificabile anche come il minimo assoluto dell’energia lib-

era.

Le simulazioni condotte danno risulati in accordo con le aspettative

del modello per quanto riguarda la presenza di una configurazione

preferenziale, ossia l’esistenza di uno stato nativo.

75

T χ<0.05 χ<0.1

0.05

0.1

0.15

18 % 25 %

25 % 29 %

30 % 63 %

χ>0.7

7 %

10 %

11 %

0.15 3 % 3 % 81 %

beta21

sequenzaidrofoba

Figura 6.10: Percentuale delle simulazioni con Beta21 che al ter-mine del processo (t = 200τ) sono in una configurazione identica(χ 6 0.05), molto simile (χ 6 0.1, zona metastabile) o completamentedissimile (χ > 0.7) allo stato nativo. Nell’ultimo riquadro in grigioi valori per le simulazioni compiute con la sequenza completamenteidrofoba (B)21.

Dalla tabella riportata in Figura 6.10 emerge una valutazione quan-

titativa della bonta del ripiegatore Beta21 tramite le percentuali di

successo nel processo di ripiegamento divise in fasce di similarita in

contrapposizione con quelle per la sequenza idrofoba (B)21.

Ricordiamo che per ∆χ 6 0.1 due strutture appaiono indistinguibili.

Per quanto riguarda l’aspetto termodinamico della simulazione in questo

lavoro non si e compiuta un’analisi approfondita sulla cruciale influen-

za della temperatura dell’ ambiente. Osservando i dati nella tabelle

6.10 ci limitiamo a concludere che il processo di ripiegamento appare

molto sensibile al cambio di T indice dell’incidenza di un fattore en-

tropico sulla scelta della configurazione di minimo.

Il valore della temperatura per la quale nelle nostre simulazioni si ha il

piu alto successo (T = 0.15 corrispondente circa alla temperatura am-

biente) risulta coerente con gli esperimenti reali compiuti sui filamenti.

76

Conclusioni

La modellizzazione della proteina necessita la costruzione di un algo-

ritmo per la dinamica di una catena di punti connessi che identificano

gli amminoacidi.

La struttura semi-rigida di questo oggetto si presta ad essere parametriz-

zata mediante angoli e la sua dinamica fatta evolvere tramite rotazioni

nello spazio.

L’utilizzo dei quaternioni che si fa all’interno dell’algoritmo proposto

rappresenta il contributo originale dato in questo lavoro di tesi.

La topologia della catena composta di elementi unidimensionali (aste)

fa sı che una completa parametrizzazione in senso quaternionale sia

ridondante; essa potrebbe rivelarsi particolarmente utile per una trat-

tazione degli amminoacidi come corpi rigidi non puntiformi.

Le variabili dinamiche, nel nostro caso, restano i vettori che individ-

uano i segmenti della catena e le relative velocita angolari; i quater-

nioni intervengono all’interno dello schema producendo l’avanzamento

rotazionale delle posizioni.

Il secondo contributo di questo lavoro nella creazione dell’integratore

sta nel sistema di calcolo delle forze vincolari che si propagano lungo

la catena e la mantengono connessa.

L’integratore e strutturato in modo da raggiungere localmente (sulle

posizioni relative) un errore al terzo ordine preservando implicitamente

le distanze tra elementi contigui.

I test di prova compiuti per il caso di una catena di N pendoli connes-

si mostrano una deriva energetica lineare dell’ordine di una parte su

centomila dopo un periodo di 100τ con h = 0.005τ (N 6 10). L’errore

cresce linearmente con l’aumento del numero di elementi.

77

I tempi di calcolo in presenza di campo costante crescono linearmente

con N , quadraticamente nel caso di interazioni a coppie per le quali

e necessaria la misura delle distanze reciproche, raddoppiano nel pas-

saggio dal secondo al terzo ordine nello schema di integrazione.

La non triviale introduzione di una componente stocastica nella di-

namica (rumore o bagno termico) e gestita sull’impronta della dinam-

ica di Langevin nel limite di bassa frizione.

Opportune modifiche sono apportate nel trasporto della trattazione

dalle variabili cartesiane (x , v) a quelle della catena (e, ~ω).

Come verifica del procedimento i test sul moto browniano di un fila-

mento immerso in un solvente a temperatura costante riproducono con

buona fedelta le caratteristiche statistiche attese teoricamente (dis-

tribuzione maxwelliana delle velocita, equipartizione dell’energia).

Utilizzando tale integratore viene costruito il modello per la proteina.

Gli esperimenti condotti sulla sequenza campione Beta21 rivelano la

presenza di un’alta probabilita di comparsa della caratteristica β-sheet

corrispondente allo stato nativo della sequenza.

Tale probabilita e modulata dalla temperatura del bagno termico.

Ad esempio per una simulazione isoterma a temperatura T = 0.15

(circa la temperatura ambiente) il 30% dei processi di ripiegamento

conducono a configurazioni che differiscono per meno del 5% 2 dallo

stato nativo, il 63% per meno del 10%, solo il 10% per piu del 70%.

Per la sequenza (B)21 tipicamente associata ad un cattivo ripiegatore,

identificando lo stato nativo con il minimo energetico assoluto si ot-

tiene che circa il 3% dei processi conducono a strutture che differiscono

per meno del 10%, almeno l’81% per piu del 70%.

Il modello che si e costruito da risultati che confermano l’affidabilita

dell’integratore accoppiato ad una componente stocastica (dinamica

di Langevin) e che si dimostrano coerenti e in linea con le simulazioni

effettuate in altri lavori di ricerca e con gli esperimenti reali.

2La percentuale e ricavata dalla funzione di similarita χ.

78

Appendice A

Metodo Verlet L’algoritmo rappresenta la soluzione diretta dell’e-

quazione mr = f per la trattazione di moti di oggetti puntiformi.

L’evoluzione della posizione r(t+ h) si basa sulla posizione r(t), l’ac-

celerazione a(t) e la posizione r(t− h) del passo precedente

r(t+ h) = 2r(t) − r(t− h) + h2a(t) (6.2)

Come si puo notare la velocita non compaiono poiche sono state elim-

inate utilizzando lo sviluppo di Taylor

r(t+ h) = r(t) + hv(t) +h2

2a(t) + ...

r(t− h) = r(t) − hv(t) +h2

2a(t) − ...

(6.3)

La velocita al tempo corrente v(t) potrebbe essere necessaria per cal-

colare la traiettoria nel caso di forze dipendenti da essa (magnetiche,

per esempio) o per stimare l’energia cinetica.

Essa puo essere ottenuta

v(t) =r(t+ h) − r(t− h)

2h(6.4)

L’equazione 6.2 e corretta con una stima locale dell’ordine O(h4) men-

tre le velocita lo sono con O(h2).

La caratteristica del metodo Verlet e quella di essere centrato; r(t+h)

e r(t− h) giocano un ruolo simmetirco e rendono l’algoritmo tempo-

ralmente reversibile.

Lo schema logico della procedura e illustrato in Figura 6.11.

Leap-frog I limiti del metodo Verlet risiedono nell’imprecisione con

la quale vengono maneggiate le velocita e dalla deriva numerica che

79

t−h t+h t t−h t+h t t−h t+h t

rv

a

Figura 6.11: Schema rappresentativo del metodo Verlet. Nei riquadriposizione r, velocita v, accelerazione a relative ai tempi t−h, t, t+h.Le frecce indicano il calcolo a partire da grandezze note (in grigio).

emerge dall’equazione 6.2 dove un piccolo termine dell’ordine h2 e

sommato ad una differenza di grandi termini dell’ordine h0.

Nel 1970 Hockney propose la variante denominata leap-frog, salto della

rana.

Lo schema si struttura

r(t+ h) = r(t) + hv(t+h

2) (6.5a)

v(t+h

2) = v(t− h

2) + ha(t) (6.5b)

Le grandezze archiviate sono la posizione e l’accelerazione correnti

r(t), a(t) e la velocita v(t− h2).

L’equazione 6.5b aggiorna la velocita al semi-passo successivo che

viene utilizzata in 6.5a per il calcolo della nuova posizione.

La velocita corrente puo essere calcolata

v(t) =1

2[v(t+

h

2) + v(t− h

2)] (6.6)

L’eliminazione di v dal sistema 6.5a,6.5b mostra come l’algoritmo sia

algebricamente equivalente al metodo Verlet.

Velocity-Verlet L’equazione 6.6 mostra ancora come la velocita

siano maneggiate in modo non completamente soddisfacente.

Nel 1982 (Swope-Andersen-Berens-Wilson) viene proposta la variante

detta velocity-Verlet

r(t+ h) = r(t) + hv(t) +h2

2a(t) (6.7a)

80

r

v

a

t−h t t+h t−h t t+h t−h t t+h

Figura 6.12: Schema rappresentativo della variante Leap-frog anal-ogo a quello in Figura 6.11. I riquadri interposti rappresentano legrandezze ai midstep.

t t+h

rv

a

t+h/2t t+h t+h/2 t t+h t+h/2

Figura 6.13: Schema rappresentativo della variante Velocity-Verlet.

v(t+ h) = v(t) +h

2[a(t) + a(t+ h)] (6.7b)

Ancora una volta eliminando v ci si riconduce al comune Verlet.

L’algoritmo richiede l’archiviazione di r(t), v(t), a(t) e si svolge in due

tappe 3.

Da 6.7a la posizione e fatta evolvere a r(t+h) e la velocita al semipasso

e calcolata

v(t+h

2) = v(t) +

h

2a(t) (6.8)

Si ricavano cosı le velocita aggiornate

v(t+ h) = v(t+h

2) +

h

2a(t+ h) (6.9)

La stabilita numerica di questa variante fa si che essa sia la piu utiliz-

zata per opportune simulazioni computazionali.

3in modo simile all’algoritmo Gear Predictor-Corrector

81

82

Appendice B

Metodo di integrazione quaternionale Considerando il fatto che

non si possono ottenere equazioni prive di singolarita numeriche uti-

lizzando solo 3 variabili (angoli di Eulero, per esempio) Evans [12]

nel 1977 propose l’uso delle 4 variabili quaternionali come coordinate

generalizzate.

Come e stato discusso nel Capitolo 1 il quaternione

q = (q0, q1, q2, q3)

e soggetto al vincolo

q20 + q2

1 + q22 + q2

3 = 1 (6.10)

Riconducendoci ai convenzionali angoli di Eulero, le coordinate qi sono

definite

q0 = cosθ

2cos

φ+ ψ

2

q0 = sinθ

2cos

φ− ψ

2

q0 = sinθ

2sin

φ− ψ

2

q0 = cosθ

2sin

φ+ ψ

2

(6.11)

e soddisfano l’equazione differenziale

q0q1q2q3

=1

2

q0 −q1 −q2 −q3q1 q0 −q3 q2q2 q3 q0 −q1q3 −q2 q1 q0

0ωx

ωy

ωz

(6.12)

83

scritta in forma quaternionale

q =1

2q · ω (6.13)

dove ω rappresenta il vettore velocita angolare nel riferimento sol-

idale con il corpo.

Le equazioni 6.12 vanno accoppiate con le equazioni di Eulero per la

dinamica del corpo rigido

ωx =τ bx

Ixx

+

(

Iyy − Izz

Ixx

)

ωx ωz x→ y → z (6.14)

dove τ b e il momento delle forze nel sistema solidale.

Insieme rappresentano un sistema differenziale al primo ordine risolu-

bile numericamente.

Come abbiamo gia detto la dipendenza posizionale della derivata po-

sizionale in 6.12 rende inapplicabile il comune metodo Verlet (vedi

Appendice A); in secondo luogo vi e la questione del mantenimento

del vincolo sulla norma di q.

Una formulazione stile leap-frog per l’integrazione delle equazioni 6.12

e 6.14 e stata avanzata nel 1981 da Fincham [13].

~Ls(t) = ~Ls(t− h2) + h

2~τ s(t)

q(t+ h) = q(t) + hq(t+ h2)

(6.15)

dove la prima equazione coinvolge il momento angolare ~Ls riferito

al sistema esterno per il quale vale la legge dinamica Ls = τ s.

Le difficolta annidate nel sistema 6.15 sono principalmente la necessita

di operare le trasformazioni tra sistema solidale ed esterno

Ls = q Lb q

e tra momento e velocita angolari

ωb = I−1Lb

dove I e il tensore costante di inerzia nel sistema solidale.

Cio da luogo ad una equazione implicita da risolversi iterativamente e

84

con l’ausilio di una estrapolazione4.

Il vincolo inoltre e preservato forzatamente con una riscalatura delle

coordinate qi al termine di ogni passo di integrazione.

Recentemente [7] Omelyan presenta un metodo di integrazione

dove tale vincolo e preservato implicitamente e che fa uso solo di in-

terpolazioni e non di estrapolazioni.

In modo simile a 6.15 il momento angolare viene fatto avanzare al

semipasso Ls(t + h2) quindi si ricava la velocita angolare nel sistema

solidale

ωb(t+h

2) = I−1q(t+

h

2) Ls(t+

h

2) (6.16)

L’orientazione aggiornata del corpo si ottiene con una trasfor-

mazione ortonormale nello spazio dei quaternioni Q.

q(t+ h) =

[

1 − h

2ωb(t+

h

2)

]−1 [

1 +h

2ωb(t+

h

2)

]

q(t) (6.17)

dove 1 = (1, 0, 0, 0) e ωb = (0, ωx, ωy, ωz) nello spazio Q.

Le equazioni 6.16 e 6.17 si risolvono iterativamente considerando che

in 6.16 ci si avvale dell’interpolazione

q(t+h

2) =

1

2[q(t) + q(t+ h)] (6.18)

L’algoritmo di Omelyan raggiunge una precisione e una stabilita

confrontabili a quelle ottenute tramite un approccio atomico.

Matrici di rotazione (metodo R-SHAKE) In contrapposizione

con l’utilizzo dei quaternioni o degli angoli di Eulero per i quali la

4predizione di una variabile ad un tempo avanzato ottenuta con errori di ordinepiu basso

85

matrice di rotazione A viene parametrizzata 5

A =

q20 + q2

1 + q22 + q2

3 2(q1q2 − q0q3) 2(q1q3 + q0q2)2(q1q2 + q0q3) q2

0 − q21 + q2

2 − q23 2(q2q3 − q0q1)

2(q1q3 − q0q2) 2(q2q3 + q0q1) q20 − q2

1 − q22 + q2

3

=

cosψ cosφ− cos θ sinφ sinψ − sinψ cosφ− cos θ sinφ cosψ sin θ sin φcosψ sinφ+ cos θ cos φ sinψ − sinψ sinφ+ cos θ cosφ cosψ − sin θ cosφ

sin θ sinψ sin θ cosψ cos θ

(6.19)

in questo metodo le 6 varaibili indipendenti di A vengono fatte

evolvere direttamente sotto il vincolo

AAT = 1 (6.20)

Scriviamo la Lagrangiana per un corpo rigido rotante

L =M

2rcmrcm +

1

2Tr[AIA

T]−V (A, rcm) +Tr[Λ(AAT − 1)] (6.21)

dove i primi due termini rappresentano le energie cinetiche traslazionale

e rotazionale, rcm e la posizione del centro di massa, M e la massa

totale, I e il tensore di inerzia, il terzo termine e il potenziale totale e

l’ultimo termine e stato aggiunto per assicurare il mantenimento del

vincolo 6.20 tramite la matrice simmetrica Λ di moltiplicatori di La-

grange (6 condizioni di vincolo).

Per costruire le equazioni del moto definiamo il momento coniugato

(matrice 3 × 3)

Π ≡ ∂L∂A

= AI (6.22)

Le equazioni di Hamilton risultano

A = Π · I−1

Π = −∂AV + 2A · Λg(A) = AAT − 1 = 0

(6.23)

5la matrice A trasporta le coordinate dal sistema solidale al sistema esternoxs = Axb.

86

A questo punto si applica l’algoritmo R-SHAKE (derivato da SHAKE)

Π(t+h

2) = Π(t− h

2) + h[2A(t)Λ(t) − ∂A(t)V ]

A(t+ h) = A(t) + h Π(t+h

2) · I−1

g[A(t+ h)] = 0

(6.24)

L’errore locale ad ogni passo di integrazione e dell’ordine O(h3).

Con questo metodo il vincolo e imposto forzatamente mentre e sarebbe

possibile [10] costruire una evoluzione esponenziata di A che lo man-

tenga implicitamente con lo svantaggio del calcolo dell’esponenziale

matriciale ad ogni passo per ogni molecola.

L’algoritmo R-SHAKE, pur presentando un errore locale maggiore per

brevi simulazioni per un fissato passo di integrazione rispetto ad un

algoritmo quaternionale, mostra maggiore stabilita e assenza di derive

energetiche per simulazioni di lunga durata a causa del suo carattere

simplettico.

87

88

Bibliografia

[1] Anfinsen C. Science 181, 223 (1973)

[2] Anfinsen C. Scheraga H.A. Adv.Prot.Chem. 29, 205 (1975)

[3] Branden C. Tooze J. Introduzione alla struttura delle proteine

Zanichelli 1993

[4] Clementi C. Effective potentials for protein folding models SISSA

1998

[5] Hamilton W.R. On Quaternions Proceedings of the Royal Irish

Academy 3, 1-16 (1847)

[6] Kuipers J.B. Quaternions and rotation sequences, Princeton

University Press, 1999

[7] Omelyan I.P. Algorithm for numerical integration of the rigid-body

equations of motion Physical Review E 58, 1169 (1998)

[8] Ryckaert J.P. Ciccotti G. Berendsen H.J.C. J. Comput. Phys. 23,

327 (1977)

[9] Andersen H.C. J. Comput. Phys. 52, 24 (1983)

[10] Ahlrichs R. Brode S. Comput. Phys. Commun. 42, 59 (1986)

[11] Dullweber A. Leimkuhler B. McLachlan R. J. Chem. Phys. 107,

5840 (1997)

[12] Evans D.J. Mol.Phys. 34 , 317 (1977)

[13] Fincham D. Mol.Simul. 8 , 165 (1992)

89

[14] Omelyan I.P. A new leapfrog integrator of rotational motion, Jan

1999

[15] Kol A. Laird B. Leimkuhler B. A symplectic Method for Rigid-Body

Molecular Simulation, March 1997

[16] Kneller G.R. Hinsen K. Phys.Rev. E 50, 1559 (1994)

[17] Kalra D. A formulation of Rigid-Body Assemblies for Computer

Graphics Modeling

[18] Honeycutt J.D. Thirumalai D. Metastability of the folded states of

globular proteins Proc.Natl.Acad.Sci.USA 87 3526 (1990)

[19] Baumketner A. Hiwatari Y. Phys.Rev. E 66, 11905 (2002)

[20] Sorenson J.M. Head-Gordon T. Proteins 37, 582 (1999)

[21] Shea J.E. et al. J. Chem. Phys. 109, 2895 (1998)

[22] Klimov D.K. Thirumalai D. Phys.Rev.Lett. 79, 317 (1997)

[23] Stillinger F.H. Head-Gordon T. Hirshfeld C.L. Phys.Rev. E 48,

1469 (1993)

[24] Kalra D. A general formulation of rigid-body assemblies for computer

graphics modeling

[25] Risken H. The Fokker-Planck Equation Springer-Verlag (1989)

[26] Turchetti G. Dinamica Classica dei sistemi fisici Zanichelli 1998

[27] Landau L.D. Lifsits E.M. Fisica statistica Editori Riuniti 1978

[28] Guo Z. Thirumalai D. Kinetics of protein folding: nucleation

mechanism, time scales and pathways Biopolymers, 36,83 (1994)

[29] Levinthal C. J. Chem. Phys. 65, 44 (1968)

[30] Nymeyer H. Garcıa A. Onuchic J.N. Proc. Natl. Acad. Sci. USA

95, 5921 (1998)

90

[31] Rapaport D.C. Phys. Rev. E 66, 011906 (2002)

[32] Baumketner A. Hiwatari Y. Phys.Rev. E 66, 011905 (2002)

[33] Micheletti C. Seno F. Maritan A. Banavar J.R. Computational

Materials Science 20, 305 (2001)

[34] Micheletti C. Seno F. Maritan A. Banavar J.R. Learning effective

amino acid interactions through iterative stochastic techniques

[35] Klimov D.K. Thirumalai D. Phys.Rev. Lett. 79, 317 (1997)

[36] Veitshans T. Klimov D.K. Thirumalai D. Folding and design 2,

1 (1997)

[37] Zhou Y. Karplus M. J. Mol. Biol. 293, 917 (1999)

[38] Sebastian K.L. Alok K.R. Paul Phys. Rev. E 62, 927 (2000)

91