Dispense del corso Modellistica applicata a molecole

di interesse biologico

Corso di Laurea Magistrale in Biotecnologie

Molecolari

Piero Procacci

14 novembre 2017

Indice

1 I Force Fields atomistici per Sistemi di interesse biologico 2

1.1 Potenziale bonded (legato o di valenza) per i sistemi biologici . . . . . . . 6

1.2 Potenziale non bonded . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2 Elementi di Meccanica Statistica 30

2.1 Equilibrio termodinamico, funzione di distribuzione e spazio delle fasi . . 35

2.2 Formulazione microscopica di primo e secondo principio . . . . . . . . . . 38

2.3 Insieme microcanonico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

2.4 Insieme Canonico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3 Funzione Partizione ed Equilibrio Chimico 48

3.1 Gas perfetti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.2 Funzione di partizione molecolare . . . . . . . . . . . . . . . . . . . . . . 51

3.3 Equilibrio chimico per una miscela ideale di gas . . . . . . . . . . . . . . 59

4 Costante di Affinita Farmaco-Proteina 69

4.1 Teoria RRHA per il calcolo dell’energia libera di dissociazione . . . . . . 74

4.2 Cenni sul calcolo dell’energia di interazione farmaco-proteina in solvente

implicito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

1

Capitolo 1

I Force Fields atomistici per Sistemi

di interesse biologico

Nel capitolo che segue descriveremo le metodologie per il calcolo delle energie configu-

razionale dei sistemi di interesse biologico a livello atomistico. La conoscenza dell’energia

configurazionale del sistema costituisce, come vedremo, un elemento fondamentale per

individuare le strutture molecolari assunte dal sistema in condizioni fisiologiche.

Per “force field ”si intende l’insieme dei parametri o parametrizzazione che con-

corrono alla determinazione dell’energia di interazione tra le particelle (atomi) che

costituiscono il sistema. Il generale possiamo dunque scrivere:

V (r1, r2, ..., rN) = E (1.1)

dove r1, r2, ..., rN sono i vettori1 posizione della particella (atomo) 1, della particella

2, ..., della particella ennesima. La funzione energia definita dall’equazione 1.1 e dunque

1N.B. I vettori sono definiti nello spazio cartesiano da tre componenti, lungo x, y e z. Quando ci si

riferisce ad una quantita vettoriale (per esempio velocita, forza, posizione etc) la si scrive in grassetto

per distinguerla da una quantita scalare. Dunque con r indico il vettore posizione (tre coordinate x,y,z)

rispetto all’origine, mentre con r indico la distanza dall’origine (un numero)

2

una funzione a 3N variabili perche il vettore posizione della particella i-esima si definisce

specificando le sue componenti x, y, z.

Il sistema biologico sul quale ci concentreremo in questo corso e costituito da un

recettore (tipicamente una proteina), da un farmaco (una molecola molto piu piccola)

e dal solvente circostante (normalmente acqua). Se Np e il numero di atomi che costi-

tuiscono la proteina, Nf il numero di atomi che costituiscono il farmaco e Ns il numero

di atomi che costituiscono il solvente, tali che N = Np + Nf + Ns, l’equazione 1.1 puo

dunque essere scritta come:

V (r1, r2, ..., rN) = Vp(r1, r2, ..., rNp) + Vf (r1, r2, ..., rNf ) +

+ Vpf(r1, r2, ..., rNp , r1, r2, ..., rNf) +

+ Vps(r1, r2, ..., rNp , r1, r2, ..., rNs) +

+ Vfs(r1, r2, ..., rNf , r1, r2, ..., rNs) +

+ Vs(r1s, r2s, ..., rNs) (1.2)

Sulla base della Eq. 1.2, la forma potenziale 1.1 puo essere scritta come una somma

di termini corrispondenti a:

1) Vp(r1, r2, ..., rNp) Potenziale proteina-proteina che dipende dalle coordinate dei soli

Np atomi delle proteina.

2) Vf (r1, r2, ..., rNf ) Potenziale farmaco-farmaco che dipende dalle coordinate dei soli Nf

atomi del farmaco.

3) Vpf(r1, r2, ..., rNp , r1, r2, ..., rNf ) Potenziale farmaco-proteina che dipende dalle Npf =

Np +Nf coordinate della proteina e del farmaco

4) Vps(r1, r2, ..., rNp , r1, r2, ..., rNs) PotenzialeNps = Np+Ns proteina-solvente che dipende

dalle coordinate della proteina e del solvente.

5) Vfs(r1, r2, ..., rNf , r1, r2, ..., rNs) Potenziale farmaco-solvente che dipende dalle Nfs =

Nf +Ns coordinate del farmaco e del solvente.

3

6) Vs(r1s, r2s, ..., rNs) Potenziale solvente-solvente che dipende dalle Ns coordinate del

solo solvente.

Nella reazione farmaco-recettore, ovvero ai fini della determinazione dell’affinita’ di

un farmaco verso un recettore, ciascuno di questi termini del potenziale di interazione e

importante.

P + F PF −→ [P][F]

[PF]= ki (1.3)

La funzione 1.2 rende il valore dell’energia per una configurazione del sistema costi-

tuito dal farmaco, dalla proteina e dal solvente. Per configurazione del sistema si intende

l’insieme dei 3N valori delle coordinate che costituiscono il sistema stesso composto da N

particelle (atomi). Come vedremo le configurazioni ”scelte” dal sistema biologico (cioe

quelle che e’ piu probabile osservare in condizioni fisiologiche) sono quelle a energia

inferiore. Vedremo meglio in seguito che nei sistemi a temperatura finita, viene ”popo-

lata” maggiormente la configurazione con l’energia libera piu bassa. L’energia libera si

definisce come:

G = E − TS (1.4)

Essa e una funzione delle coordinate di ordine sistema. Nel caso del sistema farmaco-

recettore descritto dalla Eq. 1.2, essa e’ una funzione delle N = Np+Nf+Ns coordinate.

G(r1, r2, ..., rN) = E(r1, r2, ..., rN)− TS(r1, r2, ..., rN) (1.5)

Il senso della notazione G(r1, r2, ..., rN) e’ quello di indicare G come una funzione

non esplicita di tutte coordinate del sistema. Come vedremo meglio in seguito, G e’ una

funzione integrale che viene cioe’ calcolata integrando su tutte le coordinate. Quanto al

significato della notazione E(r1, r2, ..., rN), essa e direttamente collegabile al potenziale

4

definito nelle equazioni 1.1 e 1.2:

E(r1, r2, ..., rN) =< V (r1, r2, ..., rN) > (1.6)

Nella Ee. 1.6, la notazione con le parentesi angolari, < V (r1, r2, ..., rN) >, sta ad

indicare un un integrale su tutte le coordinate del sistema, corrispondente al valor medio

del potenziale (in senso termodinamico statistico) nelle condizioni del sistema, come

vedremo in seguito. Per quanto riguarda il solvente, esistono due modi per tenerlo in

conto. 1) Il modo esplicito, nel quale il solvente e semplicemente costituito da un insieme

di molecole fatte a loro volta di atomi. 2) Il metodo implicito, nel quale il solvente e

considerato un mezzo dielettrico continuo in cui sono immersi il farmaco e la proteina.

I potenziali Vp e Vf nella Eq 1.2 vengono, in questo ultimo approccio, “rinormalizzati”.

L’ equazione 1.2 si riferisce al modello a solvente esplicito costituito da Ns particelle

(atomi). L’energia complessiva nel modello implicito si scrive come

V ′(r1, r2, ..., rN) = V ′p(r1, r2, ..., rNp) + V ′f (r1, r2, ..., rNf ) +

+ V ′pf(r1, r2, ..., rNp , r1, r2, ..., rNf) (1.7)

Il simbolo ’ nell’equazione 1.7 come soprascritto di Vp, Vf e Vpf , sta ad indicare che

queste funzioni sono opportunamente modificate rispetto al caso esplicito per la presen-

za del solvente continuo. Prima di affrontare il problema farmaco-recettore da un punto

vista meccanico-statistico (cioe dal punto di vista delle equazione 1.5 e 1.6), bisogna

vedere come fare a calcolare il valore dell’energia data una configurazione del sistema,

ovvero noto l’insieme delle coordinate x, y, z di tutti i suoi atomi. Per fare questo oc-

corre definire esplicitamente una forma analitica della funzione empirica V (r1, r2, ..., rN).

Affronteremo il problema a ”pezzi”. Vedremo prima come sono definiti i potenziali

Vp(r1, r2, ..., rNp) e Vf (r1, r2, ..., rNf ) cosiddetti ”bonded”, cioe che definiscono la strut-

tura topologica delle molecole proteina e farmaco. Vediamo prima il caso della proteina:

5

Vp(r1, r2, ..., rNp) = V bondedp (r1p , r2p , ..., rNp) + V non−bonded

p (r1p , r2p , ..., rNp) (1.8)

Il potenziale V bondedp (r1, r2, ..., rNp) e detto termine “bonded ”“di valenza ”e coinvolge

tutti gli Np atomi delle proteine. Il termine V nonbondedp (r1, r2, ..., rNp) e detto potenziale

“non bonded ”e di nuovo coinvolge le coordinate di tutti gli atomi della proteina. Esso

e responsabile delle interazioni tra atomi non legati e lo vedremo meglio in seguito.

Possiamo naturalmente scrivere una equazione analoga alla equazione 1.8 per il farmaco

o per il solvente :

Vf (r1, r2, ..., rNf ) = V bondedf (r1, r2, ..., rNf ) + V non−bonded

f (r1, r2, ..., rNf )

Vs(r1, r2, ..., rNs) = V bondeds (r1, r2, ..., rNs) + V non−bonded

s (r1, r2, ..., rNs) (1.9)

Sulla base dell’equazione 1.8 (o della equazione 1.9) abbiamo assunto che il potenziale

intra-proteina (o intra-farmaco/solvente) possa essere scritto come una somma di due

termini, il potenziale bonded e quello non-bonded (I termini misti nella Equazione 1.2,

cioe i termini Vpf , Vps, Vfs in solvente esplicito o il termine V ′pf in solvente implicito,

come vedremo nelle sezioni seguenti, sono descritti dal solo potenziale “non-bonded”).

Concentriamoci ora sulla forma del potenziale bonded.

1.1 Potenziale bonded (legato o di valenza) per i

sistemi biologici

Per capire questo potenziale la cosa migliore e fare un esempio pratico. Supponiamo di

aver a che fare con una “mini proteina ”costruita da una sola alanina (vedi Figura 1.1)

Il fatto che questo esempio sia sufficiente a spiegare la forma funzionale del potenziale

“bonded ”, implica il concetto di trasferibilita dei parametri associati a tale potenziale.

Cosı in generale possiamo scrivere che :

6

V bonded(r1, r2, ..., rNp) = V bonded(r1, r2, ..., rNp , αb1, α

b2, ..., α

bp) (1.10)

Dove gli αb1, αb2, ..., α

bp sono i “parametri ”del potenziale bonded. Questi parametri

costituiscono una parte del cosiddetto “force field”. L’insieme di tutti i parametri riferiti

ai diversi termini, bonded, non-bonded, Vf , Vp, Vps, etc, etc. costituiscono il “force field”.

Figura 1.1: Mini-proteina costituita dal solo residuo Alanina con due unita terminali metiliche.

La mini-proteina ha tutti gli ingredienti di una proteina, il legame peptidico (evi-

denziato in rosso), gli angoli diedri psi e phi, etc. etc.

7

Vedremo ora, nel caso dell’interazione bonded per la proteina, le struttura della

funzione V bonded(r1, r2, ..., rNp , αb1, α

b2, ..., α

bp). Abbiamo che:

V bonded(r1, r2, ..., rNp , αb1, α

b2, ..., α

bp) = V stretch(r1, r2, ..., rNp , α

st1 , α

st2 , ..., α

stp ) +

+V bend(r1, r2, ..., rNp , αben1 , αben2 , ..., αbenp ) +

+V tors(r1, r2, ..., rNp , αtors1 , αtors2 , ..., αtorsp ) (1.11)

Una formula del tutto simile si applica al potenziale di tipo bonded del farmaco che

in questo caso dipendera dalle r1, r2, ..., rNf coordinate del farmaco e dagli αb1, αb2, ..., α

bf

parametri di tipo bonded per il farmaco.

Vediamo il primo termine nella equazione 1.11:

V stretch(r1, r2, ..., rN , αst1 , α

st2 , ..., α

stp ) =

1

2

Nstrech∑i=1

Ki(ri − r0i)2 (1.12)

La sommatoria nella equazione 1.12 va su tutti gli stretching presenti nel sistema,

ovvero su tutti i legami. Nelle mini proteine di Figure 1.1, r0i contiene 21 “stretching’,

di cui 10 di tipo C-H, 2 di tipo C=O, 2 di tipo N-H, 4 di tipo C-N, 1 di tipo C-C,

2 di tipo C-C(=O). Io ho quindi raggruppato gli stretching possibili nelle molecole di

Figura. 1. in “classi” Nella 1.12, teoricamente ci sono tante coppie di parametri roi e Ki

quanti sono gli stretching del sistema. In realta il numero di r0i e Ki e pari al numero

di classi. Nel caso della nostra proteina di Figura 1.1 le classi sono solo 6 (contro 21

stretching). Cio significa che i parametri relativi ad uno stretching C-H sono trasferibili

a tutti gli stretching appartenenti alla stessa classe. Cio che caratterizza una “classe” di

stretching e il tipo (atomo) e ibridazione dei due atomi coinvolte nello stretching: tutti

e 9 gli stretching C-H coinvolgono un carbonio eh un idrogeno, nei quali il carbonio

e ibridato sp3. Veniamo ora al significato dei parametri r0i e Ki nelle equazione 1.12

Come detto vi sono 6 classi di stretching nella proteina di Figura. 2. e dunque vi

sono 12 parametri, α1, ..α12, nella nostra molecola di Figura 1.1. Per capire a cosa

sono legati questi parametri di stretching, facciamo il grafico della funzione potenziale

8

Figura 1.2: Destra: Grafico della funzione dell’energia di uno stretching generico Vstr =

12K(r−r0)2 in funzione della distanza r tra i due atomi legati.Sinistra: Definizioni di coordinate

vettoriali di due atomi

associate ad uno solo di questi parametri. ro e la distanza di equilibrio, per la quale

l’energia di stretching vale 0. Per tutti gli altri valori di r, V e maggiore di 0. Dunque la

configurazione piu favorita energeticamente e proprio quella per cui r = r0, cioe quella

per la quale i due atomi coinvolti nello stretching si trovano proprio alla distanza r0. r

e la distanza tra i due atomi coinvolti nello stretching. Essa dipende dalle coordinate di

questi atomi (supponiamo 1 e 2) in modo tale che:

r = |r1 − r2| =√

(x1 − x2)2 + (y1 − y2)2 + (z1 − z2)2 (1.13)

La costante K nell’espressione per lo stretching 1.12 e legata alla curvatura della fun-

zione energia, cioe alla forza del legame tra atomi coinvolti nello stretching. Piu alto e K

e maggiore e la forza del legame. In realta l’espressione 1.12 e solo un’approssimazione

della realta fisica di un legame chimico. La funzione energia “vera” in funzione della

distanza di legame assomiglia alla forma funzionale V = 12K(r − r0)2 − Dc soltanto

9

Figura 1.3: Funzione di Morse vera per un legame chimico (tratto solido) e sua approssi-

mazione armonica (curva tratteggiata) V = 12Ki(ri − r0i)

2 − Dc dove Dc e’ un off-set legato

all’energia dissociazione

per valori di distanze vicini al valore di equilibrio r0. Questo accade (cioe la distan-

za di legame oscilla attorno a r0) soltanto quando la temperatura del sistema non e

troppo elevata (tipicamente T=300K). Dunque l’espressione 1.12 costituisce una buona

approssimazione della realta quantomeccanica (rappresentata dalla curva solida di Figu-

ra 1.3 solo a temperature ordinarie, cioe quelle alle quali i legami chimici non si rompono

mai: In una proteina, se non intervengono modifiche post-transazionali indotte da agen-

ti esterni (tipicamente enzimi), l’integrita dei legami chimici e preservata. Il termine

Dc in Figura 1.3 (corrispondente alla ordinata del minimo della funzione) corrisponde

all’energia di dissociazione ed e un termine costante, che viene normalmente non in-

10

cluso nell’espressione energia di stretching: se invece si vuole confrontare l’energia di

stretching di due diverse configurazioni atomiche delle molecole, questo termine (Dc) si

elide e non conta nulla. La costante K e legata alla frequenza di oscillazione del legame.

maggiore e la frequenza di vibrazione del legame e piu elevato e K secondo la relazione:

K = µω2 (1.14)

µ =m1m2

m1 +m2

ω =

√K

µ(1.15)

dove µ e la “massa ridotta”. la frequenza di vibrazione. La K puo essere derivata, per

esempio, dagli spettri infrarossi di molecole organiche. La K relativa allo stretching C-H

puo essere dedotta dall’analisi degli spettri degli idrocarburi aromatici

ω = 2πν −→ 1

λ(cm−1) =

ν

c−→ ν = 3000 cm−1 × 3× 1010cm sec−1

ν ' 9× 1013sec−1 −→ ω ' 3.14× 9× 1013 radianti sec−1 (1.16)

Dunque la frequenza ω dello stretching C-H e pari a ' 3×1014 sec−1. Cio implica che

la K e data da K ' 9×1028 ( u.m.a. sec−2), visto che la massa ridotta del legame C-H e

circa 1, pari cioe a quella dell’atomo piu leggero (in u.m.a.). In realta le unita di misura

della K, “costante di forza”, sono in Kcal mol−1 A−2, dove 1A = 10−10 m. In questa

unita la K del legame C-H desunta dagli spettri infrarossi e data da c.a 200 KcalA−2.

Infatti: 9· 1028 ∼ 1029 g mol−1sec−2 ∼ 1026 Kg mol−1sec−2 ∼ 1026 Joule mol−1metri−2 ∼

106 Joule mol−1A−2 ∼ 103 kJoule mol−1 A−2 ∼ 200 kcal mol−1A−2

Veniamo ora al termine detto di bending espresso dal termine V bend(r1, r2, ..., rNp)

nella Eq. 1.12. Abbiamo che:

V bend(r1, r2, ..., rN , αben1 , αben2 , ..., αbenp ) =

1

2

bendings∑1=1

K(b)i (α− α0i)

2 (1.17)

Al solito K(b)i e la costante di forza (legata alla frequenza) dell’i-esimo bending, α0i e

l’angolo di bending di equilibrio (cioe quello per il quale l’energia di bending e 0) e

11

Figura 1.4: Spettro infrarosso del butano in fase gassosa. L’assorbimento a ca 3000 cm−1

corrisponde ai modi di stretching della molecola.

la sommatoria e’ estesa a tutti i possibili bendings (angoli) della molecola. Nel caso

di Figura 1.1, questa sommatoria si estende a 36 possibili bending (c.a. 1.5 volte gli

stretching) puri. Per ogni carbonio ibridato sp3 noi possiamo definire 6 possibili bending

e per ogni carbonio/azoto ibridato sp2 noi contiamo 3 bending, per un totale appunto

di 36. Le “classi” in questo caso sono solo 10 (ovvero i parametri da definire nella

Equazione 1.17 sono solo 20, i 10 angoli di equilibrio e le 10 costanti di forza legate alla

frequenza di bending per le 10 classi ). Le classi sono definite dal tipo di atomi e dalla

loro ibridazione. Nuovamente la K e legata alla frequenza di bending (che nel butano

[vedi Figura 1.4] e circa 1000 cm−1). La relazione K ↔ ω e un po’ piu complessa che

nel caso dello stretching: ω2 = 4K/Ib dove Ib e il momento di inerzia attorno all’asse

di bending. Il valore dell’angolo α dipende dalle coordinate di tre atomi, ovvero dei tre

12

Figura 1.5: Classi di bending per la mini-proteina di Figura 1.1

atomi coinvolti nel bending. Dunque:

α = (x1, y1, z1, x2, y2, z2, x3, y3, z3) (1.18)

Nella figura 1.1 r12· r13 e il prodotto scalare dei due vettori:

cosα =r12· r13

|r12|· |r13|

r12 = (r2 − r1) = [(x2 − x1), (y2 − y1), (z2 − z1)] (1.19)

r13 = (r3 − r2) = [(x3 − x2), (y3 − y2), (z3 − z2)] (1.20)

Ricordiamo che il prodotto scalare tra due vettori a e b e definito dalle:

a·b = axbx + ayby + azbz

13

Dove ax, ay, az sono le componenti x, y, e z del vettore a e bx, by, bz sono le componenti

x, y e z del vettore b. Nella Figura 1.1, con la notazione|r12| si indica il modulo o norma

del vettore r12 che va dall’atomo 1 all’atomo 2. Esso coincide con la distanza ordinaria

(euclidea) tra l’atomo 1 e l’atomo 2, e cioe r12 = [(x2− x1)2 + (y2− y1)2 + (z1− z2)2]1/2.

Dunque in accordo alla formula riportata a destra nella figura 1.1 e sulla base delle

definizioni di prodotto scalare e modulo abbiamo che l’angolo di bending viene scritto

esplicitamente in termini delle 9 coordinate cartesiane dei tre atomi coinvolti come:

cosα =(x2 − x1)(x3 − x1) + (y2 − y1)(y3 − y1) + (z1 − z2)(z3 − z1)√

(x2 − x1)2 + (y2 − y1)2 + (z1 − z2)2√

(x3 − x1)2 + (y3 − y1)2 + (z3 − z1)2

α = arccos(x2 − x1)(x3 − x1) + (y2 − y1)(y3 − y1) + (z1 − z2)(z3 − z1)√

(x2 − x1)2 + (y2 − y1)2 + (z1 − z2)2√

(x3 − x1)2 + (y3 − y1)2 + (z3 − z1)2

E veniamo ora all’ultimo termine nell’equazione 1.11, ovvero

V tors((r1, r2, ..., rNp , αtors1 , αtors2 , ..., αtorsp ).

Questo termine e dato da:

V tors =1

2

ntors∑k=1

V0i(1 + cos(niϕi − γi)) (1.21)

Nell’equazione 1.21, ϕi e l’angolo diedro per la torsione i-esima e dipende da 4 coordinate

atomiche, mentre V0, ni e γ1 sono dei parametri il cui significato fisico (e valore) dipende

dall’ibridazione degli atomi coinvolti. La sommatoria 1.21 e’ di nuovo estesa a tutte

le possibili Ntors torsioni che si possono contare nella nostra molecola. Tali torsioni

sono circa quattro/cinque volte il numero di atomi e dunque nel caso di Figure 1.1

esse sono circa 60. Per esempio in un sistema tipo etano si contano 9 possibili angoli

diedri, cioe torsioni che sono definite da tutte le possibili quadruplette di atomi che le

definiscono. ϕ e l’angolo tra i piani definiti dagli atomi (punti) 1,2,3 e 2,3,4. Nelle

figure a destra sto guardando le quadruplette lungo l’asse di legame 2-3. ϕ si dice angolo

diedro. Nella molecola di etano: posso definire 9 di queste quadruplette (vedi Figura

14

Figura 1.6: Angolo diedro: angolo fra i piani definito dagli atomi 1,2,3 e 2,3,4 (tratteggiato)

1.7). L’angolo diedro φ dipende dalle coordinate cartesiane dei quattro atomi coinvolti

secondo la Equazione:

cosϕ =(r21 × r23)(r32 × r34)

|r21 × r23||r32r34|(1.22)

dove (r21 × r23) e il prodotto vettoriale del vettore che va dall’atomo 1 all’atomo 2

e del vettore che va dall’atomo 2 all’atomo 3. Il prodotto dei due vettori A e B da

luogo ad un vettore (A×B) perpendicolare sia ad A che a B. Dunque il significato del

Figura 1.7: Possibili torsioni nella molecola di etano.

15

prodotto r21 × r23 e quello della normale al piano 1,2,3. Il senso o verso di A × B si

ottiene applicando la cosiddetta regola della mano destra (vedi Figura 1.8). Nota che in

Figura 1.8: Regola della mano destra per il prodotto vettoriale A×B.

accordo alla definizione 1.22, si ottiene che, nel caso in cui tutti gli atomi 2,3,4 stessero

su uno stesso piano con 1 e 4 in posizione cis, l’angolo diedro corrispondente a questa

configurazione sarebbe pari a zero. Quando 1 e 4 stanno su uno stesso piano ma in

posizione trans l’angolo diedro e pari a 180 (cosϕ = −1). Nella equazione 1.22 per r21

intendiamo il vettore che va da 2 a 1 e cioe r21 = r1−r2 = (x1−x2), (y1−y2), (z1−z2) e

analogamente per gli altri 4 vettori. Il prodotto vettoriale a×b e definito dall’ equazione:

a× b = det

∣∣∣∣∣∣∣∣∣i j k

ax ay az

ax ay az

∣∣∣∣∣∣∣∣∣ = i(aybz − azby) + j(azbx − axbz) + k(axby − aybzx) (1.23)

Dove i ≡ (1, 0, 0), j ≡ (0, 1, 0), k ≡ (0, 0, 1),sono rispettivamente i versori lungo x, y e z

(un versore e un vettore con lunghezza unitaria).

16

Figura 1.9: Configurazioni cis (φ = 0) e trans (φ = 180)

Veniamo adesso al significato fisico dei parametri V0ini e γi nell’ equazione 1.21.

Vediamoli con un esempio: Nel caso delle torsioni che coinvolgano un legame centrale

23 ibridato sp3, abbiamo che ni = 3 e γi = 0 mentre V0i ha un valore nell’ordine delle

frazioni di kcal/mole. V0i e infatti una energia e rappresenta, come vedremo, l’altezza

massima delle barriere torsionali energetiche tra configurazioni eclissata e sfalsata. In

questo specifico caso la torsione relativa e’ data da (cfr Eq. 1.21)’:

V =1

2V0(cos 3y + 1) (1.24)

L’andamento dell’energia 1.24 in funzione dell’angolo diedro e rappresentato in Figura

1.10, dove si riporta a titolo di esempio il potenziale torsionale della molecola di normal-

butano attorno al legame C2-C3. I massimi corrispondono alle configurazioni eclissate

(sfavorite energeticamente) mentre i minimi corrispondono alle configurazioni sfalsate

(favorite energeticamente). L’altezza della barriera torsionale e tale che alle temperature

17

Figura 1.10: Andamento dell’energia torsionale per una torsione di tipo sp3 nella molecola

di n-butano; N.B. I minimi e massimi non sono tutti alla stessa altezza come previsto dalla

Eq. 1.24 per effetto del potenziale non-bonded(vedi sezione successiva) ed tra i gruppi metilici

CH3

ordinarie, essa viene agevolmente superata quando l’ibridazione del legame centrale e

sp3. Nel caso in cui l’ibridazione sia sp2 noi abbiamo che ni = 2, γi = 180 e V0 e’ molto

maggiore rispetto al caso sp3, ed e’ dell’ordine di circa 2 kcal/mol. In questo caso il

potenziale torsionale e’ dato dalla

V (ϕ) =1

2V0[cos(2ϕ− 180) + 1] (1.25)

L’andamento dell’energia in funzione dell’angolo e quella descritta in figura 1.11 In

questo caso V0 e l’altezza della barriera energetica che separa le configurazioni cis e trans.

Sulla base della Eq. 1.25 le configurazioni cis e trans sono iso-energetiche, cioe hanno

la stessa energia. Il Potenziale torsionale (ed in generale il solo potenziale “bonded”

o di valenza) non e in grado di differenziare energeticamente le configurazione cis e le

18

Figura 1.11: Andamento dell’energia torsionale per una torsione di tipo sp2

trans. Come vedremo in seguito sara l’aggiunta del potenziale “non bonded” a quello

di tipo “bonded” che differenziera le due configurazioni cis e trans (o le configurazioni

sfalsate ed eclissate nel caso dell’ibridazione sp3 del legame torsionale). Riassumendo il

potenziale bonded del mio sistema biologico e dato da:

V bonded(r1, r2, ..., rN , αb1, α

b2, ..., α

bρ) =

1

2

Nstrech∑i=1

Ki(ri − roi)2

+1

2

bendings∑1=1

K(b)i (α− α0i)

2 +1

2

n−tors∑k=1

V0i(i+ cos(niϕi − γi)) (1.26)

dove le sommatorie sono estese a tutti i possibili stretching, bending e torsioni in-

dipendentemente dalla loro appartenenza (farmaco, proteina o solvente): in altre parole

19

possiamo sempre scrivere (nel caso dello stretching):

1

2

stretch∑i=1

Ki(ri − r0i)2 =

1

2

stretchp∑i=1

Kpi (rpi − r

p0i)

2 +1

2

stretchf∑i=1

Kfi (rfi − r

f0i)

2

+1

2

stretchs∑i=1

Kfi (rsi − rs0i)2 (1.27)

Dove stretchp, stretchf , stretchs, sono rispettivamente i numeri di legami (o stretching)

presenti sulla proteina, sul farmaco e sul solvente. Cioe il numero totale di stretching nel

sistema si ottiene sommando gli stretching dei suoi componenti (farmaco e proteina e

solvente). Nel caso in cui il solvente sia implicito l’ultimo termine a destra nella Eq. 1.27

non compare. Considerazioni del tutto analoghe possono essere fatte per il potenziale

complessivo di bending o torsionale.

Chiudiamo questa sezione con una considerazione importante: la parte bonded del

potenziale riferito al farmaco e alla proteina non subisce alcuna modifica se si usa o

meno il solvente implicito. La parte che cambia a seconda del metodo usato per trattare

il solvente e quella “non bonded”. Veniamo di seguito al potenziale non bonded nella

equazione 1.8 relativa al potenziale proteina-proteina.

1.2 Potenziale non bonded

L’interazione non bonded tra atomi non direttamente legati o non coinvolti in bendings

impedisce che le catene laterali in una proteina possano andare l’una sull’altra, oppure

e responsabile della formazione dei ponti salini o dei legami ad idrogeno nella struttura

secondaria. In altre parole, se il potenziale bonded mantiene inalterata la struttura

primaria della proteina conservando l’integrita dei legame, degli angoli e delle torsioni (la

cui topologia resta invariata per tutte le possibili conformazioni del sistema), il potenziale

non bonded determina la forma tridimensionale nello spazio di una proteina (molto piu

di quanto non accada in un farmaco che e’ quasi sempre una molecola relativamente

20

rigida). Il potenziale non-bonded determina infatti la struttura secondaria e terziaria

nella proteina. Questa affermazione si riferisce, come vedremo, al potenziale non-bonded

totale e non al solo potenziale non bonded proteina-proteina che compare nella Eq. 1.8

Il bilancio delle interazioni non-bonded solvente-proteina e infatti responsabile della

stabilita delle formazioni terziarie delle macromolecole. Prima di vedere tutto questo

veniamo alle definizioni del potenziale non-bonded per la sola proteina:

Vnon bonded(r1, r2, ..., rNp) =

Np−1∑i=1

Np∑j=i+1

v(ij) |i− j| > 3 (1.28)

Il potenziale non-bonded e dato quindi dalla somma di termini di potenziale a coppia,

che dipendono cioe dalla sola distanza r tra l’i-esimo ed il j-esimo atomo del sistema. In

una proteina ci sono N(N−1)2

possibili coppie atomiche. Per esempio, riferendoci sempre

al caso delle mini proteine di figura Figura 2., abbiamo 22 atomi e dunque 22· 21/2 = 231

possibili coppie di atomi. Se indichiamo gli atomi con dei numeri interi da 1 a 22: Le

possibili interazioni a coppie si ottengono considerando la parte superiore della matrice

di interazione. Nella equazione 1.28 non vengono considerate le coppie che coinvolgono

atomi direttamente legati (stretching) oppure le coppie con atomi in posizione 13 in

un bending. Dunque alle 23 coppie, occorre sottrarre le coppie relative agli stretching

(che sono 21) e quelle relative ai bending (che sono 36) mentre si considerano tra le

interazione non legate le coppie di atomi in posizione 14 nelle torsioni. Dunque le coppie

totali nella equazione 1.28 per il caso della mini-proteina di 1.1 sono 231−21−36 = 174.

L’interazione non-bonded tra atomi legati (1-2) e atomi coinvolti in un bending (1-3)

vengono escluse (vedi Figura 1.12).

Veniamo ora alla forma della funzione v(rij) nella Eq. 1.28 relativa ad una coppia

generica tra le 174 possibili tra atomi non legati nella nostra mini proteina. Tale funzione

e’ data dalla somma di due termini, tale che

Vnon−bonded(r1, r2, ..., rNp) =∑ij

vLJ(rij) + vel(rij) |i− j| > 3 (1.29)

21

Figura 1.12: Matrice delle interazioni a coppia per il potenziale non bonded per la mini-

proteina di Figura 1.1

in cui

vLJ(rij) = 4εij

[(σijrij

)12

−(σijrij

)6]

(1.30)

vel(rij) =qiqj

4πε0rij(1.31)

Nelle Eqs. 1.31, 1.30, 1.29 rij e semplicemente la distanza relativa della coppia di

atomi ij, ovvero: rij = [(xi − xj)2 + (yi − yj)2(zi − zj)2]1/2. rij e’ la quantita che varia

al variare delle disposizioni degli atomi nello spazio tridimensionale. Le quantita σij e εij,

qi e qj sono parametri del potenziale non-bonded, ovvero le quantita fisse αnb1 , αnb2 , ..., α

nbp

nella definizione del force field non bonded:

V non−bonded = V (r1, r2, ..., rN , αnb1 , α

nb2 , ..., α

nbρ ) (1.32)

22

Il termine vLJ(rij) si dice potenziale di Lennard-Jones ed e responsabile delle inter-

azioni idrofobiche. Esso garantisce che le catene laterali non si “compenetrino” ed e per

esempio responsabile della struttura di tipo “stacking” negli aggregati intermolecolari

(vedi il caso del DNA). Il potenziale di Lennard-Jones (equazione 1.30) e dato a sua

volta della somma di due funzioni della distanza rij, una positiva e short-ranged (corto

raggio) e l’altra negativa (vedi Figura 1.13).

Figura 1.13: Potenziale di interazione di Lennard-Jones

A lunghe distanze prevale la parte negativa 4ε(σ/r)6 mentre a corte distanze prevale

la parte positiva 4ε(σ/r)12. Il risultato complessivo della funzione Vij(r) e quello di

un andamento a minimo (la linea solida nella Figura 1.13). Veniamo ora al significato

fisico dei parametri. σij e pari al valore di rij quando l’energia vale 0, quando cioe il

potenziale V(rj) intercetta l’asse delle ascisse. Al minimo, la derivata del potenziale di

23

Lennard-Jones vale 0. Dunque avremo che: [omettiamo per semplicita gli indici ij]

dV (r)

dr= 0 = 4ε0

[σ12

r13MIN

· (−12)− σ6

r7MIN

· (−6)

]0 =

12σ12

r13− 6σ6

r7−→ 12σ12 = 6σ6r6

MIN

r6MIN = 2σ6 −→ rMIN = 2

16σ (1.33)

Dunque le distanza alla quale il potenziale di Lennard-Jones assume il valore minimo

(cioe quelle per le quali la derivata e’ nulla) e pari a rMIN = 216σ. A questa distanza

l’energia vale:

vMIN = v(rMIN) = 4ε0

[σ12

4σ12− σ6

2·σ6

]= 4ε0·

(−1

4

)= −ε0 (1.34)

−ε0 corrisponde dunque al valore dell’energia di Lennard-Jones al minimo. Questi

parametri sono definiti, σ ed ε, su base atomica. In particolare σ e legato al raggio

di Van-der-Waals.

σ12 = σ1 + σ2 (1.35)

La 1.35 e la prima “regola di combinazione ”per il potenziale di Lennard-Jones. La

seconda regola si ottiene dalla:

ε12 = (ε1ε2)1/2 (1.36)

L’energia al minimo e legata alla polarizzabilita degli atomi che legata a sua volta al

numero di elettroni presenti nel sistema. Il termine attrattivo nel potenziale di Lennard-

Jones corrisponde ad una interazione dipolo-dipolo indotto. La nube elettronica dei

due atomi si organizza in media in modo da allineare due dipoli di verso uguale in modo

da minimizzare l’energia elettrostatica di polarizzazione (vedi Figura 1.14). Sulla base

24

Figura 1.14: Origine delle forze attrattive di dispersione. Gli atomi di neon non hanno alcuna

tendenza a formare un legame. I dipoli istantanei presenti sugli atomi producono un campo

elettrico che induce un momento di dipolo sugli atomi vicini, generando delle forze attrattive.

delle regole di combinazione 1.35 1.36 possiamo definire, per il potenziale non-bonded

di Lennard-Jones delle classi “atomiche”. Sulla base della regola di combinazione 1.35, i

parametri di coppia σij ed εij vengono definiti da quantita puramente atomiche tenendo

conto del contesto chimico. Cosı per esempio: σc = 1.6 e εc = 0.2 kcal/mole per un

atomo di carbonio ibridato sp3. I raggi di Van der Waals per gli atomi inseriti in una

proteina dipendono come detto dal contesto chimico cosı come i parametri ε0. Per

esempio un idrogeno legato ad un ossigeno avra un σ piu piccolo di quello relativo ad un

idrogeno legato ad un carbonio alifatico. Questo perche l’H dell’OH ha in media meno

elettroni intorno a se rispetto ad un atomo H legato ad un C (sp2) che e un atomo meno

elettronegativo di O. Nel caso della nostra mini proteina di Figura 1.1 abbiamo una

classe di H legata a C alifatici per quanto riguarda i parametri σ ed ε [in altre parole gli

25

atomi 1, 2, 3, 10, 12, 13, 14, 20, 21 e 22 hanno tutti gli stessi σ ed ε], una classe di atomi

C (sp3) [gli atomi 4, 9, 19 in Figura 1.1], una classe di atomi N [atomi 7, 17], una classe

di atomi C (sp2) [atomi 5 e 15], una classe di idrogeni HN [atomi 8 e 18], una classe di

atomi di ossigeno [atomi 6 e 16], per un totale di 5 classi per i paramatri atomici σ ed

ε. Veniamo ora al potenziale elettrostatico (o coulombiano):

Figura 1.15: Distribuzione di carica elettronica sulla molecola di acqua: la molecola ha attorno

a se 10 elettroni (8 provenienti dall’ossigeno e due dagli idrogeni), ma sull’ossigeno vi sono in

media circa 8.8 elettroni contro gli 0.6 elettroni presenti sugli idrogeni. Cio conferisce all’atomo

di H una carica parziale di 0.4 e, ed all’atomo di O una carica parziale di -0.8 e

Vel(rij) =qiqj

4πε0rij= 332

qiqjrij

(1.37)

Di nuovo rij e la variabile che dipende dalla 6 coordinate cartesiane degli atomi e qi e

qj sono paramatri. ε0 e la permittivita dielettrica nel vuoto che vale circa 8· 10−12 nel

sistema SI (quello senza i coulomb!). qi e qj sono le cariche elettriche parziali presenti

sugli atomi i e j. Tali cariche si originano a causa dalle differenze nelle elettronegativita

degli atomi. Per esempio nella molecola di acqua (vedi Figura 1.15), ci sono in media

piu elettroni sull’ossigeno che non sugli idrogeni: L’ossigeno nell’acqua, a causa della sua

maggiore elettronegativita, ha intorno a se “in media” circa 8.8 elettroni in luogo degli

26

8 dell’atomo isolato, mentre gli idrogeni ne hanno in “media ”0.4 in luogo 1 dell’atomo

isolato. Il numero 332 e’ una costante di conversione, tale che Vel risulta espressa in

Kcal/mol se le cariche parziali sono espresse in elettroni e le distanze in A. A titolo di

esempio, l’energia coulombiana o elettrostatica di “una mole” di “legami ad idrogeno”

O–H e data:

Vel =0.4· 0.8

1.8· 332 ∼ 59 kcal/mole (1.38)

che e in effetti elevatissima. Le energie elettrostatiche sono molto maggiori delle energie

relative al potenziale di Lennard-Jones. Il potenziale Vel(rij) e attrattivo e di segno

negativo se le cariche sono di segno opposto e repulsivo e di segno positivo se le cariche

sono dello stesso segno. Due atomi che si attraggono (ad esempio ossigeno ed idrogeno)

e che non sono legati producono un andamento dell’energia elettrostatica Vel del tipo

indicato in Figura.: Il potenziale di Lennard-Jones al minimo vale pochissimo rispetto

all’energia elettrostatica (vedi Figura 1.16). Tuttavia a corte distanze e il potenziale

repulsivo di Lennard-Jones che prevale. Il minimo del potenziale somma VLJ + Vel

(linea tratteggiata in Figura 1.16) capita per distanze leggermente inferiori rispetto alla

distanza del minimo del potenziale di Lennard-Jones. Il potenziale non-bonded del tipo

1.29 e il solo termine che contribuisce i termini misti Vpf , vps e Vfs nell’equazione 1.2. La

forma funzionale dell’energia e dunque:

V (r1, r2, ..., rN) = V bondedp (r1, r2, ..., rNp) + V non−bonded

p (r1, r2, ..., rNp) +

+ V bondedf (r1, r2, ..., rNf ) + V non−bonded

p (r1, r2, ..., rNf )

+

Np∑i=1

Nf∑j=1

vpf(rij) + +

Np∑i=1

Ns∑j=1

vps(rij) +

Nf∑i=1

Ns∑j=1

vfs(rij)

+Ns∑i=1

Ns∑j=1

v(rij) +

NMOL(solv)∑i=1

V bondeds (ri1, r

i2, ...rNs) (1.39)

I primi quattro termini si riferiscono all’energia bonded e non-bonded del farmaco e

proteina. Il quinto termine e l’interazione non-bonded farmaco-proteina (costruita at-

27

2.4 2.8 3.2 3.6 4 4.4 4.8 5.2 5.6 6Distance (Angstroem)

-40

-20

0

20

40

Ener

gy (

kca

l m

ol-1

)

Coulomb Energy

Lennard-Jones Energy

Total Energy

Figura 1.16: Energia Coulombiana ed Energia di dispersione-repulsione per una interazione

O...H

traverso l’assegnazione delle classi atomiche e le regole di combinazione 1.35 e 1.36 per

la parte Lennard-Jones). Il sesto termine e l’interazione non-bonded proteina-solvente

(se esplicito). Il settimo termine e l’interazione non-bonded farmaco-solvente (se esplic-

ito). L’ottavo termine e l’interazione non-bonded solvente-solvente (sempre costruita

con le regole di combinazione atomo-atomo 1.35 e 1.36) e il nono termine e l’interazione

bonded del solvente del tipo: (nel caso in cui questo fosse acqua).

V bondedsolv =

NMOL∑i=1

1

2KOH(ri,1 − r0)2 +

1

2KOH(ri,2 − r0)2

+ +

NMOL∑i=1

KHOH(αi − α0)2 (1.40)

28

Dove r0 = r0OH = 1A e α0 = αHOH = 105 sono rispettivamente la distanza in-

tramolecolare OH e l’angolo di bending di equilibrio e con ri,1, ri,2 si indicano le distanze

effettive delle lunghezze dei due legami O-H nella molecola i-esima.

Una volta assegnati tutti i paramatri del force field, ovvero stabiliti i valori delle

costanti k, r0, α0, V t0 , µ, γ, σ, ε e q per tutte le classi presenti nel sistema, siamo in

grado dalla espressione generale 4.62 di calcolare il valore dell’energia del sistema per

qualunque configurazione degli atomi che lo compongono, ovvero per qualunque set di

3N (con N = Np +Nf +Ns) coordinate cartesiane che definiscono il sistema complessivo

proteina, farmaco, solvente.

29

Capitolo 2

Elementi di Meccanica Statistica

La meccanica statistica costituisce la base microscopica della termodinamica chimica.

Le proprieta macroscopiche dei sistemi termodinamici. come temperatura, pressione,

densita etc. etc, hanno tutte una origine che discende dalla loro natura molecolare.

Per esempio la temperatura e’ direttamente connessa alla velocita media degli atomi

costituenti le molecole, la pressione e’ collegata alla forza media esercitata dagli urti

continui delle molecole sulle pareti del recipiente contenente il sistema, la densita e’

legata all’impaccamento molecolare e alla natura delle forze tra le molecole costituenti

il sistema massivo.

A partire dalla conoscenza microscopica/molecolare e’ possibile, con gli strumenti

della meccanica statistica, ricostruire e spiegare con rigore il comportamento macroscop-

ico dei sistemi. La termodinamica classica e’ una scienza puramente fenomenologica.

Dati alcuni postulati o principi, frutto dell’osservazione sperimentale, la termodinamica

costruisce una serie di relazioni che permettono di spiegare nel dettaglio il comportamen-

to fenomenologico di sistemi come gas, liquidi, solidi, soluzioni, sistemi multicomponenti,

plurifase e cosı via. I fondamenti basilari della termodinamica classica sono: il concetto

di equilibrio, il primo principio ed il secondo principio.

Il concetto macroscopico di equilibrio: Il concetto di equilibrio da un punto di vista

30

macroscopico e’ intuitivamente semplice. In uno stato termodinamico di equilibrio non

ci sono flussi (movimenti) di materia o di energia, non ci sono cioe cambiamenti rilevabili

sperimentalmente nel sistema. Tutto appare stabile. Per esempio in un sistema monofase

all’equilibrio non sottoposto ad alcun campo esterno, tutte le grandezze termodinamiche

intensive (pressione, temperatura, densita etc) ed estensive (volume, energia, quantita

di materia.. ) non variano nel tempo e nello spazio. Se per esempio io misurassi con un

densimetro la densita in un campione di acqua in un becker troverei (a meno di errori

sperimentali) lo stesso valore in ogni punto del liquido. Se cosı non fosse, io assisterei

ad un flusso di acqua dalle zone ad alta densita a quelle a bassa densita in modo da

equilibrare il sistema. Stessa cosa per la temperatura. In un liquido all’equilibrio, la

temperatura e’ la stessa in qualunque punto del campione io scelga di misurarla. Se cosı

non fosse, si assisterebbe ad un flusso di energia (calore) dalla parte piu calda a quella piu

fredda in modo da equilibrare ovunque nel campione la temperatura. Allo stesso modo

il volume e la quantita di acqua in un campione all’equilibrio appaiono stazionari nel

tempo. Da questi concetti noi deriviamo l’assunto che i sistemi, quando sono perturbati

fuori dal loro stato di equilibrio attraverso un azione esterna di durata finita, una volta

cessata la perturbazione, tendono a ripristinare la condizione di equilibrio, che puo essere

la stessa originaria o una diversa condizione, se le condizioni termodinamiche al contorno

sono state modificate per effetto della perturbazione. Per esempio, se noi modifichiamo

la temperatura della stanza nella quale si trova il nostro campione di acqua portandola

da 20 a 30 gradi centigradi, i valori di equilibrio della densita e della temperatura

nel campione saranno modificati, cosı come saranno modificati i valori del volume e la

quantita di acqua nel becker. Primo principio e secondo della termodinamica: Il primo

principio della termodinamica e’ espresso dall’equazione

dE = δQ− δW (2.1)

Dove δQ rappresenta il calore assorbito dal sistema e δW il lavoro fatto dal sistema ed E

31

e’ l’energia totale del sistema. Dunque il contenuto di energia di un sistema aumenta se si

fornisce calore al sistema oppure se si fa del lavoro sul sistema, viceversa esso diminuisce

se il calore fluisce dal sistema all’ambiente e se il sistema fa del lavoro sull’ambiente. Per

esempio in una espansione adiabatica (cioe senza possibilita che il sistema scambi calore

con l’ambiente), il gas si raffredda perche e’ il sistema che compie del lavoro sull’ambiente

esterno. Per effetto del primo principio (Eq. 2.1), la quantita di energia termica “persa”

dal gas nell’espansione corrisponde esattamente al lavoro di espansione fatto contro la

pressione esterna. Se io invece scaldassi il mio sistema a volume costante, la sua energia

salirebbe per effetto del calore trasferito dall’ambiente. Dal primo principio discende il

fatto che in un sistema isolato (cioe un sistema che non scambia ne’ calore ne’ lavoro con

l’ambiente esterno) l’energia e’ costante. L’energia da un punto di vita termodinamico

e’ una grandezza collegata alla capacita del sistema di fornire all’ambiente o ad altri

sistemi calore e lavoro. Come vedremo, la definizione meccanico-statistica di energia e’

collegata strettamente alla natura molecolare del sistema. Nella Eq. 2.1 i differenziali

δQ e δW relative a calore e lavoro sono differenziali non esatti, cioe riferiti a quantita

il cui valore finale in un dato processo dipende dal modo in cui il processo stesso e’

effettuato. Q e W non sono funzioni di stato, grandezze cioe che dipendono solo dalle

variabili termodinamiche relative agli stati finali ed iniziali e non dal particolare percorso

seguito per passare da uno stato all’altro. Per un sistema che e’ in grado solo di fare

lavoro attraverso variazioni di volume (come per esempio un gas in un pistone), la Eq.

2.1 puo essere scritta in termini di funzioni di stato come

dE = TdS − pdV (2.2)

nella quale sia S (entropia) che V (Volume) sono funzioni di stato e T e p sono la

temperatura e la pressione del sistema. La fissazione (incomprensibile alla maggior parte

studenti) della termodinamica sui “gas nei pistoni” deriva dalla origine ingegneristica

di questa scienza, che venne applicata per la messa a punto e lo studio dei primi motori

a scoppio, per i quali il modello del “gas nel pistone” era fondamentale. La Eq. 2.2,

32

con l’introduzione della funzione di stato S ci porta direttamente alla definizione del

secondo principio della termodinamica.

Il secondo principio della termodinamica e’ assai piu sfuggente e meno chiaro. Esso

ha a che fare con la direzione del tempo. Se si butta una pallina per terra, questa

rimbalzera con rimbalzi sempre meno alti fino a fermarsi a terra. La terra non trasferisce

mai energia alla pallina in modo che questa faccia rimbalzi sempre piu alti. Cosı il

calore non passera mai spontaneamente da un corpo freddo ad uno piu caldo senza

che non intervenga qualcos’altro che lo faccia succedere, per esempio che si fornisca

energia elettrica [lavoro] ad un frigorifero (definizione di Clausius). In termini di motori

(la fissazione della termodinamica classica) e’ impossibile produrre effetti meccanici

semplicemente raffreddando un corpo ad una temperatura inferiore a quella dell’oggetto

piu freddo circostante (definizione di Kelvin). E’ relativamente semplice capire che se

fosse possibile prendere calore da un colpo piu caldo e produrre una quantita equivalente

di lavoro (senza cederne una parte ad un corpo piu freddo), cio equivarrebbe a far passare

calore da un corpo piu freddo ad uno piu caldo se la macchina di Kelvin fosse accoppiata

ad un frigorifero. Se infatti il lavoro prodotto dalla macchina di Kelvin (quella cioe che

viola i secondo principio) fosse proprio quello che il frigorifero prende dal corpo freddo

per trasferirlo fuori da se, l’effetto netto sarebbe quello che il frigorifero si raffredda a

scapito dell’ambiente caldo circostante. Questi discorsi sfuggenti implicano l’esistenza di

una funzione di stato che ci dica in quale senso procede un dato processo termodinamico

“non reversibile’ o “spontaneo”. Questa funzione di stato la cui evoluzione segue quella

del tempo, e’ detta entropia. In un sistema isolato, c’e un senso univoco del tempo tale

che

∆S ≥ 0 (2.3)

Cioe in un sistema isolato avvengono processi o che non fanno variare questa funzione di

stato (processi reversibili) o che la fanno variare facendola aumentare (processi spontanei

o irreversibili). Per esempio, nel sistema di Figura 2.1, che rappresenta due gas di tipo

33

Figura 2.1: Nel sistema si e’ tolta e rimesse per due volte una parete impenetrabile. Durante

il tempo in cui il sistema e’ stato senza parete i gas si sono mescolati. Non puo MAI accadere,

sulla base del secondo principio, che a partire da gas mescolati, togliendo la parete i gas

spontaneamente si smescolino.

diverso, io ho tolto e rimesso in due eventi successivi una parete impenetrabile ai gas.

Tra i due eventi, i due gas avevano a disposizione l’intero recipiente. La nostra esperienza

ci dice che il tempo deve essere scorso da sinistra verso destra. Corrispondentemente,

34

la grandezza entropia deve essere cresciuta andando da sinistra verso destra.

Rivediamo ora questi concetti - equilibrio termodinamico, primo principio, secondo

principio - da un punto di vista microscopico.

2.1 Equilibrio termodinamico, funzione di distribuzione

e spazio delle fasi

Cosa succede se io in un corpo (un gas o liquido per esempio) all’equilibrio mi facessi

piccolo come una delle molecole che compone il sistema? Cosa vedrei intorno a me?

Vedrei una grandissima confusione, molecole che arrivano in ogni istante da tutte le

possibili direzioni, con velocita casuale e sempre diverse. Avrei cioe l’impressione di

un ambiente in moto turbinoso e disordinato, il contrario assoluto dello stato di quiete

constatato da un osservatore esterno di dimensioni normali. Cos’e che microscopica-

mente si trova all’equilibro? Cio che non varia nel tempo, per me osservatore piccolo,

e’ la distribuzione delle diversita che vedo a attorno a me. Se per esempio io contassi

quante molecole arrivano dalla mia destra in un dato lasso di tempo e se ripetessi questa

stessa misura molte volte troverei che questo numero non cambia nel tempo, o meglio

cambia poco attorno ad un valore medio costante. Non solo. Ma se ripetessi questa

misura contando stavolta le molecole che arrivano dall’alto troverei un numero medio

identico a quello trovato per le molecole che arrivavano da sinistra e cosı per ogni pos-

sibile direzione. Oppure ancora se misurassi in un certo istante le velocita di tutte N

le particelle e contassi quante molecole hanno velocita comprese tra 1 e 2 m/s, quante

tra 2 e 3 m/s, quante tra 3 e 4 m/s etc. etc.( e calcolassi cosı le frazioni n12/N , n23/N ,

n34/N ,..) e se ripetessi questa stessa misura in un qualsiasi istante successivo, troverei

che la frazione di molecole aventi velocita tra 1 e 2 m/s o tra 2 e 3 m/s o tra 3 e 4

m/s, non sono variate apprezzabilmente rispetto a prima. In un sistema all’equilibrio

termodinamico, la distribuzione delle diversita microscopiche e’ stazionaria nel tempo.

35

Nell’esempio sopra la frazione di molecole aventi velocita tra x ed x+dx m/s rappresenta

una probabilita:1

Probabilita dell′evento X =numero di eventi di tipo X

numero di tutti eventi possibili. (2.4)

In un sistema termodinamico fatto da N particelle, gli eventi possibili sono un numero

incredibilmente alto. Tali eventi si definiscono stati microscopici. Un particolare sta-

to microscopico “ν” del sistema e’ caratterizzato dalla specificazione di tutte le 3N

coordinate r1, r2, ..., rN atomiche e di tutti i 3N momenti p1,p2, ...,pN delle par-

ticelle (atomi) che costituiscono il sistema. La funzione di distribuzione delle diversita

non e’ altro che la probabilita di osservare il sistema in un certo stato microscopico ν

caratterizzato dalle 6N coordinate r1, r2, ..., rN ,p1,p2, ...,pN, ovvero:

Pν = f(r1, r2, ..., rN ,p1,p2, ...,pN) dr1dr2, ...drNdp1dp2, ...dpN =

= Probabilita′ che la particella 1 stia tra r1 e r1 + dr1

e che la particella 2 stia tra r2 e r2 + dr2

...

e che particella 1 abbia momento compreso tra p1 e p1 + dp1

e che la particella 2 abbia momento compreso tra p2 e p2 + dp2

..... (2.5)

Lo spazio degli eventi microscopici possibili e’ 6N-dimensionale ed e’ chiamato anche

spazio delle fasi. Un punto in questo spazio 6N-dimensionale rappresenta un particolare

1Per esempio quando si tirano due dadi (uno rosso e uno blu), gli eventi possibili sono 36 (ciascuna

faccia del dado rosso si combina con le sei facce del dado blu e viceversa). La probabilita che si verifichi

la somma “7 “ e’ data dagli eventi possibili (che sono 6, dado rosso tre e dado blu quattro e viceversa,

dado rosso 2 e dado blu cinque e viceversa, dado rosso uno e dado blu sei e viceversa. Dunque la

probabilita di “fare 7”’ lanciando due dadi e’ 6/36 = 1/6

36

stato microscopico del sistema caratterizzato da 6N coordinate r2, ..., rN ,p1,p2, ...,pN.

Il dominio delle spazio delle fasi e’ costituito, per ciascuna coordinata ri, dal volume

che delimita il sistema, mentre ciacuna delle componenti del momento pi puo’ variare

tra −∞ e ∞. Il sistema evolve nel tempo attraversando molti stati microscopici (le

coordinate e le velocita delle particelle cambiano disordinatamente nel corso del tempo).

Questa evoluzione temporale corrisponde ad una traiettoria nello spazio delle fasi 6N-

dimensionale.

In generale, la funzione probabilita Pν rappresenta come detto la frazione dell’evento

ν sul totale degli eventi. Per come e’ definita, essa e’ necessariamente sempre positiva.

Pν > 0 (2.6)

Inoltre la probabilita di trovare il sistema in uno qualunque degli stati microscopici

possibili si ottiene sommando le probabilita individuali di tutti gli stati microscopici.

Questa probabilita cumulativa deve ovviamente essere unitaria, cioe si deve avere la

certezza di osservare il sistema in uno dei ν stati possibili. Deve essere cioe∑ν

Pν = 1 (2.7)

La condizione Eq. 2.7 si dice condizione di normalizzazione e deve essere obbedita

dalla funzione di distribuzione delle diversita Eq. 2.5 perche questa rappresenti una

probabilita.

La funzione di distribuzione di Eq. 2.5 e’ strettamente connessa alla natura micro-

scopica del sistema e la sua conoscenza permette il calcolo di proprieta medie del sistema.

Per esempio, dato un certo stato microscopico ν caratterizzato dalle coordinate e mo-

menti r1, r2, ..., rN ,p1,p2, ...,pN possiamo valutare l’energia Eν di questo particolare

stato che (come abbiamo visto nel capitolo precedente ) e’ una funzione delle coordinate

del sistema. Il valor medio (cioe quello che si osserva sperimentalmente) dell’energia puo

essere calcolato pesando l’energia degli stati microscopici ν con la funzione Pν in modo

37

tale che contribuiscano maggiormente alla media gli stati microscopici piu probabili.

Cioe

E(media) = E =< E >=∑

PνEν (2.8)

In generale avremo che, data una proprieta Aν associabile ad uno stato microscopico

ν, la media della grandezza < A > che si osserva sperimentalmente e’ data da

< A >=∑

PνAν (2.9)

La Eq. 2.9 rappresenta il legame tra il livello microscopico e’ l’osservazione sperimentale

macroscopica. In linea di principio, a partire dalla conoscenza molecolare (e dunque dalla

conoscenza di Pν) e’ possibile ricavare grandezze sperimentali macroscopiche, come per

esempio (e come vedremo) la costante di affinita farmaco recettore. Avendo definito la

funzione di distribuzione Pν degli stati microscopici (delle diversita) nella Eq. 2.5, sulla

base della Eq. 2.9, vediamo che se la probabilita Pv e’ stazionaria (cioe non varia nel

tempo come accade per le facce di un dado) allora tutte le media macroscopiche sono a

loro volta stazionarie: il sistema macroscopicamente si trova in uno stato di quiete o di

equilibrio.

2.2 Formulazione microscopica di primo e secondo

principio

Per quanto riguarda il principio 2.1, il suo equivalente microscopico e’ presto detto:

l’energia totale, come somma dell’energia potenziale espressa da una formulazione atom-

istica del tipo 1.2 e dell’energia cinetica, in un sistema isolato, e’ costante. In termini

matematici si ha che

N∑i

p2i

2mi

+ V (r1, r2, ..., rN) = E = costante (2.10)

38

dove pi ≡ mivi rappresenta il momento della quantita di moto dell’atomo i-esimo.

La sommatoria a sinistra nella Eq. 2.10 e’ l’energia cinetica del sistema. Ogni stato

microscopico e’ caratterizzato da una sua energia cinetica. La media <∑N

ip2i

2mi>,

come vedremo, e’ collegata alla temperatura del sistema.

La formulazione microscopica del secondo principio e’ forse piu comprensibile delle

definizioni macroscopiche basate sulle macchine termiche. Come abbiamo detto sopra,

nei sistemi all’equilibrio esiste una funzione stazionaria delle diversita microscopiche

del sistema, Pv che dipende dalle caratteristiche delle singole particelle quali velocita,

posizione, energia, oltre che dalle condizioni termodinamiche in cui si trova il sistema. In

generale, i sistemi evolvono spontaneamente verso l’equilibrio in modo da rendere questa

funzione delle diversita piu distribuita possibile. In altre parole tale funzione evolve in

modo tale da diminuire la nostra conoscenza del sistema.2 Per capire meglio questa

affermazione, facciamo un esempio. Nello stato 1, abbiamo un gas che e’ costretto a stare

in una meta’ di un recipiente. Nello stato 1 la funzione di distribuzione delle diversita

ci avrebbe detto che nella meta’ di destra del recipiente non vi era alcuna probabilita

di trovare particelle. Quando si toglie la barriera, la funzione di distribuzione delle

diversita comincia ad evolvere (sappiamo che evolve in modo da aumentare la nostra

ignoranza circa lo stato microscopico del sistema) e la probabilita di trovare particelle

2Questa affermazione deve intendersi riferita ad i sistemi isolati; per esempio l’acqua sottoraffreddata

contenuta in un bicchiere in contatto termico con il laboratorio, alla temperatura di -1 spontaneamente

diventa solida. E’ difficile sostenere che in questo processo la nostra conoscenza sia diminuita (o il dis-

ordine aumentato). In un solido la nostra conoscenza delle posizioni degli atomi e’ molto maggiore di

quella relativa al liquido. Il punto e’ che il processo non e’ avvenuto in un sistema isolato e la nostra

conoscenza del sistema e’ localizzata al solo recipiente contenente il liquido che e’ cristallizzato. Se

la nostra conoscenza locale del liquido che e’ diventato solido e’ sicuramente aumentata, e’ diminuito

certamente il nostro grado di conoscenza dell’ambiente circostante. Verso quell’ambiente infatti e’ stato

trasferito calore, per effetto della solidificazione, calore che ha reso maggiormente distribuita la fun-

zione delle diversita relativa all’ambiente esterno. La funzione di distribuzione complessiva del sistema

recipiente +ambiente esterno in accordo al secondo principio deve essere sicuramente aumentata.

39

nella parte destra diviene diversa da zero, fino a quando, nello stato di equilibrio 2,

essa sara identica alla probabilita di trovare particelle nella parte sinistra del recipiente.

Mentre noi nello stato 1 con certezza assoluta avremmo potuto affermare che non c’erano

atomi nella parte destra, nello stato 2 questa conoscenza e’ andata persa e non possiamo

dire altro che le molecole possono stare sia di qua che di la, con uguale probabilita. Nella

Figura 2.1 vediamo un altro esempio di diminuzione di conoscenza per effetto del secondo

principio. Prima di rimuovere la parete impenetrabile noi sapevamo che gli atomi rossi

stavano con certezza nella parte superiore del recipiente. Dopo la rimozione della parete,

questa certezza/conoscenza non c’e piu visto che la sola cosa che possiamo dire e’ che gli

atomi rossi, cosı come quelli blu possono essere trovati con eguale probabilita nella parte

superiore ed inferiore. E’ chiaro che l’entropia deve essere legata in qualche modo alla

funzione di distribuzione delle diversita. Vedremo meglio come nella prossima sezione.

2.3 Insieme microcanonico

L’evoluzione della funzione delle diversita per il sistemi fuori dall’equilibrio, che tende a

diminuire la nostra conoscenza al massimo grado possibile compatibilmente con le con-

dizione esterne (al contorno) di volume, pressione temperatura etc etc, ci suggerisce come

dovrebbe essere fatta la funzione di distribuzione delle diversita per un sistema isolato.

Il sistema isolato in meccanica statistica si chiama ensemble microcanonico. Questo

ensemble e’ caratterizzato da un numero N di particelle costante, un volume V fisso, ed

una energia totale E costante. Queste sono le condizioni termodinamiche al contorno che

caratterizzano l’ensemble microcanonico detto anche ensemble NVE. L’ensemble micro-

canonico costituisce l’insieme fondante della meccanica statistica. Vediamo di seguito

perche.

Se il nostro sistema isolato (ad energia costante E) si trovasse fuori dall’equilibro, la

sua funzione di distribuzione delle diversita Pv non sarebbe stazionaria ed evolverebbe

40

nel tempo verso una distribuzione finale, compatibile con le condizioni macroscopiche

N V E. Questa funzione di distribuzione delle diversita finale (di equilibrio), in accordo

al secondo principio, sara necessariamente quella che presuppone la conoscenza minima

dello stato microscopico del sistema. Il sistema, nello stato stazionario finale come in

ogni stato intermedio, puo assumere un numero enorme di stati possibili tutti con la

stessa energia E. Indichiamo questo numero astronomico di stati possibili con Ω. Nella

Eq. 2.10 ci sono infatti infiniti modi di combinare le posizioni e le velocita degli atomi

in modo da dare sempre lo stesso numero E corrispondente all’energia macroscopica

(costante) del sistema. La nostra ignoranza circa lo stato microscopico del sistema e’

massima quando le probabilita di osservare il sistema in uno qualsiasi di questi stati

microscopici - tutti ad energia totale E (cinetica piu potenziale) - sono tutte identiche.

Sulla base del secondo principio (o di massima ignoranza), questa (Pv tutte uguali ∀ν)

deve necessariamente essere la forma finale (di equilibrio) della funzione di distribuzione

delle diversita. Le probabilita di osservare un dato stato microscopico nel sistema NVE

sono pertanto tutte identiche e pari

Pν =1

Ω(2.11)

La quantita Ω e’ un fattore di normalizzazione che corrisponde come detto al numero

totale di stati microscopici (combinazioni di velocita e posizioni) del sistema microcanon-

ico tutti ad energia E (stati degeneri). Questi stati degeneri, in accordo alla Eq. 2.11,

sono tutti equiprobabili3. La Pv e’ la funzione di distribuzione delle diversita dell’insieme

microcanonico. Essa e’ detta semplicemente funzione di distribuzione microcanonica e

rappresenta la probabilita di trovare un certo stato microscopico ν. La probabilita di

3nell’esempio di nota 1 dell’acqua nel bicchiere che solidifica, la distribuzione finale nel sistema isolato

(NVE) bicchiere + ambiente esterno non e’ chiaramente quella di Eq. 2.11. Perche? Semplicemente

perche non si e’ atteso abbastanza. Il sistema non e’ in equilibrio prefetto ma in un stato di quasi

equilibrio o stato metastabile. Nello stato finale del sistema isolato mondo+bicchiere la distribuzione

degli stati microscopici sara proprio quella di 2.11 ed il bicchiere con l’acqua dentro sara sparito.

41

trovare il sistema NVE in uno qualunque degli Ω stati possibili si ottiene sommando le

probabilita individuali di tutti gli stati (che sono tutte pari ad 1/Ω). Questa probabilita

deve ovviamente essere unitaria, cioe si deve avere la certezza di osservare il sistema in

uno dei ν stati possibili. Deve essere cioe

∑v

Pv =Ω∑v=1

1

Ω=

Ω

Ω= 1 (2.12)

In meccanica statistica i fattori di normalizzazione delle funzioni di distribuzione sono

sempre legati a proprieta termodinamiche. Nel caso di Ω, tale quantita deve chiaramente

avere a che fare con l’entropia. In un sistema isolato, l’evoluzione del sistema e’ tale

che esso cerca di massimizzare i modi possibili di avere la stessa energia E (vedi ad

esempio l’esempio del gas che si espande nella parte vuota di un recipiente adiabatico).

Questa massimizzazione della degenerazione degli stati microscopici deve corrispondere

alla riduzione della nostra conoscenza con massimizzazione dell’entropia in accordo alla

seconda legge 2.3. Abbiamo quindi che

S = kB ln Ω (2.13)

dove kB e’ detta costante di Boltzmann ed e’ data dalla costante dei gas perfetti R divisa

il numero di Avogadro NA. La definizione 2.13, sfruttando dalla 2.2 che (∂E/∂S)V = T

implica che (∂kB ln Ω

∂E

)V

=1

T−→

(∂ ln Ω

∂E

)V

= β =1

kBT(2.14)

La Eq. 2.14 costituisce una definizione meccanico-statistica di temperatura. Essa

ci dice infatti che la temperatura e’ legata al tasso di variazione del logaritmo della

degenerazione (= numero di stati microscopici possibili) rispetto all’energia. Siccome

quando io aumento l’energia di un sistema (e poi lo isolo nuovamente), io aumento sicu-

ramente la sua degenerazione, ne consegue che anche il logaritmo (che e’ una funzione

monotona) aumenta con E e dunque la quantita β deve essere necessariamente positiva.

Non si danno cioe temperature negative. La Eq. termodinamico-statistica 2.13 e’ in

42

accordo con la proprieta di estensivita’ della grandezza macroscopica di stato nota come

entropia. Consideriamo due sistemi microcanonici A e B tali che il primo sia caratter-

izzato da NA, EA, VA ed il secondo da NB, EB, VB. Chiaramente il sistema somma sara

pure microcanonico ad energia EA+B = EA + EB con un numero di particelle pari a

NA+B = NA +NB ed un volume VA+B = VA + VB. L’entropia macroscopica del sistema

somma sara analogamente data da SA+B = SA + SB. D’altro canto se ΩA e’ il numero

di stati microscopici del sistema microcanonico A ed ΩB e’ il numero di stati micro-

scopici del sistema microcanonico B, il numero di stati totali del sistema microcanonico

somma sara dato da ΩA+B = ΩA × ΩB, poiche ogni microstato di B si combina con

tutti i possibili ΩA stati di A e viceversa. Per le proprieta dei logaritmi, abbiamo che

kB ln ΩA+B = kB ln ΩAΩB = kB ln ΩA + kB ln ΩB = SA + SB.

2.4 Insieme Canonico

L’insieme microcanonico e’ l’insieme fondante della meccanica statistica poiche la sua

funzione di distribuzione Eq. 2.11 e’ nota a priori. Tuttavia nell’insieme microcanonico

non c’e alcun controllo delle variabili termodinamiche intensive. I sistemi che interessano

di piu dal punto di vista sperimentale sono quelli dove noi possiamo imporre dall’esterno

il valore di grandezze intensive come temperatura e pressione o potenziale chimico. Per

esempio una cellula vivente e’ un sistema a temperatura (e non energia) costante. L’in-

sieme a temperatura, volume e numero di particelle costanti, e’ detto insieme canonico

o insieme NVT. La funzione di distribuzione degli stati microscopici nell’insieme canon-

ico puo’ essere derivata da quella microcanonica, considerato un insieme canonico in

contatto termico (cioe che puo’ scambiare solo energia) con un enorme insieme isolato

circostante che funge da “bagno” termico. L’insieme che contiene il sistema canonico

(vedi Figura 2.2) ha il compito di fornire o togliere calore al sistema canonico NVT

inscritto in modo da mantenere la temperatura di questo costante. Nell’insieme canon-

43

Figura 2.2: Sistema canonico in contatto termico con un bagno. Il sistema

ico l’energia non e’ costante come nell’insieme somma, ma fluttua attorno ad un valore

medio costante (deve essere cosı, altrimenti il sistema canonico non si troverebbe all’e-

quilibrio). La probabilita di osservare un certo stato microscopico nell’insieme canonico

con energia Ev deve essere proporzionale al numero di stati microscopici del bagno com-

patibili con questo valore. Se l’insieme canonico ha in un certo istante energia Ev,

significa che il bagno in quello stesso istante ha energia E − Ev poiche l’energia del

sistema somma deve essere sempre pari ad E. Se il bagno ha energia E − Ev, il nu-

mero di microstati (equiprobabili) del bagno corrispondenti a questo valore sara dato

da Ω(E − Ev). Maggiore sara questo numero di microstati del bagno compatibili con

il valore Ev del sistema canonico e maggiore sara la probabilita di osservare l’evento

energia Ev nell’insieme canonico. Siccome gli stati del bagno sono tutti equiprobabili,

44

deve necessariamente essere

Pv ∝ Ω(E − Ev). (2.15)

Considerando che E Ev (il bagno e’ enorme rispetto al sistema NVT inscritto),

abbiamo che

Pv ∝ Ω(E − Ev) = e(ln Ω(E−Ev))

= eln Ω(E)−(∂ ln(E)/∂E)Ev

∝ e−βEv (2.16)

Nella Eq. 2.16, si e’ prima sviluppato in serie di Taylor la funzione ln(E −Ev) attorno

a valore E e poi si e’ utilizzata la Eq. 2.14. La Pv e’ la funzione di distribuzione (delle

diversita) nell’insieme canonico a meno di un fattore 1/Q di normalizzazione. Questo

fattore si ottiene imponendo la condizione che la probabilita di osservare uno qualunque

degli stati microscopici ν dell’insieme NVT debba essere unitaria, ovvero

1 =∑v

Pv = (1/Q)∑v

e−βEv (2.17)

da cui otteniamo che

Q =∑v

eβEv ≡ QNV T (2.18)

Lo stato energetico del sistema Eν dipende dallo stato microscopico del sistema in

accordo alla equazione

Eν = H(r1, r2, ..., rN ,p1,p2, ...,pN) =N∑i

p2i

mi

+ V (r1, r2, ..., rN) (2.19)

La funzione delle 6N coordinate H(r1, r2, ..., rN ,p1,p2, ...,pN), corrispondente all’ener-

gia totale istantanea del sistema, e’ nota come Hamiltoniano. La probabilita canonica

Pν puo’ dunque essere scritta esplicitando la dipendenza dallo stato microscopico del

sistema in termini di Hamilitoniano come

Pν =e−βH(r1,r2,...,rN ,p1,p2,...,pN )dr1dr2, ...drNdp1dp2, ...dpN∫e−βH(r1,r2,...,rN ,p1,p2,...,pN )dr1dr2, ...drNdp1dp2, ...dpN

(2.20)

45

La quantita QNV T =∫e−βH(r1,r2,...,rN ,p1,p2,...,pN )dr1dr2, ...drNdp1dp2, ...dpN (che nor-

malizza la probabilita canonica) si chiama funzione di partizione del sistema e, analoga-

mente al caso dell’ensemble microcanonico, essa deve essere collegata ad una qualche

grandezza termodinamica caratteristica del sistema termostatato. Nei sistemi termod-

inamici a temperatura costante, la grandezza basilare e’ l’energia libera di Hellmotz

A.4 E’ relativamente semplice fare vedere che il logaritmo della funzione di partizione

canonica QNV T =∑

ν e−βEν e’ uguale a −βA. Infatti dalla termodinamica classica noi

sappiamo che

A = E − TS (2.21)(∂A

∂T

)N,V

= −S (2.22)

Dove E, S, T sono l’energia e l’entropia e la temperatura del sistema. Dalla Eq. 2.22

e dalla definizione Eq. 2.14, noi possiamo scrivere che(∂A

∂T

)N,V

= −(∂A

∂β

)N,V

β

T= −S(

∂A

∂β

)N,V

β = TS

A+

(∂A

∂β

)N,V

β = A+ TS(∂(βA)

∂β

)N,V

= E (2.23)

4In un sistema termostatato, l’equilibrio termodinamico e’ raggiunto quando l’energia libera del

sistema raggiunge il valore minimo. E’ possibile mostrare che la minimizzazione dell’energia libera

A del solo sistema termostato corrisponde alla massimizzazione dell’entropia del sistema + ambiente

esterno.

46

Ma dalla definizione di media 2.8, con Pν = e−βEν/QNV T abbiamo che

E =< E > =

∑ν Eνe

−βEν∑ν e−βEv

= −(∂ ln(

∑ν e−βEv)

∂β

)= −

(∂ lnQNV T

∂β

)(2.24)

la Eq. 2.24 e’ una tipica relazione meccanico-statistica in cui il valore medio macroscop-

ico di una grandezza termodinamica (l’energia) e’ espresso in funzione di una quantita,

QNV T , che dipende dagli stati microscopici del sistema. Sostituendo il valore di E

trovato nella Eq. 2.24 nella Eq. 2.23, otteniamo finalmente che

−βA = lnQNV T (2.25)

Di nuovo, come visto nel microcanonico, anche nell’insieme canonico il logaritmo della

funzione partizione (cioe del fattore che normalizza la probabilita Pν) e’ collegato a

proprieta macroscopiche dell’ensemble, l’entropia nel caso del microcanonico (NVE) e

l’energia libera nel caso del canonico (NVT). Per il caso che ci interessa, un farmaco ed

un recettore che interagiscono in condizioni fisiologiche (cioe a temperatura costante),

la relazione Eq.2.25 sara la nostra base di partenza.

47

Capitolo 3

Funzione Partizione ed Equilibrio

Chimico

Nel capito precedente abbiamo visto la funzione di partizione del sistema canonico, a N,

V, T costanti, Eq. 2.18 Abbiamo inoltre visto che tale funzione e’ collegata all’energia

libera di Hellmotz (vedi Eq. 2.25). In questa sezione vedremo che nel caso di un gas

di particelle che non interagiscono, la funzione di partizione canonica puo essere calco-

lata esattamente. Cio consente, come vedremo, di valutare la costante di un equilibrio

chimico sulla base delle proprieta molecolari di reagenti e prodotti di reazione.

3.1 Gas perfetti

La Eq. 2.18 viene scritta classicamente come un integrale sullo spazio degli stati

microscopici o spazio delle fasi:

Q =1

h3N· 1

N !

∫e−βH(rN !pN )drNdpN (3.1)

Nella Eq. 3.1 abbiamo utilizzato la notazione abbreviata per indicare le coordinate

48

delle particelle

rN ≡ r1, r2, ....rN

pN ≡ p1,p2, ....pN

Il fattore h3N viene introdotto per rendere la funzione di partizione Q classica adi-

mensionale. h e la costante di Planck ed ha le dimensioni di una energia per tempo

[ricordare che E = hν]. L’integrale nella 3.1 ha le dimensioni di un momento elevato

alla 3N moltiplicato per una distanza alla 3N , ovvero, per l’appunto di una energia

per tempo elevata alla 3N . Dunque il drNdpN ' h3N rappresenterebbe il “volume”

minimo nello spazio delle fasi in accordo al principio di indeterminazione secondo cui

∆x∆px = h. Il fattore 1N !

tiene conto della indistinguibilita delle particelle costituenti il

sistema.1 La funzione H(rN ,pN) e l’Hamiltoniano del sistema (vedi anche Eq. 4.72) :

H(rNpN) = K(pN) + V (rN) (3.2)

L’Hamiltoniano in coordinate cartesiane e dato dalla somma dell’energia cinetica K(pN)

che dipende solo dai momenti delle particelle che costituiscono il sistema e da V (rN)

che e il potenziale (force-field) tra le particelle che costituiscono il sistema. Il potenziale

interatomico e’ stato ampiamente discusso nel capitolo 1 e la forma funzionale di V (rN)

e’ riportata nella Eqs. 1.2 e 4.62.

Sfruttando la relazione 3.2 la funzione partizione classica, Eq. 3.1, puo essere scritta

come il prodotto di due fattori:

Q =1

h3N

1

N !

∫e−βK(pN )dpN

∫e−βV (rN )drN (3.3)

Q = Qgas·Qexcess (3.4)

1L’introduzione del fattore di volume h3N per rendere la funzione di partizione adimensionale e di

1/N ! (che e’ essenziale per rendere l’energia libera A estensiva come deve) sono manifestazioni della

natura quantistica della materia.

49

dove nella Eq. 3.4 si sono definite le quantita

Qgas =V N

h3NN !

∫e−βK(pN )dpN (3.5)

Qexcess =1

V N

∫e−βV (rN )drN (3.6)

La Qgas e’ detta funzione di partizione del gas ideale e si chiama cosı perche in un gas

perfetto di particelle monoatomiche il potenziale interatomico e nullo ( V (rN) = 0 ) e

dunque:

Qgas =1

h3NN !

∫e−βH(rN !pN )drdpN =

1

h3NN !

∫e−βK(pN )dpNV

N (3.7)

c. v. d.2

La Qexcess = 1V N

∫e−βV (rN )drN e detta funzione di partizione di eccesso. Per un dato

insieme termodinamico, −kBT lnQexcess rappresenta infatti l’energia libera di eccesso

(dovuta alle interazioni intermolecolari) rispetto all’energia libera del gas perfetto:

A = −kBT lnQ = −kBT lnQgasQexcess

= −kBT lnQgas − kBT lnQexcess = Agas + Aexcess

Nel caso in cui si abbia a che fare con un gas perfetto monoatomico (cioe costituito

da atomi), la Eq. 3.5 puo essere integrata direttamente. Infatti noi sappiamo che:

K(pN) =3N∑i=1

p2i

2m

Dove m e la massa degli atomi che costituiscono il gas. la sommatoria e’ estesa a tutte

le 3N componenti dei momenti. Dunque:

2Nel derivare la Eq. 3.7, si e’ sfruttato il fatto che∫V

∫V...∫Vdr1dr2....drN = V N poiche

l’integrazione va fatta sul dominio delle coordinate, ciascuna delle quali varia entro il volume V.

50

Qgas =V N

h3NN !

∫ ∫..

∫exp

(−β

3N∑i=1

p2i

2m

)dpN

=V N

h3NN !

3N∏i=1

∫exp

(−β p

2i

2m

)dpi (3.8)

=V N

h3NN !

(∫exp

(−β p

2

2m

)dp

)3N

(3.9)

=

√(2πmkbT )3N

h3NN !V N =

√(2πmkbT

h2

)3NV N

N !(3.10)

Nella Eq. 3.8, sfruttando le proprieta degli esponenziali, si e’ fattorizzato l’integrale 3N

dimensionale in un prodotto di 3N integrali ciascuno dipendente da una sola compo-

nente del momento. Nel passaggio successivo, si e’ sfruttato il fatto che le componenti

dei momenti sono definite tutte sullo stesso dominio e dunque gli integrali definiti ar-

gomento della produttoria sono tutti identici. Nell’ultimo passaggio infine si utilizzata

la definizione 2.14 per β e sfruttato il fatto che l’integrale gaussiano nella Eq. 3.9 vale∫e−αx

2dx = (π/α)1/2.

3.2 Funzione di partizione molecolare

La 3.10 e la funzione di partizione traslazionale per un gas perfetto monoatomico. Ve-

diamo adesso come e’ fatta la funzione di partizione in un gas di molecole identiche

che non interagiscono tra loro, ovvero in un gas perfetto molecolare. Come vedremo

la conoscenza della funzione di partizione molecolare e’ importante per il calcolo delle

costanti di equilibrio nelle reazioni chimiche. Per un gas perfetto molecolare, l’energia

totale del sistema e’ semplicemente data dalla somma delle energie delle singole molecole

e non sono presenti termini di accoppiamento (interazioni non bonded) di alcun tipo.

Dunque un microstato del sistema e’ tale che

Eν = Eν1 + Eν2 + Eν3 + ...EνN (3.11)

51

in cui ν1, ν2, ... νN rappresentano gli stati energetici delle singole molecole. Nota bene

che se Vnon−bonded(rN) = 0 l’equazione 3.11 e’ esatta (non costituisce cioe un approssi-

mazione). La funzione di partizione di un gas di molecole non interagenti (attraverso

interazioni di tipo non-legato) e’ quindi data da

Q =∑ν

e−βEν =1

N !

∑ν1,ν2,ν3...νN

e−β(Eν1+Eν2+Eν3+...EνN ) (3.12)

Il fattore N ! deriva come abbiamo visto per la Eq. 3.1 dall’indistinguibilita’ delle

molecole. Dalla 3.12 otteniamo che:

Q =1

N !

∑ν1,ν2,ν3...νN

(e−βEν1

)·(e−βEν2

)· ...·

(e−βEνN

)Q =

1

N !

(∑ν1

e−βEν1

)·

(∑ν2

e−βEν2

)· ...·

(∑νN

e−βEνN

)

Q =1

N !

(∑µ

e−βEµ

)N

(3.13)

Nella Eq. 3.13 si e’ sfruttato il fatto che le somme sono tutte uguali essendo estese sugli

stati energetici molecolari relativi a molecole identiche e di utilizzato l’indice µ (in luogo

di ν ) per rimarcare il fatto che questo indice va sui soli stati energetici molecolari. Se

definiamo

qMOL =∑µ

e−βEµ (3.14)

come funzione di partizione molecolare, utilizzando la Eq. 3.13 e la definizione 3.14

abbiamo che la funzione di partizione del sistema costituito dal gas perfetto di molecole

identiche puo essere scritta come

Q =qNMOL

N !(3.15)

Per poter proseguire nella nostra derivazione e calcolare espressamente la funzione di

partizione del gas di molecole, assumeremo che l’energia molecolare sia data dalla somma

dell’energia elettronica, vibrazionale, rotazionale e traslazionale:

Eµ = E(el)µ1

+ E(vib)µ2

+ E(rot)µ3

+ E(trasl)µ4

(3.16)

52

Nella equazione sopra, gli indici µ1, µ2, µ3, µ4 corrono rispettivamente sugli stati

elettronici, vibrazionali, rotazionali e traslazionali della molecola.

La Eq. 3.16, al contrario delle Eqs. 3.11, 3.12, 3.13 che sono esatte nell’ipote-

si che le molecole non interagiscano, costituisce una approssimazione doppia. Nella

separazione tra energia dei nuclei (e cioe vibrazionale, rotazionale, traslazionale) ed

elettronica, abbiamo utilizzato l’approssimazione di Born-Oppenheimer, mentre nella

separazione dei gradi di liberta nucleari della molecola e’ possibile far vedere che solo

l’energia traslazionale si puo effettivamente separare rigorosamente dalle restanti en-

ergie vibrazionali e rotazionali. Queste ultime sono in realta accoppiate (per esempio

dal termine di Coriolis). La Eq. 3.16 (sia per quanto riguarda l’approssimazione di Born-

Oppenheimer, sia per quanto riguarda il termine di accoppiamento vibro-rotazionale) e’

tanto piu vera quanto piu’ la temperatura e’ bassa. Per le temperature ordinarie (300

K), la Eq. 3.16 costituisce un ottima approssimazione per una larga parte di molecole

di interesse biologico. La Eq. 3.16 implica che la funzione di partizione molecolare, Eq.

3.14 possa essere fattorizzata nel seguente modo

qMOL =∑µ

eβEµ =∑

µ1µ2µ3µ4

e−β(E(el)µ1

+Eµ2(vib)+Eµ3(rot)+Eµ4)(trasl)

=∑µ1

e−βE(el)µ1

∑µ2

e−βEµ2(vib)∑µ3

e−βEµ3(rot)∑µ4

e−βEµ4(trasl)

= qel · qvib · qrot · qtrasl (3.17)

in cui abbiamo definito

qel =∑µ1

e−βE(el)µ1 funzione di partizione elettronica (3.18)

qvib =∑µ2

e−βE(vib)µ2 funzione di partizione vibrazionale (3.19)

qrot =∑µ3

e−βE(rot)µ3 funzione di partizione rotazionale (3.20)

qtrasl =∑µ4

e−βE(trasl)µ4 funzione di partizione traslazionale (3.21)

53

Queste quattro funzioni di partizione possono essere calcolate indipendentemente.

Vediamo prima il calcolo della funzione di partizione elettronica. In un insieme di

molecole non interagenti allo stato gassoso, in ogni istante praticamente tutte le molecole

si trovano nello stato elettronico fondamentale. Abbiamo infatti visto che e−βEν/QNVT

rappresenta la probabilita canonica di trovare il sistema nello stato ν alla temperatura

β. Dati due stati elettronici molecolari E0 (fondamentale) ed E1 primo stato eccitato,

il rapporto tra le probabilita e dato da:

ρ0

ρ1

=e−βE0/qel

e−βE1/qel= e−β(E0−E1)

La differenza (E0−E1) per gli stati elettronici vale circa 200 nm (nanometri), cioe circa

4 eV (elettronvolt), pari a 6.5· 10−22 kJ ovvero a 6.5· 10−22· 6.023· 1023 kJ/mol ' 4.0· 102

kJ /mol dunque:

ρ0

ρ1

= e400RT = e

4002.4 ∼ e150 = 1050 (3.22)

dove a 300 K abbiamo che RT = kBT · 6.023· 1023 ∼ 2, 4 Kj/mole. Sulla base della Eq.

3.22, troviamo che lo stato fondamentale e circa 1050 volte piu probabile alla temperatura

ordinaria di essere osservato rispetto allo stato eccitato (al primo!). Da cio consegue che

l’unico termine che conta nella sommatoria 3.18 e’ il primo, quello cioe corrispondente

all’energia dello stato elettronico fondamentale E0. Dunque possiamo tranquillamente

scrivere che

qel = e−βE0 (3.23)

Proseguiamo il calcolo della funzione di partizione molecolare qMOL, Eq. 3.17, con il

calcolo della qvib, definita nella Eq. 3.19, ovvero della funzione di partizione vibrazionale.

Se M e’ il numero di atomi che compongono la molecola, il numero di modi normali di

vibrazione di una molecola e’ pari a 3M-6. Per esempio la molecola di acqua ha 3×3−3 =

54

3 modi normali di vibrazione, due di tipo stretching ed un bending. Nell’approssimazione

armonica, l’energia vibrazionale e data dalla somma delle energie vibrazionali dei modi

normali di vibrazione. Dunque abbiamo in generale che

qvib =∑

ν1,ν2,...,ν(3M−6),

e−β

(Eν1+Eν2+...+Eν(3M−6)

)(3.24)

dove la sommatoria e’ estesa a tutti gli stati energetici di tutti i modi normali di

vibrazione. Il modo normale i-esimo ha energia quantistica pari ha

Eνi = hνi

(n+

1

2

)con n = 0, 1, 2... (3.25)

dove νi rappresenta la i-esima frequenza normale di vibrazione. Sostituendo la 3.25

nella Eq. 3.24, sfruttando le proprieta degli esponenziali e fattorizzando le sommatorie,

otteniamo che

qvib =∑n

e−βhν1(n+ 12)·∑n

e−βhν2(n+ 12)· ...·

∑n

e−βhν(3M−6)(n+ 12)

qvib = q1vib· q2

vib· ...· q3M−6vib (3.26)

In cui si e’ definito la funzione di partizione dell’-i-esimo modo normale come

qivib =∑n

e−βhνi(n+ 12) (3.27)

N.B: le sommatorie nella Eq. 3.26 non sono tutte uguali poiche le frequenze vibrazionali

νi sono in generale diverse tra loro. Il calcolo della funzione di Eq. 3.27 puo essere

effettuato sfruttando la serie geometrica. Poniamo

e−βhνi = x (3.28)

Poiche sia le frequenze di vibrazione che β sono numeri positivi, l’esponente nella Eq.

3.28 sara un numero minore di 1 (e maggiore di zero) e dunque 0 < x < 1. Possiamo

scrivere allora ∑n

e−βhνi(n+ 12) = e−

βhνi2

∞∑n=0

xn = e−βhνi

2 · 1

1− x

= e−βhνi

2 · 1

1− e−βhνi(3.29)

55

Dunque la funzione vibrazionale sara data da:

qvib =3M−6∏i=1

e−βhνi

2 · 1

1− e−hνikBT

(3.30)

Se definiamo la quantita En0

En0 =1

2

3M−6∑i

hνi (3.31)

come energia vibrazionale di punto zero, la Eq. 3.30 puo essere riscritta come

qvib = e−βEn0

3M−6∏i=1

1

1− e−hνikBT

(3.32)

Il termine En0 si incorpora spesso nell’energia elettronica, definendo E ′0 = E0 + En0.

Abbiamo pertanto che

qelqvib = e−βE′0

3M−6∏i=1

1

1− e−hνikBT

(3.33)

L’energia di punto zero e’ una energia di natura quantistica. A temperatura zero l’ener-

gia di un oscillatore classico vale infatti zero, mentre quella di un oscillatore quantistico

vale hν/2. Ad alta temperatura, cioe quando kBT hν, oscillatori quantistici e clas-

sici si comportano allo stesso modo. Alle alte temperature, l’energia di un oscillatore

classico e’ data infatti da kBT , indipendentemente dalla sua frequenza di vibrazione.

L’energia di un oscillatore quantistico di frequenza hν puo essere ricavata dalla funzione

di partizione vibrazionale utilizzando la Eq. 2.24:

EQM =1

2hν − ∂ ln(1− e−

hνkBT )−1

∂β

=1

2hν +

hv

(1− e−hνkBT )

(3.34)

Nel grafico di Figura 3.1 sono riportati gli andamenti dell’energia classica e dell’energia

quantistica di un oscillatore armonico. Ad alta temperatura, abbiamo che (1−e−hνkBT ) '

hν/kBT e dunque EQM ' 12hν + kBT . Come si vede, gli effetti quantistici relativi al-

l’energia vibrazionale di punto zero, sono importanti soltanto per basse temperature

56

e cominciano a svanire gia quando kBT/hν > 2. In particolare, i modi normali di vi-

brazione con frequenze al di sotto dei 200 cm−1 si comportano in maniera essenzialmente

classica a temperatura fisiologica. Veniamo adesso al calcolo di qrot (Eq. 3.20 ) nella

Figura 3.1: Andamento quantistico e classico dell’energia di un oscillatore armonico in

funzione della temperatura

funzione di partizione molecolare Eq. 3.17, ovvero calcoliamo la funzione di partizione

rotazionale. Nel caso rotazionale gli effetti quantistici sono ancora meno importanti

che nel caso vibrazionale. Le “frequenze” rotazionali sono infatti centinaia di volte piu

basse di quelle vibrazionali e gli effetti quantistici si possono osservare solo alle bassis-

sime temperature e/o per molecole a bassissimo momento di inerzia (come H2). E’

possibile mostrare3 che la funzione di partizione (classica) rotazionale di una molecola

diatomica e’ data da

qrot =8π2IkBT

h2(per una molecola diatomica)

3La dimostrazione parte dalla definizione di energia cinetica rotazionale come erot = 12Iω

2 e prosegue

con l’integrazione della funzione di partizione sulle coordinate rotazionali del sistema.

57

nella quale I = mr2 e il momento di inerzia attorno all’asse perpendicolare al legame.

Nel caso di una molecola poliatomica noi avremo che

qrot = π1/2

(8π2I1kBT

h2

)1/2(8π2I2kBT

h2

)1/2(8π2I3kBT

h2

)1/2

(3.35)

dove I1, I2 e I3 sono i momenti di inerzia attorno agli assi principali. E veniamo infine

al calcolo della funzione di partizione traslazionale. Questo calcolo lo abbiamo gia fatto

per il gas ideale monoatomico, quando abbiamo derivato la funzione di partizione, Eq.

3.10. Sulla base di quella equazione possiamo facilmente derivare che la funzione di

partizione molecolare traslazionale e’ data da

qtrasl =

√(2πmkbT

h2

)3

V (3.36)

dove m e chiaramente la massa della molecola e V il volume del sistema canonico che

contiene il gas di molecole non interagenti.

Raccogliendo tutti i termini derivati per le funzioni di partizione elettronica, vi-

brazionale, rotazionale e traslazionale (vedi Eqs. 3.23, 3.30, 3.32, 3.35, 3.36) abbiamo

l’espressione finale generale per la funzione di partizione molecolare definita nella Eq.

3.17:

qMOL = e−βE′0·

(3M−6∏i=1

1

1− e−hνikBT

)·√π(8π2kbT )3I1I2I3

h3·

√(2πmkbT

h2

)3

V (3.37)

In realta la Eq. 3.38 dipende da volume V del sistema ed e’ dunque una quantita

estensiva, mentre il nome “funzione di partizione molecolare” si riferisce ad una quantita

intensiva (cioe per molecola). La funzione di partizione molecolare propriamente detta

si ottiene dalla Eq. 3.37 dividendo per il volume V del sistema:4

qMOL = e−βE′0·

(3M−6∏i=1

1

1− e−hνikBT

)·√π(8π2kbT )3I1I2I3

h3·

√(2πmkbT

h2

)3

(3.38)

4E’ importante notare (per il futuro, quando si parlera di stati standard) che la funzione di partizione

molecolare definita nella 3.38 ha le dimensioni di un volume inverso attraverso il fattore traslazionale.

58

La funzione di partizione molecolare qMOL puo essere calcolata se si conosce la strut-

tura della molecola, di qualunque molecola; la Eq. 3.38 puo essere calcolata, per esempio,

per la molecola “farmaco” e la molecola “proteina”. Vediamo ora perche la qMOL e cosı

importante in termodinamica statistica chimica e come essa sia connessa all’attivita’

(efficacia) di un farmaco.

3.3 Equilibrio chimico per una miscela ideale di gas

A tale scopo consideriamo una reazione chimica del tipo:

A + B AB (3.39)

A+B − AB = 0

Nell’equilibri sopra la specie A potrebbe riferirsi ad un farmaco e la specie B ad una

proteina. La stechiometria impone il seguente bilancio di massa:

δnA = δnB = −δnAB (3.40)

in cui δnA, δnB, δnAB, rappresentano le variazioni del numero di moli delle specie A, B

e AB. La variazione di energia libera per effetto di variazioni della composizione chimica

del sistema (cioe del numero di moli) e’ data da

dA = µAδnA + µBδnB + µABδnAB (3.41)

all’equilibrio, quando cioe l’energia libera di Hellmotz e’ minima (dA = 0 ) avremo che

µAδnA + µBδnB + µABδnAB = 0 (3.42)

Utilizzando la Eq. 3.40 nella Eq. 3.42, otteniamo che

µAδnA + µBδnA − µABδnA = 0

59

e cioe:

µA + µB = µAB (3.43)

La Eq. 3.43 stabilisce che i potenziali chimici a sinistra nell’equilibrio chimico 3.39

devono (sommati) essere uguali ai potenziali chimici (sommati) a destra. Ora, se il

sistema nel quale avviene la reazione 3.39 costituisce un insieme canonico, sappiamo

che in questo caso l’energia libera di Hellmotz e’ legata alla funzione di partizione del

sistema dalla relazione A = −kBT lnQ (vedi Eq. 2.25). Nel caso in cui il nostro sistema

sia costituito da una miscela gassosa ideale di A, B, AB all’equilibrio, abbiamo che la

funzione partizione puo essere scritta come5

Q =1

NA!(qA)NA· 1

NB!(qB)NB · 1

NAB!(qAB)NAB (3.44)

Dove:

qA = funzione di partizione molecolare della specie A

qB = funzione di partizione molecolare della specie B

qAB = funzione di partizione molecolare della specie AB

Utilizzando la relazione Eq. 2.25 tra energia libera e funzione di partizione canonica Q,

e sfruttando le proprieta dei logaritmi abbiamo che

A = −kBT lnQ

= −kBT ln

(1

NA

· qNAA1

NB

· qNBB1

NAB

· qNABAB

)= −kBT

[ln

(1

NA

· qNAA)

+ ln

(1

NB

· qNBB)

+ ln

(1

NAB

· qNABAB

)]= AA + AB + AAB (3.45)

5La Eq. 3.44 si riferisce ad un sistema di molecola A, B ed AB ideale. Dunque le molecole non

interagiscono tra loro. Ora dentro il recipiente noi abbiamo supposto che possa avvenire una reazione

chimica del tipo 3.39 e questa reazione avviene a seguito dell’interazione di A con B contraddicendo

l’assunzione di idealita. In realta per miscela “ideale” noi intendiamo una miscela di A, B e AB in cui

il potenziale di interazione non bonded tra tutte le particelle presenti sia nullo e la reazione avvenga

solo per urti causali tra le molecole che non si attraggono ne’ si respingono.

60

in cui AA, AB, AAB sono le energia libera della specie A, B e AB nel nostro insieme

canonico posto alla temperatura T. La quantita AA, sulla base dell’ultimo passaggio

nella Eq. 3.45, e’ definita come

AA = −kBT ln

(1

NA!· qNAA

)Il potenziale chimico (per molecola) della specie A pio’ essere calcolato derivando l’e-

spressione sopra rispetto al numero di di molecole della specie A ed utilizzando ancora

una volta le proprieta dei logaritmi:

µA =∂AA∂NA

= −kBT∂

∂NA

ln

(1

NA!· qNAA

)=

∂

∂NA

[kBT ln (NA!)− kBTNA ln qA] (3.46)

Usiamo l’approssimazione di Stirling6 per il logaritmo di un numero fattoriale:

lnN ! ∼= N lnN −N (3.47)

Con questa approssimazione Eq. 3.47, la Eq. 3.46 puo essere riscritta come

µA =∂

∂NA

[kBT (NA lnNA −NA)− kBTNA ln qA] ,

da cui:

µA = kBT (lnNA − 1 + 1)− kbT ln qA

= kBT lnNA − kBT ln qA (3.48)

6Questa approssimazione e’ pessima per piccoli valori di N, ma diventa sempre migliore quando N

diventa grande. Per esempio per N=5, ln 5! vale 3.4 mentre 5 ln 5− 5 vale 3.04. Per N=10, l’approssi-

mazione di Stirling e’ gia migliore. Abbiamo infatti ln 10! = 13.7 e 10 ln 10 − 10 = 13.0. Nel nostro

caso N e’ grandissimo (dell’ordine del numero di Avogadro, e dunque l’approssimazione di Stirling e’

praticamente esatta

61

Analogamente per le energia libere AB ed AAB relative alle specie B ed AB, otterremo

che

µB = kBT lnNB − kBT ln qB (3.49)

µAB = kBT lnNAB − kBT ln qAB (3.50)

Nelle Eqs. 3.48, 3.50, 3.50 siamo riusciti a scrivere il potenziale chimico (per moleco-

la) relativo alle specie A, B e AB in funzione della funzione di partizione molecolare

corrispondente qA, qB, qAB. Ma dalla relazione generale per la funzione partizione

molecolare (Eq. 3.37) possiamo ricavare che

qMOL =

√(2πmkBT

h2

)3

V · qint (3.51)

in cui abbiamo definito la funzione di partizione molecolare per i gradi interni qint

qint = qrot· qvib· qel (3.52)

Le equazioni 3.51, 3.52 si applicano a ciascuna delle specie presenti nella nostra miscela

ideale in equilibrio chimico. Per esempio per la specie A, possiamo utilizzare la funzione

di partizione molecolare del tipo Eq, 3.52 relativa alla molecola di tipo A nella espressione

per il potenziale chimico di A:

µa = kBT lnNA − kBT ln

√(2πmAkBT

h2

)3

V · qint

µAkBT

= lnNA − lnV − ln

√(2πmAkBT

h2

)3

· qint

µAkBT

= ln ρA + ln

[(h2

2πmAkBT

)1/2]3

· 1

qint

βµA = ln

(ρA·λ3

A

qint

)(3.53)

62

in cui β = 1/kBT , abbiamo sfruttato il fatto che la densita per la specie A e’ definita

dalla ρA = NA/V ed abbiamo infine posto che

λA =

√h2

2πmAkBT(3.54)

La quantita λA e detta anche lunghezza d’onda termica (della specie A nel caso specifico)

perche ha proprio le dimensioni di una lunghezza d’onda.7

In modo del tutto analogo ricaviamo che

βµB = ln

(ρB·λ3

B

qint

)(3.55)

βµAB = ln

(ρAB·λ3

AB

qint

)(3.56)

Moltiplicando la Eq. 3.43 a destra e a sinistra per β e sfruttando le relazioni 3.53, 3.55,

3.56 noi abbiamo che:

βµA + βµB = βµAB

ln

(ρA·λ3

A

qAint

)+ ln

(ρB·λ3

B

qBint

)= ln

(ρAB·λ3

AB

qABint

)(3.57)

Sfruttando le proprieta dei logaritmi e considerando che nella miscela ρA = NA/V ≡ [A]

(e cioe la densita della specie A e’ equivalente alla sua concentrazione espressa in molecole

per litro, per esempio), la Eq. 3.57 puo essere riarrangiata nel modo seguente

ln

(ρAρBρAB

)= ln

(qAintq

Bint

qABint

λ3AB

λ3Aλ

3B

)= ln

(qAqBqAB

)qA· qBqAB

=[A][B]

[AB]= Keq (3.58)

in cui si e’ tenuto conto, sulla base della Eq. 3.38, che qA = qAint/λ3A. La relazione

Keq =[A][B]

[AB]=qA· qBqAB

(3.59)

7Piu piccola e’ la lunghezza termica e meno “quantistica” e’ la particella (molecola) in questione.

Nel caso della nostra miscela, quando la lunghezza d’onda termica e’ molto minore della distanza media

tra le particelle (che e’ circa uguale a (V/N)1/3) allora il gas ideale segue la statistica classica.

63

e’ uno dei piu grandi successi della meccanica statistica e consente, in linea di principio,

di calcolare le costanti di equilibrio delle reazioni chimiche se si conoscono le funzioni

di partizione molecolari di reagenti e prodotti senza conoscere alcun dettaglio circa il

meccanismo molecolare attraverso il quale la reazione chimica avviene.8

Riscriviamo adesso l’equazione 3.57 nella seguente forma

µA + µB − µAB = kBT ln

(ρAρBρAB

)− kBT ln

(qAqBqAB

)(3.61)

Dalla termodinamica classica, noi sappiamo che il potenziale chimico di un gas ideale e’

dato da

µ = µ0 + kBT ln ρ (3.62)

nella quale µ0 e’ il potenziale chimico del gas nello stato standard (concentrazione 1M

a T=25 e P=1 atm). Se utilizziamo l’espressione Eq. 3.62 per la nostra miscela ideale

sostituendola nella Eq. 3.61, otteniamo che

µ0A + µ0

B − µ0AB = −kBT ln

(qAqBqAB

)(3.63)

Definendo la quantita

∆Areaz0 = µ0

A + µ0B − µ0

AB (3.64)

come la differenza di energia libera standard per la reazione A+B = AB , la Eq. 3.63

puo essere riscritta nella forma

e−β∆Areaz0 =

qA· qBqAB

= Keq (3.65)

8Seguendo una procedura del tutto analoga a quella adottata per ricavare la Eq. 3.59, e’ possibile

estendere l’equazione 3.59 alla reazione generica aAr + bBr + cBr + ... uAp + vBp +wBp + ..., nella

quale Ar, Br.. sono i reagenti e Ap, Bp.. i prodotti. In questo caso si perviene all’espressione generale

qaApqbBp

qcCp...

quArqvBr

qwCr...

=[Ap]a[Bp]b[Cp]c..

[Ar]u[Br]v[Cr]w..(3.60)

64

La equazione sopra mette in relazione l’energia libera di reazione con il rapporto tra

le funzioni di partizione molecolare. Tale rapporto e’ a sua volte uguale, secondo la

relazione Eq. 3.59, alla costante dell’equilibrio chimico.

Concludiamo con una ultima importante considerazione. Nella Eq. 3.59, la costante

di equilibrio per la nostra reazione A+B=AB, che esemplifica l’interazione di una specie

A (per esempio un farmaco) con una specie B (una proteina), ha le dimensioni dell’in-

verso di un volume, cosı come il rapporto delle funzioni di partizione molecolari (volume

indipendenti). Gli argomenti dei logaritmi nella Eq. 3.57 sono in effetti adimensionali,

come dovrebbe essere. Quando la Eq. 3.57 viene riarrangiata nella forma 3.59, gli ar-

gomenti dei logaritmi a destra e a sinistra hanno le dimensioni entrambi di un volume

inverso e le unita di misura - se si vuole che le cose tornino - devono essere scelte in

maniera coerente. Anche nella Eq. 3.63 bisogna stare attenti. Se si vuole che le cose

tornino, bisogna che l’argomento del logaritmo a destra sia adimensionale. La Eq. 3.65

per come e’ stata ricavata dalla Eq. 3.61 deve essere intesa come riferita agli stati stan-

dard che come detto, per le concentrazioni, corrispondono ad 1 M per tutte le specie.

La Eq. 3.65 deve essere dunque intesa come

e−β∆Areaz0 = Keq

=

(ρAρBρAB

)(ρ0AB

ρ0Aρ

0B

)(3.66)

Se per le concentrazioni ρ nella costante di equilibrio si usano le unita di moli/litro allora

la concentrazione dello stato standard di tutte le specie e 1 M, il che implica che

ρ0AB

ρ0Aρ

0B

=1

1× 1= 1 (3.67)

Dunque il fattore relativo alle concentrazioni degli stati standard delle tre specie vale 1.

Ma se invece delle concentrazioni in moli/litro, si usassero le densita ρ = N/V in termini

di volume per molecola per definire la costante di equilibrio (per esempio esprimendo

il volume in A3) le dimensioni della Keq sarebbero quelle di un volume. In quel caso il

65

fattore relativo alle concentrazioni degli stati standard sarebbe pari a

ρ0AB

ρ0Aρ

0B

=6.0232023/1027

6.0232023/1027 × 6.0232023/1027=

1

6.0232023/1027=

1

1661(3.68)

ovvero bisognerebbe dividere la costante di equilibrio (espressa in A3) per 1661 A3.

Vediamo ora una applicazione semplice del calcolo della costante di equilibrio per

una reazione chimica in fase gassosa, rimandando al capitolo successivo la discussione

dell’applicazione della Eq. 3.59 al problema che ci sta piu a cuore, ovvero il calcolo

dell’affinita’ di un farmaco A per una proteina B. Consideriamo il seguente equilibrio

relativo alla formazione di monossido di azoto a partire da azoto ed ossigeno.

N2 + O2 2NO (3.69)

Sulla base delle Eqs. 3.57 e 3.60, possiamo scrivere che

keq =qN2intq

O2int

qNOint qNOint

λ3NOλ

3NO

λ3N2λ3O2

(3.70)

La funzione di partizione interna qint relativa alle diverse specie e’ definita nella Eq. 3.52

in cui i tre termini rotazionale,9 vibrazionale ed elettronico sono a loro volta definiti nelle

Eqs. 3.35, 3.33. Per la specie N2, per esempio, essa vale

qN2int = e−βE

′N2

(1

1− e−hνN2kBT

)2IN2kBT

h2(3.71)

in cui IN2 e’ il momento di inerzia della molecola attorno all’asse perpendicolare al

legame, νN2 e’ la frequenza di vibrazione di stretching ed E ′N2e’ l’energia di stato fonda-

mentale comprensiva dell’energia di punto zero hνN2/2. Assumiamo ora che la reazione

avvenga ad elevata temperatura, quando cioe kBT hν per tutte le specie presenti.

In questo caso (vedi Figura 3.1 e discussione relativa) abbiamo visto che la funzione di

partizione vibrazionale e’ data da

1

1− e−hνkBT

' kBT

hν(3.72)

9N.B.: la funzione di partizione rotazionale per le molecole N2, O2 ed NO e’ ovviamente quella per

molecole diatomiche.

66

valida per tutte le specie N2, O2, NO. Utilizzando questa approssimazione nella Eq.

3.71, ricavando in maniera del tutto analoga le qint per le altre specie e sostituendo il

tutto nella Eq. 3.70, si ottiene la seguente espressione per la costante di equilibrio:

Keq = e−β∆E

(IN2IO2

INOINO

)(m

3/2N2m

3/2O2

m3/2NOm

3/2NO

)(3.73)

dove mN2 , mO2 , mNO sono le masse molecolari delle tre specie coinvolte e

∆E = 2E ′NO − E ′N2− E ′O2

(3.74)

e’ definita come energia (o entalpia se siamo a P costante) della reazione. In questo caso

non ci sono problemi di stati standard o dimensioni perche la costante di equilibrio e’

adimensionale per stechiometria. I reagenti ed i prodotti della reazione chimica hanno

piu o meno momenti di inerzia e masse molecolari simili in modo tale che i fattori

dipendenti da momenti e masse valgono all’incirca 1 e la costante di equilibrio vale

all’incirca

Keq ' e−β∆E (3.75)

Dunque la costante di equilibrio per la reazione Eq. 3.69 dipende in ottima approssi-

mazione solo dalla differenza tra le energie di stato fondamentale delle specie (inclusa

l’energia vibrazionale di punto zero). Queste energie, nel caso di molecole piccole come

quelle in questione, possono essere calcolate facilmente mediante programmi ab initio.

Come vedremo nel prossimo capitolo, questa approssimazione, in base alla quale la Keq

dipende in sostanza dalle energie di stato fondamentale delle specie coinvolte, contin-

ua ad essere piuttosto buona anche per alcune classi di interazioni farmaco-proteina.

Nel caso di molecole complesse come farmaci e proteine, non si possono tuttavia usare

metodi ab initio per valutare la differenza

∆E = EAB − EA − EB. (3.76)

che risulterebbero troppo costosi dal punto di vista computazionale. Piuttosto, bisogna

ricorrere ad una parametrizzazione della superficie di Born-Oppenheimer relativa allo

67

stato elettronico fondamentale delle specie. Questa parametrizzazione e’ costituita, come

illustrato nel capitolo 1, dal cosiddetto force field (Vedi Eqs. 1.2 e 4.62)

68

Capitolo 4

Costante di Affinita

Farmaco-Proteina

Nel capitolo precedente (vedi Eq. 3.63) abbiamo visto che per una miscela ideale di gas

con specie A, B e AB, in equilibrio chimico secondo la Eq. 3.39, l’energia libera standard

di dissociazione(per mole)1 per la reazione AB A+B vale :

∆GAB = −RT ln

(qAqBqAB· 1661A3

)(4.1)

Nella Eq. 4.1 si e’ assunto le funzioni di partizione molecolari volume dipendenti siano

misurate in A3 in modo tale che l’argomento del logaritmo sia adimensionale. Nota

l’energia libera standard (riferita, lo ricordiamo, a reagenti e prodotti a T=25 e P=1

Atm in concentrazioni 1 M), la costante di dissociazione (in M/litro) per la reazione

AB A+B puo essere calcolata in accordo alla

Kd = e−∆GAB/RT ≡ [A][B]

[AB](4.2)

1N.B. Per passare da energia libera per molecola [ potenziale chimico] come nella Eq. 3.63 ad energia

libera per mole basta moltiplicare a destra e a sinistra per NA=numero di Avogadro in modo tale che

R = Costante dei gas perfetti = NAkB .

69

Le due equazioni 4.2 e 4.1 valgono come detto per una miscela di gas ideali, ma possono

essere utilizzate anche per soluzioni diluite delle specie A B ed AB. In cio che segue

assumeremo che A sia la molecola farmaco, B sia la molecola proteina e AB sia il

complesso farmaco-proteina. Assumeremo pertanto che la reazione di formazione del

complesso avvenga in un solvente (generalmente acqua) con A, B e AB in concentrazioni

tali che si possano trascurare le interazioni non legate (ancorche mediate dal solvente)

tra le varie specie. In queste condizioni, con alcuni caveat che verranno discussi nel

capitolo seguente, Le Eqs. 4.1, 4.2 continuano a valere.

Nell’ambito della chimica farmaceutica, la costante di dissociazione farmaco (A)

proteina (B) si misura normalmente in nanomoli o micromoli per litro (simbolo nM o

µM). Farmaci potenti hanno costanti di dissociazione nell’ordine delle nanomoli/litro.

Cio significa che il denominatore nella 4.2 e molto piu grande (all’equilibrio quando

∆µ = 0, vedi Eq. 3.42) del numeratore. Corrispondentemente l’energia di dissociazione

standard ∆G = G0A +G0

B −G0AB e positiva e grande per farmaci potenti.2

Come e’ noto dalla termodinamica classica, a temperatura e pressione costanti

abbiamo che (vedi anche Eq. 2.21):

∆G0AB = ∆H0

AB − T∆S0AB (4.3)

Il ∆G0AB ( che e legato alla costante di dissociazione secondo la Eq. 4.2) e dato dalla

somma di due termini di segno generalmente opposto: l’entalpia standard ∆H0AB di

dissociazione (positiva) e l’entropia moltiplicata per T (negativa). Cio implica che tanto

l’energia quanto l’entropia relative allo stato associato AB siano minori dell’energia e

dell’entropia dello stato dissociato A + B. Dunque, normalmente, la formazione del

2In generale con energia libera di stato standard ci si riferisce all’energia libera di Gibbs, che coinvolge

la variazione di entalpia ∆H = E + p∆V (a pressione p = 1 Atm costante) in luogo della variazione

di energia ∆E. Tuttavia in fase liquida (dove avviene l’interazione farmaco-proteina) il termine p∆V

vale qualche frazione di caloria per mole ed e’ solitamente trascurabile rispetto all’energia libera di

dissociazione (dell’ordine delle kcal/mole).

70

complesso farmaco-proteina e favorita dal fattore energetico (positivo) e sfavorita dal

fattore entropico (negativo). Nell’esempio di figura, per la reazione di formazione del

Figura 4.1: Bilancia entropia/entalpia per l’affinita di un farmaco

complesso AB a partire dai prodotti A, B, l’affinita di un farmaco e’ tanto maggiore

quanto minore e’ l’energia libera di formazione, cioe quanto piu bassa e’ l’energia libera

del complesso in relazione a quello dello stato dissociato. Tale affinita si puo aumentare

o aumentando il “peso” dell’entalpia a destra nella bilancia (in modo tale che l’energia

libera complessivamente scenda) o diminuendo il peso dell’entropia di formazione a

sinistra nella bilancia. Non potendo agire sulla proteina, si diminuisce la perdita di

entropia di formazione del complesso, per esempio, rendendo i farmaci piu “rigidi”. Ne

consegue che, dal punto di vista energetico, il farmaco migliore e’ quello che si associa con

piu forza alla proteina (normalmente attraverso interazioni di tipo non legato), mentre,

dal punto di vista entropico, il farmaco migliore e quello con maggiore rigidita nello stato

dissociato.3 Vedremo adesso come, attraverso la relazione 4.1 e la fattorizzazione delle

3E’ questa una delle ragioni per le quali gran parte dei farmaci efficaci sono costituiti da gruppi

71

funzioni di partizione molecolari, l’energia libera di dissociazione 4.1 venga espressa come

somma di contributi che possono, in molti casi, essere calcolati esattamente a partire

dalla conoscenza della struttura molecolare del farmaco e della proteina. Lasciando per

adesso perdere il fattore 1661 A3 relativo alla “adimensionalizzazione” del rapporto tra

le funzioni partizione quando queste sono espresse in A3, sfruttiamo la relazione 3.17

per ciascuna delle funzioni di partizione molecolari delle specie A,B e AB scrivendo

qAqBqAB

=qelAq

elB

qelAB

· qtraslA qtrasl

B

qtraslAB

· qrotA qrot

B

qrotAB

· qvibA qvib

B

qvibAB

(4.4)

Inserendo la Eq. 4.4 nella equazione 4.1 otteniamo

∆G0AB = −RT ln

(qelAq

elB

qelAB

)−RT ln

(qtraslA qtrasl

B

qtraslAB

)−

− RT ln

(qrotA qrot

B

qrotAB

)−RT ln

(qvibA qvib

B

qvibAB

)∆G0

AB = ∆GelAB +∆Gtrasl

AB +∆GrotAB +∆Gvib

AB (4.5)

Di nuovo, tutto cio vale per specie non interagenti o soluzioni diluite.

Nella Eq. 4.5 abbiamo definito

∆GelAB = −RT ln

qelAq

elB

qelAB

= GelA +Gel

B −GelAB (4.6)

∆GtraslAB = −RT ln

qtraslA qtrasl

B

qtraslAB

= GtraslA +Gtrasl

B −GtraslAB (4.7)

∆GrotAB = −RT ln

qrotA qrot

B

qrotAB

= GrotA +Grot

B −GrotAB (4.8)

∆GvibAB = −RT ln

qvibA qvib

B

qvibAB

= GvibA +Gvib

B −GvibAB (4.9)

aromatici variamente sostituiti e/o variamente connessi tra di loro. Un elevato tasso di aromaticita

diminuisce infatti il disordine torsionale nelle molecole riducendo le perdite potenziali di entropia per

effetto della formazione del complesso.

72

Le funzioni di partizione molecolari Eqs. 3.35, 3.27 e 3.23 sono adimensionali e solo la

funzione di partizione traslazionale, Eq. 3.36, ha le dimensioni di un volume. Dunque

lo stato standard (il fattore 1660A3, volume per molecola) entra solo nel calcolo di

∆GtraslAB . Tutte le altre dipendono solo dalla temperatura. Sulla base della fattorizzazione

dell’equazione 4.5 e delle equazioni 4.6, 4.7, 4.8 e 4.7 che ne sono derivate, e possibile

scrivere:

∆GelAB +∆Gtrasl

AB +∆GrotAB +∆Gvib

AB = −RT lnKeq (4.10)

Tenendo conto della Eq. 4.3, possiamo affermare che la Keq che noi misuriamo diretta-

mente negli esperimenti in vivo ed in vitro dipende da otto addendi “molecolari”ovvero

4 energie di dissociazione (elettronica, traslazionale, rotazionale e vibrazionale) e 4

corrispondenti entropie:

∆G0AB = ∆Hel

AB +∆HtraslAB +∆Hrot

AB +∆HvibAB −

− T∆SelAB − T∆Strasl

AB − T∆SrotAB − T∆Svib

AB

= −RT lnKeq (4.11)

Come vedremo, le espressioni analitiche di ciascuno di questi 8 termini sono in linea di

principio note se e nota la struttura molecolare e dunque possiamo determinare, dalla

conoscenza della struttura molecolare di farmaco e proteina la costante di equilibrio. Se

si puo fare questo, ne consegue che e possibile, solo attraverso il calcolo, concepire delle

variazioni della molecola farmaco in modo da rendere il ∆G0AB complessivo il piu alto

possibile, agendo su uno o piu degli 8 termini “molecolari ”presenti nella 4.11. Questo

processo si chiama “progettazione razionale di farmaci ”.

73

4.1 Teoria RRHA per il calcolo dell’energia libera

di dissociazione

Il calcolo dell’energia libera di dissociazione basato sul calcolo analitico del rapporto

fra le funzioni di partizione molecolari di farmaco, proteina e complesso prende il nome

di teoria “rigid rotor harmonic approximation” (RRHA). In tale approccio si assume

infatti che il farmaco e la proteina si comportino come rotatori essenzialmente rigidi

soggetti a vibrazioni intramolecolari di tipo armonico. Vediamo ora nel dettaglio il

contributo di ciascuno degli 8 termini nella 4.11. Cominciamo dal termine elettronico

da ∆GelAB = ∆Hel

AB − T∆SelAB :

∆GelAB = −RT ln

qelAq

elB

qelAB

= −RT lne−E

0A/RT e−E

0B/RT

e−E0AB/RT

(4.12)

Dove E0A, E0

B e E0AB sono le energie (per mole) elettroniche delle specie A, B e del

complesso AB. Nella derivazione della Eq. 4.12 abbiamo utilizzato la relazione 3.23

nella quale le energie per molecola sono state sostituite da energia per moli dunque

utilizzando il fattore 1/RT in luogo4 di β. Come accennato nel capitolo precedente (vedi

discussione sulla reazione N2 + O2 = 2NO), queste energie non si calcolano ab initio,

ma utilizzando un force-field (descritto al Capitolo 1), ovvero una parametrizzazione (in

funzione delle coordinate nucleari) dell’energia di stato fondamentale delle specie, A,

B e AB in equilibrio nella miscela. Sfruttando le proprieta degli esponenziali, la 4.12

diviene banalmente:

∆GelAB = −RT (−E0

A/RT − E0B/RT + E0

AB/RT ) = EA + EB − EAB (4.13)

Notare che nel ∆GelAB non c’e il termine entropico perche c’e un solo stato possibile

per le specie dissociate (lo stato fondamentale) e quelle associate, ovvero: SA + SB −4Di fatto si e’ moltiplicato nella Eq. 3.23, nell’argomento dell’esponenziale βE ≡ E/kBT sopra e

sotto per NA.

74

SAB = 0 − 0 = 0. Dal punto di vista elettronico non c’e disordine di nessun tipo.5

Dunque abbiamo che:

∆HelAB = EA + EB − EAB (4.14)

T∆S0AB = 0 (4.15)

Consideriamo ora al termine successivo nell’equazione 4.5, ovvero:

∆GtraslAB = −RT ln

qtraslA qtrasl

B

qtraslAB

·V0 (4.16)

dove V0 e il volume molecolare standard pari a 1661A3. Come abbiamo visto in prece-

denza, l’introduzione del fattore V0 e’ necessario per rendere l’argomento del logaritmo

adimensionale. Il termine individuale dell’energia libera (standard) traslazionale per la

specie A vale:

GtraslA = −RT ln qtrasl

A = −RT ln

√(2πmAkBT

h2

)3

V0

(4.17)

dove mA e’ la massa molecolare della specie farmaco (indicata con il label A).

Ora, noi sappiamo che dG = −TdS − pdV +∑

A,B,AB µdn. Dunque, essendo dG un

differenziale esatto, avremo necessariamente che(∂Gtrasl

A

∂T

)VnB,nAB

= −SA. (4.18)

Sulla base dell’equazione sopra, derivando la 4.17 a sinistra rispetto alla temperatura T ,

5Le cose cambiano un po’ se il farmaco dovesse avere lo stato fondamentale degenere; allora per

effetto dell’associazione questa degenerazione sarebbe rimossa e si riavrebbe una perdita di entropia

per effetto della formazione del complesso.

75

si ottiene facilmente6 l’entropia traslazionale:

−SA =∂

∂T

(−RT ln

[(2πmAkBT

h2

)3/2

V0

)]

= −R ln

[(2πmAkBT

h2

)3/2

V0

]− 3

2R

= −3

2R ln

(2πmAkBT

h2

)−R lnV0 −

3

2R

SA =3

2R ln

(2πmAkBT

h2

)+R lnV0 +

3

2R (4.20)

Osservando infine che

GtraslA = −3

2RT ln

(2πmAkBT

h2

)−RT lnV0 = Etrasl

A − TStraslA (4.21)

otteniamo l’energia standard per mole della specie A:

HtraslA = Gtrasl

A + TStraslA =

3

2RT (4.22)

In modo del tutto analogo, possiamo ricavare le energie standard(considerando nullo

il termine di volume) traslazionali per la specie B proteina e per il complesso farmaco

proteina AB. Nell’espressione 4.22 relativa alla specie A, osserviamo che l’energia tral-

sazionale e’ sempre a 3/2RT indipendentemente dalle caratteristiche molecolari della

specie. Ne consegue che

∆EtraslAB = Etrasl

A + EtraslB − Etrasl

AB =3

2RT +

3

2RT − 3

2RT =

3

2RT, (4.23)

e cioe la differenza di energia traslazionale ai fini della formazione del complesso e sempre

32RT (positiva) indipendentemente dalla struttura molecolare di farmaco e proteina e

6Nella derivata della Eq. 4.17 rispetto a T , occorre tenere in conto che T compare sia nel fattore

RT che moltiplica il logaritmo che nell’argomento del logaritmo stesso e che

RT∂

∂T

ln

[(2πmkbT

h2

)3/2

V0

]= RT

3

2· ∂∂T

lnT =3

2R (4.19)

76

dalla natura dell’interazione. Il risultato 4.22 e una conseguenza del cosiddetto teorema

di equipartizione dell’energia: l’energia (cinetica) si ripartisce per 12RT per ogni grado

di liberta. Il farmaco ha tre gradi di liberta traslazionale esattamente come la proteina

e dunque la loro energia traslazionale (cinetica) sara 32RT in ambedue i casi.

Come abbiamo visto, l’entropia della specie A e’ data dalla

SA =3

2R ln

(2πmakBT

h2

)+R lnV0 +

3

2R (4.24)

e analogamente quella per le specie proteina B e complesso farmaco-proteina AB. Sta-

volta, al contrario di quanto visto nel caso dell’energia traslazionale identica per tutte

le specie, l’entropia traslazionale e’ differente da specie a specie. Tuttavia, le caratteris-

tiche molecolari entrano nella 4.24 soltanto attraverso la massa molecolare mA riferita

nella fattispecie alla specie farmaco A. Possiamo dunque scrivere per le tre specie :

SA =3

2R ln

mA

h2+

3

2R ln(2πkBT ) +R lnV0 +

3

2R

SB =3

2R ln

mB

h2+

3

2R ln(2πkBT ) +R lnV0 +

3

2R

SA+B =3

2R ln

mA +mB

h2+

3

2R ln(2πkBT ) +R lnV0 +

3

2R

(4.25)

Se utilizziamo le Eqs. 4.25 per il calcolo del ∆StraslAB , otterremo

∆StraslAB = Strasl

A + StraslB − Strasl

AB

=3

2R ln

M

h2+

3

2R ln(2πkBT ) +R lnV0 +

3

2R

(4.26)

Nella 4.36, la quantita M ,

M =mAmB

mA +mB

(4.27)

e la massa ridotta del sistema farmaco-proteina.

77

Riassumendo, per i soli contributi traslazionali all’energia libera standard di disso-

ciazione ∆GAB (Eq. 4.2) abbiamo ottenuto che

∆GtraslAB = ∆Etrasl

AB − T∆StraslAB

=3

2RT − 3

2RT ln

M

h2− 3

2RT ln(2πkBT )−RT lnV0 −

3

2RT

= −3

2RT ln

M

h2− 3

2RT ln(2πkBT )−RT lnV0 (4.28)

nella quale abbiamo utilizzato la Eq. 4.36 e la Eq. 4.22. La equazione 4.28 stabilisce

un primo risultato generale importante. Osservando che la massa ridotta del sistema

farmaco-proteina puo essere approssimata con la massa M ' mA del solo farmaco,7 la

Eq. 4.28 puo’ infatti essere scritta come

∆GtraslAB = −RT ln

[(mA

h22πkBT

)3/2

V0

](4.30)

A causa della piccolezza della costante di Planck, h, l’argomento adimensionale com-

plessivo del logaritmo nella Eq. 4.30 e sempre maggiore di 1 e dunque, per quanto

riguarda il solo contributo traslazionale, l’energia libera di dissociazione e sempre neg-

ativa, il che significa che la reazione di dissociazione AB → A + B e sempre termodi-

namicamente favorita. Possiamo facilmente stimare l’ordine di grandezza del contributo

traslazionale all’energia libera. Sostituendo i valori delle costanti h = 6.62×10−34 m2 kg

s−1 k = 1.38×10−23 m2 kg s−2 K−1 ed assumendo una massa tipica per un farmaco pari

a mA = 300 Da = 5× 10−25 kg, otteniamo che il contributo traslazionale all’energia lib-

era di dissociazione vale circa -10 kcal mol−1. Tale significativo contributo, che favorisce

lo stato dissociato e dovuto, in larghissima parte al guadagno entropico che si ottiene

nella dissociazione del complesso che porta a due specie a partire da una, aumentando

7Questo accade perche la massa della proteina mB mA e molto maggiore di quella del farmaco e

dunquemAmB

mA +mB' mAmB

mB' mA (4.29)

78

pertanto il disordine traslazionale nel sistema. Questo contribito deve essere compen-

sato dagli altri contributi perche possano esistere in natura complessi farmaco-proteina

stabili. Vedremo in seguito che i contributi positivi all’energia libera di dissociazione,

cioe quelli che favoriscono la formazione del complesso, possono derivare soltanto dal

termine elettronico gia considerato (vedi Eq. 4.13) e da quello vibrazionale che vedremo

in seguito.

Il discorso per il contributo rotazionale all’energia libera di dissociazione, Eq. 4.1 (e

dunque alla costante di dissociazione farmaco-proteina) e simile a quello sviluppato per

il termine traslazionale. Per una specie generica , abbiamo che (vedi l’equazione 3.35

per la funzione di partizione rotazionale per molecole poliatomiche) e’ data da:

qrot = π1/2

(8π2I1kBT

h2

)1/2(8π2I2kBT

h2

)1/2(8π2I3kBT

h2

)1/2

(4.31)

in cui I1, I2, I3 sono i momenti di inerzia principali della specie in oggetto [N.B. Qui

non ci sono problemi di stato standard: l’argomento del logaritmo e adimensionale].

I momenti di inerzia sono ovviamente proprorzionali alla massa molecolare, farmaco,

proteina o complesso. E’ conveniente, per motivi che saranno piu chiari in seguito,

definire un volume Va, relativo ad esempio alla specie a, tale che

I1I2I3 = m3aV

2a (4.32)

Il volume Va e’ una quantita che si puo’ calcolare esattamente conoscendo la struttura

della molecola, che dipende dalla composizione e forma della specie molecolare a e che

risulta essere circa proporzionale al suo volume. Applicando l’analoga della Eq. 4.18

alla energia libera rotazionale −RT ln qrot ed utilizzando le Eqs. 3.35 e 4.32 otteniamo,

dopo qualche semplice passaggio, l’entropia rotazionale riferica alla specie a

SrotA =

3

2R ln

(ma

h2

)+

3

2R ln

(π1/38π2kBT

)+R lnVa +

3

2R (4.33)

Da cui, applicando l’analoga della Eq. 4.22, otteniamo l’energia rotazionale della specie:

Erot =3

2RT (4.34)

79

Dunque anche il questo caso, esattamente come nel caso traslazionale, avremo che

l’energia di dissociazione relativa al termine rotazionale e data da:

∆ErotAB =

3

2RT (4.35)

Con una procedura del tutto analoga a quella gia vista per l’entropia di dissociazione

traslazione (vedi Eqs. 4.25 e 4.36), otteniamo l’entropa rotazionale di dissociazione come

∆SrotAB = Srot

A + SrotB − Srot

AB

=3

2R ln

M

h2+

3

2R ln(π1/38π2kBT ) +R ln

(VaVbVab

)+

3

2R

(4.36)

In cui Va,Vb e Vab sono i volumi delle specie definiti dall Eq. 4.32 ed M e la massa ridotta

del sisetma farmaco-proteina. Il contributo rotazionale all’energia libera standard e

quindi dato dato da

∆GrotAB = ∆Erot

AB − T∆SrotAB

=3

2RT − 3

2RT ln

M

h2− 3

2RT ln(kBT )−RT ln

(29/2π7/2VaVb

Vab

)− 3

2RT

= −3

2RT ln

M

h2− 3

2RT ln(kBT )−RT ln

(29/2π7/2VaVb

Vab

)(4.37)

Analogamente a quanto visto per il contributo traslazionale, Eq. 4.30, anche il ternine

di rotazione relativo alla dissociazione del complesso, Eq. 4.37, e fortemente negativo.

Sfruttando il fatto che il volume/massa della proteina e in generale molto maggiore del

volume/massa del farmaco, possiamo assumere che Vab ' Vb ed M ' ma. La Eq. 4.37

puo quindi essere scritta come

∆GrotAB = −RT ln

[(mA

h2π1/38π2kBT

)3/2

Va

](4.38)

Come nel caso traslazionale, anche in questo caso, il contributo rotazionale alla dissoci-

azione, ∆GrotAB, puo essere stimato ponendo mA ' 5×10−25 kg e Va ' 5×10−28 m3, valori

tipici per un farmaco del peso molecolare di ca 300 Da, ottenendo una energia libera

80

rotazionale di dissociazione dell’ordine delle -10 kcal/mole che favorisce fortemente la

dissociazione del complesso, sempre a causa del contributo entropico, molto maggiore per

il sistema dissociato nel quale farmaco e proteina possono ruotare indipendentemente.

Combinando le due equazioni 4.37 e 4.28, otteniamo infine il contributo roto-traslazionale

all’energia libera di dissociazione del complesso farmaco proteina AB:

∆Grot−traslAB = −3RT ln

M

h2− 3RT ln(kBT )−RT ln

(V026π5

)−RT ln

(VaVbVab

)(4.39)

Complessivamente il ∆Grot−traslAB e’ dunque fortemente negativo, per un valore intorno

ai -20 kcal mole−1 per farmaci potenziale con masse comprese tra i 300 e 400 Da. Se

dunque i contributi al ∆G di dissociazione fossero solo quelli traslazionali e rotazionali

(cioe ∆Gel = ∆Gvib = 0), la formazione del complesso farmaco-proteina non avverrebbe

mai in condizioni fisiologiche. Da notare infine che i contributi rotazionali e traslazionali

teste’ calcolati sono esatti anche per soluzione diluite.

Con la Eq. 4.39 e con la Eq. 4.13 abbiamo per il momento stimato a partire dalle

conoscenze molecolari ben 6 termini degli otto termini che contribuiscono nella Eq. 4.11,

e cioe, l’energia elettronica di dissociazione (derivabile dal force field come vedremo nel

prossimo capitolo), l’entropia elettronica di dissociazione (che e’ nulla per stati fonda-

mentali non degeneri come avviene per la quasi totalita delle coppie farmaco-proteina),

l’energia roto-traslazionali di dissociazione pari a 3RT (positiva) e il termine entropico

dominante −T∆S roto-traslazionale di dissociazione, negativo e dell’ordine delle decine

di kcal/mole. Abbiamo visto che le energie roto-traslazionali sono sempre le stesse per

tutti le coppie farmaco-proteina. Abbiamo inoltre constatato come le entropie roto-

traslazionali siano anch’esse scarsamente influenzabili dalla struttura molecolare (esse

cambiano solo con il logaritmo della massa molecolare). Le costanti nella 4.39 possono

essere calcolate esattamente e dipendono, come abbiamo visto, solo dalle condizioni ter-

modinamiche (temperatura dello stato standard e dalle costanti universali h, kB, π, NA)

e dalla struttura del farmaco A, della proteina B e del complesso AB. Ripetiamo ancora

81

che l’equazione 4.39 e un risultato esatto in approssimazione classica che vale tanto per

le miscele di gas ideali, quanto per le miscele di A,B e AB in soluzione diluita.

Veniamo infine all’ultimo termine che ci resta da determinare nella Eq. 4.11, l’energia

libera vibrazionale di dissociazione ∆GvibAB. Al solito procediamo a partire dalla funzione

partizione vibrazionale definita nella 3.30. L’energia libera vibrazionale per una specie

generica e:

Gvib = −RT ln3M−6∏i=1

e− hνi

2kbT · 1

1− e−hνikbT

(4.40)

dove M e il numero di atomi che compongono la specie in questione. Abbiamo quasi

sempre che: MA MB, ovvero la proteina B ha molti piu atomi del farmaco A. Dalla

4.40 abbiamo banalmente che:

Gvib =3M−6∑i=1

gvibi =

3M−6∑i=1

Evibi − T∆Svib

i (4.41)

con gvibi definito come

gvibi = NA

hνi2

+RT ln

(1− e−

hνikbT

). (4.42)

La quantita hνi2

e’ l’energia vibrazionale di punto zero del grado di liberta i-esimo (vedi

Eq. 3.31). Nella Eq. 4.42 si e’ tenuto conto che R = KBNA, dove NA e’ il numero

di Avogadro. La espressione 4.42 ci dice che l’energia libera vibrazionale di una specie

generica si ottiene sommando sulle 3M− 6 frequenze normali le energie libere per mole

dei singoli gradi di liberta vibrazionali gvibi (in cui l’indice i e’ riferito alla frequenza di

vibrazione i-esima). La variazione dell’energia di punto zero (per mole) per effetto della

dissociazione e’ data

∆E0vib = NA

3MA−6∑i=1(1∈A)

hνi2

+NA

3MB−6∑i=1(1∈B)

hνi2−NA

3(MA+MB)−6∑i=1(1∈AB)

hνi2. (4.43)

Tale energia non dipende dalla temperatura ma solo dalle caratteristiche molecolari delle

specie e generalmente viene incorporata nella energia elettronica (vedi Eq. 3.33).

82

Valutiamo ora la parte eccedente (per mole) sull’energia di punto zero per grado

vibrazionale, cioe:

gvibi = RT ln

(1− e−

hνikbT

)(4.44)

L’entropia svibi del grado di liberta vibrazionale i-esimo puo essere calcolata derivando

la Eq. 4.44 rispetto alla temperatura tramite l’analoga della Eq. 4.18 riferita alla parte

vibrazionale. Otteniamo

∂gvibi

∂T= −svib

i = R ln

(1− e−

hνikbT

)+RT

e− hνikbT

1− e−hνikbT

·(− hνikbT 2

)svibi = −R ln

(1− e−

hνikbT

)+

NAhνi

ehνikbT − 1

· 1

T(4.45)

Dalla Eq. 4.45, sommando sui gradi di liberta vibrazionali, otteniamo l’espressione per

l’entropia vibrazionale per mole per una specie generica

Svib = −R3M−6∑

i

ln

(1− e−

hνikbT

)+

3M−6∑i

NAhνi

ehνikbT − 1

· 1

T(4.46)

Utilizzando le Eqs. 4.44 e 4.45 nella relazione termodinamica e = g+ ts, otteniamo una

espressione analitica per l’energia vibrazionale riferita al grado i-esimo:

eivib =NAhνi

ehνikbT − 1

(4.47)

Questa e l’energia termica di tipo vibrazionale per grado di liberta per mole, cioe l’energia

in eccesso rispetto all’energia di punto E0 = NAhνi/2 per mole. Vediamo cosa succede

se assumiamo kbT hνi cioe tutti i gradi di liberta vibrazionali assimilabili a gradi

di liberta classici. Quando siamo in queste condizioni, l’energia classica diviene uguale

all’energia quantistica e possiamo trascurare gli effetti dell’energia di punto zero.8 Per

kBT hνi la 4.47 diviene

eivib =NAhνi

1− 1 + hνikbT

= RT (4.48)

8Nella Figura 3.1 quando kBT hν l’energia classica per grado per molecola e’ kBT indipendente-

mente da quanto sia l’energia di punto zero. Dunque ad alta temperatura l’energia classica rappresenta

tutta l’energia vibrazionale.

83

Per (3M− 6) gradi di liberta che si comportano classicamente avremo che

Evib = (3M− 6)RT (4.49)

Dalla Eq. 4.49 (riferita ad una specie generica) possiamo speditamente calcolare l’energia

vibrazionale di dissociazione in eccesso (classica) per la reazione AB A+B , ottenendo

∆EvibAB(classico) = (3MA − 6)RT + (3MB − 6)RT − [3(MA +MB)− 6]RT

= −6RT (4.50)

L’energia vibrazionale del complesso AB e’ superiore a quella dello stato separato. Come

Figura 4.2: Perdita di gradi di liberta roto-traslazionali per effetto della formazione del

complesso farmaco-proteina.

mostrato in figura 4.2, cio avviene perche nel complesso si hanno 6 gradi vibrazionali

in piu e ciascuno di questi gradi vibrazionali contribuisce positivamente all’energia vi-

brazionale. I sei gradi in piu nel complesso sono dovuti alla contemporanea soppressione

84

di sei gradi di liberta roto-traslazionali per effetto della complessazione. In AB com-

plessati ci sono infatti solo 6 gradi di liberta roto-traslazionale (nel complesso A e B

ruotano e traslano in maniera solidale) mentre in A e B separate i gradi di liberta roto-

traslazionali sono 12. I sei gradi “persi” per effetto della associazione si sono trasformati

in gradi di liberta vibrazionali aumentando l’energia vibrazionale di eccesso del comp-

lesso. Ne consegue che, dal punto di vista vibrazionale, la formazione del complesso e’

energeticamente sfavorita di circa 1 kcal/mole (6RT ) in condizioni standard.

Per quanto riguarda l’entropia vibrazionale riferita ad una specie generica, quando

kBT hνi dalla Eq. 4.46 deriviamo che

Svib(classico) = −R3M−6∑

i

ln

(hνikbT

)+ (3M− 6)R (4.51)

L’entropia vibrazionale di una specie generica e’ dunque tanto piu grande quanto mag-

giore e’ il numero di gradi di liberta 3M− 6 presenti nel sistema. Essa e’ dunque

una quantita strettamente collegata alla struttura molecolare delle specie. L’entropia

vibrazionale della proteina e’ per esempio molto maggiore di quella del farmaco.

L’entropia vibrazionale di dissociazione sara quindi data dalla

∆SvibAB = R

3(MA+MB)−6∑i

ln

(hνikbT

)−R

3MA−6∑i

ln

(hνikbT

)−R

3MB−6∑i

ln

(hνikbT

)+

+ (3MA − 6)R + (3MB − 6)R− [3(MA +MB)− 6]R (4.52)

Come detto, nel complesso, ci sono 6 frequenze vibrazionali in piu, dovute alle rotazioni

e traslazioni vibrazionali del farmaco nel sito attivo della proteina (vedi Figura 4.2).

Se si assume che tutte le altre frequenze di vibrazione restino invariate nel processo di

formazione del complesso AB, abbiamo che:

∆SvibAB = R

6∑i∈AB

ln

(hνikbT

)− 6R (4.53)

dove le ν1, ν2, ..., ν6 sono le frequenze roto-traslazionali di vibrazione del farmaco nella

proteina. Usando per esse dei valori caratteristici tipici cm−1 o al massimo decine di

85

cm−1, la Eq. 4.52 ci dice che dal punto di vista entropico, la formazione del complesso e’

favorita. Cio avviene perche i gradi di liberta vibrazionali complessivi aumentano pas-

sando dallo stato dissociato, al complesso. Combinando le equazioni 4.50 e 4.53 possiamo

valutare infine l’energia libera di dissociazione relativa al solo contributo vibrazionale

ottenendo

∆GvibAB = ∆Evib

AB − T∆SvibAB

= −RT6∑

i∈AB

ln

(hνikbT

)(4.54)

Come detto, le frequenze νi..ν6 sono dell’ordine di qualche cm−1 o al massimo decine di

cm−1. Dunque assumendo νi = 3 : 30 cm−1 ( e quindi hνikbT∼= 0.01 : 0.1), abbiamo che il

contributo entropico vibrazionale di dissociazione in condizioni standard si puo stimare

in

∆GvibAB ∼ 10 : 20 kcal mole−1 (4.55)

Il contributo vibrazionale positivo contribuisce pertanto a compensare in parte il con-

tributo negativo roto-tralsazionale complessivamente stimato intorno alle -20 kcal mol−1.

Curiosamente, l’energia libera vibrazionale di dissociazione e’ maggiore (cioe il complesso

e piu favorito) se le frequenze di vibrazione nel sito sono piu piccole, ovvero se l’ampiezza

di oscillazione nel sito e maggiore. Nel computo complessivo di ∆GvibAB = ∆Evib

AB−T∆SvibAB

il termine entropico gioca quindi un ruolo piu importante rispetto al termine energetico

vibrazionale (di segno opposto) sempre ugiale a −6RT . Cio accade perche per il teore-

ma di equipartizione l’energia vibrazionale di dissociazione e sempre uguale ad −RT per

grado di liberta indipendentemente dalla frequenza di vibrazione nel sito attivo. Questa

situazione e illustrata in Figura 4.3

E chiaro, a questo punto, che il termine guida per l’interazione farmaco-proteina

corrisponde al termine energetico, Eq. 4.13, che deve includere l’energia vibrazionale di

86

Figura 4.3: Nel complesso a sinistra, il farmaco (pallina rossa) e’ strettamente legato (fre-

quenza di librazione nel sito maggiore) al recettore (in blu). L’ampiezza di oscillazione nel

sito attivo e’ indicato con Rv. A destra il farmaco puo oscillare con ampiezza maggiore nel

sito. I due complessi sono caratterizzati dalla stessa energia elettronica di binding, ∆Eel, dalla

stessa energia vibrazionale di binding (stabilita dal livello termico tratteggiato corrispondente

a kBT ), ma da una diversa entropia vibrazionale che favorisce il complesso a destra.

punto zero, Eq. 4.43:

∆GelAB = ∆Eel

AB −6∑i

NAhvi2

(4.56)

Nella Eq. 4.56, le quantita EA, EB ed EAB, tali che ∆EAB = EA + EB − EAB, cor-

rispondono alle energie fondamentali (per mole) del farmaco A, della proteina B e del

complesso AB nella loro configurazione di minima energia ed NA e il numero di Avo-

gadro. Per la correzione di punto zero abbiamo assunto, come in precedenza, che le

frequenze di vibrazione proprie di farmaco e proteina isolati cambino poco a seguito

87

della formazione del complesso. Il termine vibrazionale di punto zero, che puo’ essere

valutato intorno alle frazioni di kcal mol−1 per frequenze librazionali del farmaco nel

sito attivo di pochi o poche decine di cm−1, costituisce una correzione piccola rispetto

al termine elettronico propriamente detto, ∆EAB = EA + EB − EAB.

Possiamo a questo punto fornire una stima generica per l’energia libera di disso-

ciazione farmaco-proteina, inserendo i risultati relativi ai singoli contributi valutati, e

cioe Eq. 4.56 (energia libera di dissociazione elettronica), Eq. 4.39 (energia libera

di dissociazione roto-traslazionale) ed infine Eq. 4.54 ( energia libera di dissociazione

vibrazionale). Prima di procedere alla valutazione della stima conviene manipolare ulte-

riormente il termine vibrazionale. Come vedremo infatti, l’energia libera di dissociazione

complessiva non dipende dalle masse del sistema ne tantomeno dalla costante di Planck

h. Questo deve avvenire perche: 1) le masse entrano solo nell’energia cinetica (il termine

cosiddetto ideale dipendente dai momenti) e gli integrali relativi ai momenti si semplifi-

cano nella valutazione del rapporto delle funzioni di partizione Kdiss = qAqB/qAB; 2) ogni

atomo porta un fattore h alla funzione di partizione complessiva ed il numero di atomi

e lo stesso a destra o sinistra della reazione di formazione del complesso. Riscriviamo

dunque ∆GvibAB tenendo conto che

νi =1

2π

(Ki

M

)1/2

(4.57)

in cui Ki rappresenta la costante di forza relativa alla frequenza i-esima di librazione

del farmaco nel sito attivo ed M e la massa ridotta del sistema farmaco-proteina.

Utilizzando la Eq. 4.57, la Eq. 4.54 diviene

∆GvibAB = −RT

6∑i∈AB

ln

(h(KiM

)1/2

2πkbT

)

= 3RT lnM

h2+ 6RT ln(2πkBT )−RT

6∑i∈AB

ln(K1/2i ) (4.58)

88

Se definiamo

K6a =

6∏i∈AB

Ki (4.59)

in modo tale che

RT6∑

i∈AB

ln(K1/2i ) = 3RT ln(Ka) (4.60)

la Eq. 4.58 si trasforma nella

∆GvibAB = 3RT ln

M

h2+ 6RT ln(2πkBT )− 3RT ln(Ka) (4.61)

Utlizzando le Eqs. 4.56, 4.39 e 4.61 otteniamo infine una stima compatta per l’energia

libera standard di dissociazione farmaco-proteina:

∆GAB = ∆GelAB − 3RT ln

(Ka(lal0)

kBT

)(4.62)

nella quale abbiamo definito la quantita l relativa la farmaco (che ha le dimensioni di

una lunghezza)

lal0 = (πV0Va)1/3 (4.63)

La quantita la = V1/3a e’ una lunghezza caratteristica del farmaco che puo essere cal-

colata, nota la sua struttura chimica, a partire dalle Eqs. 4.63 e 4.32. La costante

l0 = 16611/3 = 11.84 A, corrisponde alla radie cubica del volume standard riferito alla

concentrazione 1M. In tabella 4.1 riportiamo i valori di la calcolati per alcuni dei farmaci

piu conosciuti. La quantita Ka misura la costante di forza geometricamente mediata sui

sei gradi di liberta librazionali del farmaco nel sito attivo ed e strettamente collegata

all’energia di binding stessa, ∆EelAB, ed dal volume del sito attivo. Se si approssima

l’energia potenziale farmaco proteina con un potenziale armonico, si puo assumere che

K =8∆Eel

AB

l2site(4.64)

89

Drug MM la (lal0)1/2

Arformoterol 348. 36.19 20.7

Aripiprazole 450. 68.06 28.38

Aspirin 184. 15.72 13.64

Fluticasone 448. 39.62 21.65

l-Thyroxine 1285 24.30 18.308

Nexiam 348 38.67 21.39

Pregabalin 161 16.54 13.99

rosuvastatin 487 48.98 24.08

Salbutamol 242 26.57 17.73

Sitagliptin 408 47.54 23.72

THC 316 39.05 21.50

Viagra 478 52.24 24.87

Tabella 4.1: valori di la (in A) per alcuni tra i farmaci piu conosciuti. MM indica il peso

molecolare del farmaco.

in cui lsite rappresenta il raggio medio nel sito di legame del recettore. L’approssimazione

implicita nella Eq. 4.64 e illustrata nella Figura 4.4.9 Inserendo la Eq. 4.64 nella Eq.

4.62, otteniamo

∆GAB = ∆GelAB − 3RT ln

(8∆Eel

AB

kBT

)− 3RT ln

(lal0l2site

)(4.65)

L’ultimo termine nella Eq. 4.65 costituisce una penalita per il binding quando (l0la)1/2 >

lsite, cioe quando le dimensioni del farmaco eccedono quelle del sito di legame. In gen-

erale, in prima approssimazione, possiamo assumere che farmaco abbia le dimensioni

ottimali per il sito di binding, ovvero (lal0)1/2 ' lsite. Possiamo inoltre trascurare l’ener-

9E importante sottolineare che la relazione 4.64 e valida solo nel caso in cui ∆EelAB sia positiva, cioe

quando il binding sia favorito dal punto di vista energetico.

90

Figura 4.4: Relazione di proporzionalita tra K ed energia di binding; l’interazione farmaco

proteina e approssimata con una parabola. lsite rappresenta il raggio del sito.

91

gia vibrazionale di punto zero del farmaco (vedi Eq. 4.56 ) nel sito di binding ponendo

∆GelAB ' ∆Eel

AB, ed ottenendo infine una espressione approssimata per l’energia libera

standard di dissociazione

∆GAB = ∆EelAB − 3RT ln

(8∆Eel

AB

kBT

)(4.66)

Nell’espressione finale del ∆G di dissociazione cui siamo pervenuti Eq. 4.66, che ci

permette di calcolare la costante di equilibrio per la reazione AB A+ B utilizzando

la 4.2, compare il solo termine relativo all’energia elettronica di binding. La dipendenza

della energia libera di dissociazione dal volume dello stato standard compare nel solo

termine di penalita −3RT ln[(lal0)1/2/lsite] che,nel passaggio dalla Eq. 4.65 alla Eq.

4.66, abbiamo in prima approssimazione trascurato assumendo (lal0)1/2 ' lsite. L’energia

libera di dissociazone corrisponde pertanto alla somma di due termini, dei quali il primo,

di natura strettamente energetica, coincide con l’energia di dissociazione ∆EAB. Il

secondo termine e di natura entropica e penalizza la formazione del complesso tanto piu

quanto piu la profondita del minimo del potenziale di legame farmaco-proteina eccede

il livello energetico termico kBT . L’importanza del secondo termine nella valutazione

complessiva dell’energia libera di dissociazione e mostrato nella Figura 4.5 in cui ∆G,

Eq. 4.66, e comparato all’energia di dissociazione ∆EelAB. Come si vede, il termine

entropico piu che compensa il termine energetico per bassi valori di ∆EAB. Inoltre,

in accordo alla Eq. 4.66, farmaci con attivita nanomolare (ovvero con energia libera

standard dell’ordine delle 11-12 kcal mol−1) devono essere caratterizzati da energia di

binding dell’ordine delle 18-20 kcal mol−1.

Come detto, l’energia di binding ∆EelAB viene valutata facendo ricorso alle parametriz-

zazione descritta del capitolo 1 riguardante il force field (con qualche accortezza per il

trattamento de solvente in soluzione diluita, come vedremo piu in dettaglio nel capi-

tolo seguente). Tale parametrizzazione rappresenta, in linea di principio, la superficie

di Born-Oppenheimer relativa alla energia di stato fondamentale del sistema farmaco-

proteina. Quando l’energia di interazione ∆EAB viene calcolata con il force field, il

92

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30

∆E (kcal mol-1

)

-4 -4

-2 -2

0 0

2 2

4 4

6 6

8 8

10 10

12 12

14 14

16 16

18 18

20 20

22 22

24 24

26 26

28 28

30 30

∆ G

(k

cal

mo

l-1)

∆E∆ G

Figura 4.5: Andamento dell’energia libera in funzione dell’energia di binding ∆E in accordo

alla Eq. 4.66

suffisso el non si riporta. L’espressione approssimata 4.66 per la differenza di energia

libera di dissociazione mette in evidenza l’importanza del calcolo dell’energia di bind-

ing per la determinazione delle costanti di equilibrio relative tra farmaci su uno stesso

target. Tale espressione e di fatto alla base delle moderne computazionale basate sul

docking, cioe sulla valutazione e confronto delle energie di binding per la identificazione

e lo sviluppo di farmaci. La valutazione dell’energia di binding ∆EAB a partire dai

force fields piu accurati, e’ dunque di fondamentale importanza per la progettazione

razionale di farmaci. Nel capitolo successivo accenneremo al calcolo di questa quan-

tita in soluzione diluite nel quale il solvente viene trattato come un mezzo dielettrico

continuo. Come deve essere calcolata l’energia di dissociazione ∆EAB? Sulla base di

quanto detto, dovremmo calcolare tutte le energie al minimo. L’energia in eccesso

vibrazionale e’ stata infatti assunta di tipo armonico.10 Per calcolare l’energia di disso-

10Nell’approssimazione armonica si assume che la configurazione del sistema si discosti leggermente da

93

ciazione nella Eq. 4.66, bisogna minimizzare tutte le specie e valutare poi la differenza

∆EAB = EA(min) + EB(min) − EAB(min). Naturalmente, l’approssimazione armonica

potrebbe essere ragionevole per farmaci rigidi, ma non lo e’ sicuramente per le proteine.

Dunque la Eq. 4.66, basata sul calcolo di ∆EAB, e’ senz’altro istruttiva ed importante

da un punto di vista qualitativo ma non costituisce una via quantitativa del tutto affid-

abile per la predizione di costanti di equilibrio. Cionondimeno essa puo’ essere utilizzata

con profitto in un primo screening di farmaci. Nell’ultimo capitolo accenneremo ad al-

cune moderne tecniche di calcolo dell’energia libera di dissociazione basate su metodi di

simulazione al computer di sistemi descritti a livello atomistico sia per quanto riguarda

le specie farmaco e proteina, sia per quanto concerne il solvente.

4.2 Cenni sul calcolo dell’energia di interazione farmaco-

proteina in solvente implicito

Come abbiamo visto, nella Eq. 4.66 abbiamo ottenuto un espressione della energia libera

di dissociazione faramco-proteina in termini della sola energia di binding

∆EAB = E0AB − E0

A − E0B (4.67)

Tale termine si valuta attraverso i force-field che parametrizzano la superficie di Born-

Oppenhainer dello stato fondamentale del mio sistema. Come abbiamo accennato in

precedenza, il solvente nel quale avvengono le reazioni chimiche, ha un ruolo pressoche’

nullo11 nella determinazione dei contributi roto-traslazioni all’energia di binding ed un

ruolo limitato nella determinazione del contributo entropico vibrazionale (Eq. 4.53).

una configuarzione di equilibrio di minima energia. L’energia per configuarzioni che differiscono da quel-

la di equilibrio la si ottiene svuluppando attorno all’unico minimo energetico. Questa approssimazione

noi l’abbiamo implicitamente utilizzata per tutta le specie A,B e AB.11In realta’, a causa dell’accoppiamento vibro-rotazionale, vi e’ un effetto solvente, in genere piccolo,

anche nei contributi roto-traslazionali all’energia libera di binding

94

Figura 4.6: Ciclo termodinamico per la determinazione dell’energia libera di binding in sol-

vente. Tale quantita puo’ essere valutata lungo il ciclo evidenziato in viola in figura, somman-

do l’energia libera di desolvatazione di A e B separati, l’energia di binding in fase gassosa e

l’energia di solvatazone del complesso AB.

L’energia di interazione inter ed intra-molecolare e, al contrario, fortemente mediata

dal solvente. La Eq. 4.67 deve pertanto essere in realta intesa come:

(E0AB)solv − (E0

A)solv − (E0B)solv = (∆E0

AB)solv (4.68)

Dove con (E0A,B,AB)solv, si intende l’energia della specie nel solvente. La reazione

chimica di formazione del complesso farmaco-recettore AB puo essere vista attraverso

un ciclo termidinamico, come mostrato in Figura 4.6

∆EsolvAB = ∆Edesolv

A +∆EdesolvB +∆E0

AB +∆EsolvAB (4.69)

Poiche G e’ una funzione di stato, la differenza (∆E0AB)solv nella Eq. 4.68 puo essere

valutata lungo il percorso evidenziato dalla freccia viola in Figura 4.6. Al termine

95

energetico calcolato in fase gassosa ∆E0AB occorre dunque aggiungere i termini ∆Edesolv

A

e ∆EdesolvB cioe le energie (libere) necessarie per desolvatare i componenti separati A e B

e portarli in fase gassosa ed il termine ∆EsolvAB , cioe l’energia di solvatazione del complesso

AB. Datosi che, per esempio, ∆EdesolvA = Egas

A −EsolvA , le energia libere di desolvatazione

di una data specie sono evidentemente identiche a le corrispondenti energie libere di

solvatazione cambiate di segno. Dunuqe avremo che l’Eq. 4.69 puo’ essere scritta come

∆EsolvAB = ∆Egas

AB −∆EsolvA −∆Esolv

B +∆EsolvAB (4.70)

Il calcolo diretto dell’energia di binding in solvente (cioe’ senza seguire la curva viola

di Figura 4.6) passa necessariamente attraverso il calcolo delle funzioni partizione delle

singole specie in presenza di solvente:

QA =

∫e−βH(rA,rsolv,pA,psolv)drAdrsolvdpAdpsolv (4.71)

dove rA , pA e rsolv, psolv sono coordinate e momenti rispettivamente della specie A e

delle molecole di solvente. L’Hamiltoniano del sistema (per la specie A solvatata nel

caso specifico) e definito come

H(rA, rsolv,pA,psolv) =∑k∈A

p2K

2mk

+∑k∈solv

p2k

2mk

+ V (rA) + V (rsolv) + V (rA, rsolv) (4.72)

ed espressioni simili per QB e QAB. Il calcolo della 4.71, a causa del termine di accop-

piamento V (rA, rsolv), e difficilissimo. Per evitare di avere che fare con l’integrazione

sulle coordinate del solvente, si assume in generale che il solvente si comporti come un

mezzo continuo (approssimazione di solvente implicito) che scherma l’interazione elet-

trostatica tra gli atomi delle specie A,B, e AB (considerate a livello atomistico) per

effetto della sua natura ”dielettrica”. Questa schermatura rappresenta un effetto ter-

modinamico medio e si riduce di fatto ad un numero (uno scalare) rappresentato dalla

costante dielettrica ε. Nella formulazione generale dell’approccio a solvente implicito, le

interazioni elettrostatiche tra due atomi con coordinate ri ed rj sono date dalla formula:

96

Eelih =

qiqjrijεij(rirj)

(4.73)

dove la costante dielettrica εij(rirj) varia da 1 a 80 a seconda che ri ed rj si riferiscano

ad atomi molto poco esposti (interni) o fortemente esposti al solvente. Esistono di-

verse tecniche per la valutazione della costante effettiva εij(rirj) che di fatto modula il

potenziale elettrico in funzione dell’esposizione al solvente (cioe al mezzo dielettrico).

Citiamo tra queste tecniche il metodo di Poisson ed il metodo di Born generalizzato. Nel

metodo Poison si risolve l’equazione di Poisson in ogni punto dello spazio a disposizione

del farmaco e della proteina per valutare, data una conformazione di essi, il potenziale

elettrostatico su una carica sonda. Questo consente di valutare le energie elettrostatiche

in gioco per la configurazione data. In pratica l’equazione di Poisson e:

Figura 4.7: Esempio di superficie di Van Der Waals o Solvent excluding surface (SES) per

il complesso 1-butyl-3-metilimidazolo tetrafluoroborato. All’interno della SES la costante

dielettrica e’ unitaria. All’esterno essa vale 80.

97

∇r· ε(r)∇rΦ(r) = −Aπρ(r) (4.74)

Nella equazione sopra, la quantita ρ(r) =∑N

i=1 qiδ(r − ri) rappresenta la densita di

carica elettrica12 qi ed ri sono le cariche atomiche e le coordinate degli atomi. l’equazione

di Poisson viene risolta come detto per una singola conformazione del complesso AB.

Dunque le quantita ρ(r), qi ed ε(r) sono note. Le qi sono infatti le cariche atomiche

specificate nel force field, mentre per ε(r) si pone ε = 1 all’interno di farmaco e proteina

ed ε = 80 fuori. In termini piu precisi, l’interno del sistema (a costante dielettrica

unitaria) e definito come lo spazio racchiuso nel volume di Van Der Waals che si ottiene

assumendo che gli atomi siano sfere con un raggio pari al raggio di Van Der Waals (vedi

Figura 4.7) L’equazione di Poisson viene risolta per il potenziale elettrostatico Φ(r).

Tale campo scalare si determina risolvendo l’equazione differenziale del secondo ordine

alle derivate parziali, Eq. 4.74, per via numerica, note ρ(r) ed ε(r). L’energia di una

data configurazione sara’ data E =∑

i qiΦ(ri).

12La funzione δ(r− ri) e detta funzione Delta di Dirac. Essa e’ tale che∫f(r)δ(r− ri)dr = f(ri)

98

Indice analitico

Angolo Diedro, 14

Approssimazione di Born OppenheimerK

Approssimazione di, 53

Approssimazione di Stirling, 61

Bagno Termico

Insieme Canonico, 43

Barriera Torsionale, 16

Bending

Potenziale di, 12

Boltzmann

Costante di, 42

Born Oppenheimer

Approssimazione di, 53

Costante di Boltzmann, 42

Costante di Equilibrio Chimico, 64, 67

Costante di Forza

Nel Potenziale di Stretching, 9

Costante di Planck, 49

Costante Dielettrica, 97

Densita di Carica Elettrica, 98

Dipolo Indotto, 24

Dirac

funzione Delta, 98

Energia

di un Oscillatore Armonico, 56

di Van Der Waals, 23

Energia Cinetica, 38, 49

Rotazionale, 57

Energia Coulombiana, 27

Energia di Solvatazione, 95

Energia Elettronica, 74

Energia Elettrostatica

In Solvente Implicito, 97

Energia Libera, 4, 46

dello Stato Standard, 64

Energia Libera di Dissociazione, 88

Energia Libera di Eccesso, 50

Energia Media, 38

Insieme Canonico, 47

Energia Potenziale, 49

Energia Roto-traslazionale di Dissociazione,

81

Energia Rotzionale, 57

99

Energia Torsionale, 17

Energia Vibrazionale, 56

di Punto Zero, 56

Energia Vibrazionale di Dissociazione, 82

Energia Vibrazionale di Punto Zero di Dis-

sociazione, 82

Entropia, 33, 39

Definizione nell’Insieme Microcanocico,

42

Rotazionale, 79

Roto-trasalzionale, 73

Traslazionale, 76

Equazione di Poisson, 98

Equilibrio Chimico

Per una Miscela Ideale di Gas, 59

Equilibrio Termodinamico, 30, 35

Etano, 14

Force Field, 2

Frequenza di Stretching, 66

Funzione di distribuzione, 37

Funzione di distribuzione Canonica, 45

Funzione di distribuzione Microcanonica,

41

Funzione di partizione

Gas ideale, 50

Funzione di partizione Canonica, 45, 46

Funzione di partizione di Eccesso, 50

Funzione di partizione Elettronica, 53

Funzione di Partizione Microcanonica, 42

Funzione di partizione Molecolare, 51, 58

Fattorizzazione, 53

Funzione di partizione Rotazionale, 57

Funzione di partizione Traslazionale, 58

Funzione di partizione Vibrazionale, 55

Funzione Hamiltoniana, 49

Gas Perfetto, 48

Funzione di partizione, 50

Hamiltoniano, 45, 49

Insieme Canonico, 43

Insieme Microcanonico, 40

Interazioni Dispersive, Repulsive, 22

Interazioni Elettrostatiche, 22

Legame ad Idrogeno, 27

Lennard-Jones

Potenziale di, 23

Media Statistica, 5, 38

Momento della Quantita di Moto, 39

Momento di Inerzia, 58

Morse

Funzione di, 9

Oscillatore Armonico, 56

100

Poisson

Equazione di, 98

Potenziale Bonded, 6

Potenziale bonded, 5

Potenziale bonded per sistemi biologici, 6

Potenziale Chimico, 60

Potenziale di Lennard Jones, 22

Potenziale Elettrostatico, 25

In Solvente Implicito, 98

Potenziale Non Bonded, 20

Potenziale non bonded, 6

Potenziale Torsionale, 14

Potenziale Torsonale

Configurazioni Sfalsate, Eclissate, 17

Primo Principio, 31

Formulazione Microscopica, 38

Processo Spontaneo, 33

Prodotto scalare tra vettori, 13

Prodotto Vettoriale, 14

Secondo Principio, 33

Formulazione Microscopica, 39

Sistema Isolato, 33, 40

Solvente Implicito, 5, 94

Stato Standard, 75

Stirling

Approssmazione di, 61

Stretching

Potenziale di, 8

Frequenza di vibrazione, 10

Temperatura

Definizione Statistica, 42

Torsione

Definizione di , 14

Valor Medio, 5

Van Der Waals

Energia di, 23

Superficie di , 98

101