Computing structural chemistry and biology
Chemistry and biology with computers
http://www.nobel.se/chemistry/laureates/1998/illpres/index.html
Why computing chemistry Snížené náklady na laboratorní experimenty Zvýšená bezpečnost při práci s toxickými látkami Předpověď chování nových látek
Why computing chemistry Interpretace experimentálních dat - NMR refinement pomocí MD, mapy elektronových hustot v RTG Pohled „za experiment“ - kam experiment nedosáhne... 50μs
4ns
NMR vnitřní pohyby
10ms
NMR nevidí...
MD MD nevidí...
Why computing Biology Rational drug design ($$$) Založený na MECHANISMU reakce: Kvantově-mechanické studium reakční cesty, popis tranzitního stavu…
Založený na RECEPTORU: Známe X-ray strukturu enzymu, komplexu. Návrhy nových inhibitorů, docking…
Založený na LIGANDU: Neznáme exp. strukturu receptoru. Konformační prohledávání ligandů, výpočet elektron. vlastností, farmakofor, QSAR
Metody výpočtu energie Kvantová mechanika
Molekulová mechanika
formaldehyd vázaný vodíkovou vazbou k fluorovodíku
Energie vypočtená řešením Schrodingerovy rovnice
Energie vypočtená z mechanického modelu Energie pnutí
Metody molekulového modelování
Ab initio kvantová mechanika
Semiempirická kvantová mechanika
Molekulová mechanika
100 atomů
1000 atomů
100 000 atomů
Exaktní řešení Schrodingerovy rovnice
Aproximativní řešení Schrodingerovy rovnice
roste výpočetní náročnost
Užití empiricky odvozené potenciálové funkce
použití empirických parametrů
Hyperplocha potenciální energie Potential Energy (hyper)Surface (PES)
tranzitní stav
Energie = funkce struktury
Globální minimum
Energie
Lokální minimum
„reakční“ bariéra
Reakční koordináta
Souřadné systémy - Methan
Dvě vnitřní souřadnice:
15 kartézských
C -- H … 1.09 Å
souřadnic
H—C—H …109.5 deg
Alanin
O
CH3
O=C-NH-CH-C-NH-CH3 CH3
ACE
ALA
NME
9
5 1
Kartézské souřadnice: X, Y, Z
Počet atomů N Počet souřadnic 3N
Vnitřní souřadnice Torzní uhly fixní valenční uhly
2
6
4 3
7
8
10 11
12
Hyperplocha potenciální energie (PES) Ramachandranova mapa
β struktura
2D plocha: torzní úhly ϕ a φ v proteinech
C-N-Cα-C
N-Cα-C-N φ ϕ
Cα
pravotočivá α šroubovice
αL ,γ-obrat
Quantum chemistry
Born-Oppenheimerova aproximace oddělení vlnové funkce pro elektrony a pro jádra
hmotnost jádra >> hmotnost elektronu ⇒ eliminace elektronů ⇒ potenciální energie (PES) je funkcí polohy jader
Kvalitativní MO teorie 2 H • → H2
2 1s atomové orbitaly (AO) (každý na 1 atomu H)
2 molekulové orbitaly (MO): antivazebný vazebný (obsazen 2 e)
Molecular Orbitals • Bonding orbitals - concentrate electron density between atoms – H2 σg (1sa + 1sb)
• Antibonding orbitals - nodal plane in electron density between atoms – H2 σu* (1sa - 1sb)
http://hyperphysics.phy-astr.gsu.edu/hbase/molecule/hmol.html
Homonuclear Diatomics
http://science.widener.edu/svb/molecule/mo.html
http://xbeams.chem.yale.edu/~batista/vvv/node33.html
MO orbitaly lineární AH2 AO (fragmentové orbitaly) → MO
H–A-H
2pz 2py 2px 2s 4 AO atomu A
2 AO atomu H
MO orbitály lomená H2A z
y x
C2 C2v bodová grupa:
A H
H
dvojčetná osa symetrie 2 roviny symetrie v ose C2
MO orbitaly lomené H2A 4 AO na atomu A
b1
a1 2px
2s
a1
b2 2py
2 AO dvou vodíků
2pz z
b2 y
a1
x
MO orbitaly H2A
Podmínky vzniku vazby: stejná symetrie fragmentálních orbitalů ⇒ co největší překryv malý rozdíl energie
Energie
6 AO → 6 MO
Polyatomic molecules • Linear combinations of atomic orbitals • MO’s delocalized over several atoms 2b2
2σu
4a1
3σg 1b1
1πu
1b2
3a1
1σu
2a1 2σg 1a1
1σg H2 O
BeH2
Electron Density • In quantum mechanics, there are two ways of describing electron location and energy Ψ: wavefunction ρ: electron density ρ = Ψ2
• Electron density: probability of locating the electron
Visualization Methods Variables: x, y, z, Ψ (or Ψ2) Contour plots
Dot plots 2σg
Surface plots 2σg
science.uwaterloo.ca/~cchieh/cact/c120/pibond.html
http://www.shef.ac.uk/chemistry/orbitron/AOs/1s/e-density.html
• Positively charged nucleus attracts electron density ρ varies with distance from nucleus ρ depends on distance, not direction • isodensity surfaces ("same density" surfaces) map of points where ρ has a particular value
http://academic.reed.edu/chemistry/roco/Density/
Electron Density
Electrostatic Potential • Visualizes electronic charge distributions • Map electronic potential onto total electron density surface using colors to represent different levels of potential
Frontier Orbitals • In a chemical reaction, the most important interaction is between the HOMO-LUMO pair • Least energetically bound ∴ most readily available for interaction with other molecules • Nucleophilic Frontier Density: measure of the susceptibility of the substrate to attack by an electrophile; based on electron distribution of active orbitals near the LUMO • Electrophilic Frontier Density: measure of the susceptibility of the substrate to attack by an electrophile; based on electron distribution of active orbitals near the HOMO
Kvantitativní teorie MO Schrodingerova rovnice Hamiltonův operátor
vlnová funkce
energie systému
Born-Oppenheimerova aproximace hmotnost jádra >> hmotnost elektronu ⇒ jádra jsou fixní
Ab initio „z počátku“ Hartree-Fockova metoda
Hamiltonián příspěvek přitažlivé interakce mezi elektrony a jádrem
jádra jsou fixní
příspěvek repulze mezi elektrony
Hartree-Fockova metoda Self-consistent field (SCF) elektron v poli ostatních
interakce mezi všemi elektrony HF zjednodušení
Báze soubor funkcí, jejichž lineární kombinací jsou tvořeny molekulové orbitaly. LCAO.
MO ← AO ← bázové funkce Bázové funkce lze považovat za atomové orbitaly (minimální báze) matematické funkce, které dávají MO potřebnou flexibilitu (rozšířené báze)
STO Slater Type Orbital
elektronová hustota v atomu
H atom
1s orbital
C atom
1s, 2s, 3x 2p orbitaly
methan CH4
4x 1s pro H, 1s, 2s, 3x 2p pro C celkem 9 funkcí
1 funkce = 1 AO
minimální báze
STO → GTO STO: e-x -x2
GTO: e
1 STO nahrazena 3 GTO
(Gaussian Type orbital) methan CH4 4x 1s pro H, 1s, 2s, 3x 2p pro C = celkem 9 funkcí STO ⇒ 27 funkcí GTO
Minimální báze: STO-3G
STO x GTO aproximace STO lineární kombinací několika (např. 3) GTO ⇒ STO-3G
STO STO-3G 1s orbital atomu H
6-31G báze (split valence) 1s orbital “core” 1 bázová funkce - 6 GTO atom C
2s, 2p orbitaly (valenční): 2 bázové funkce: 3 GTO +1 GTO atom uhlíku 1s … 1 BF … 6 GTO 2s … 1 BF (3 GTO) + 1 BF (1 GTO) 2px … 1 BF (3 GTO) + 1 BF (1 GTO) 2py … 1 BF (3 GTO) + 1 BF (1 GTO) 2pz … 1 BF (3 GTO) + 1 BF (1 GTO) 9 BF … 22 GTO
Rozšířené báze zahrnutí polarizace ⇒ přidání polarizačních funkcí (přidání funkcí odpov. vyšším orbitalům)
atom H atom C
polarizace s orbitalu přidání p orbitalu přidání d orbitalu
6-31G* - přidává sadu d funkcí (např. pro C). 6-31G(d) 6-31G** - přidává sadu d funkcí a sadu p funkcí (např. pro H) Difuzní funkce 6-31+G(d). Core LANL2DZ.
Ab initio procedury jednoduchý výpočet energie (single point) výpočet energie, vlnové funkce na fixované geometrii
optimalizace geometrie výpočet energie, vlnové funkce
změna geometrie směrem k nižší energii
výpočty frekvencí a jiné předpověď frekvencí IČ a Raman, ověření geometrie získané optimalizací. Elektrostatický potenciál. Dipólový moment.
Optimalizace energie výpočet energie - analog single point výpočet sil (1. derivace energie na každém atomu) - tzv. gradient změna struktury (proti energet. gradientu) NE
výpočet energie Jsou splněna konvergenční kritéria?
ANO
Minimum energie nalezeno Zápis nové struktury
Omezení HF metod molekula v mezihvězdném prostoru (gas phase) HF nedává exaktní řešení Schrodingerovy rovnice. Nejlepší vlnová funkce dosažitelná HF metodami - HF limita. 2 elektrony v 1 orbitalu se nepohybují nezávisle - elektronová korelace. Rozdíl mezi HF limitní energií a exaktním řešením Schrodingerovy rovnice - korelační a relativistická energie.
Korelační energie rozdíl přesné energie a HF limity
Ecorr = Eexact - EHF Výpočet korelační energie: Moller-Plesset perturbační teorie 2.řádu = MP2 Moller-Plesset perturbační teorie 4.řádu = MP4 Teorie funkcionálu hustoty = Density Functional Theory (DFT)
Relativistická korekce • • • •
Einsteinova teorie relativity Těžké atomy Elektrony na vnitřních vrstvách Rychlosti blízké rychlosti světla
Metody funkcionálu hustoty Density functional methods (DFT) Energie je vyjádřena jako funkcionál hustoty elektronového pole ⇒ elektronová hustota v okolí molekuly DFT zahrnují podstatnou část korelační energie vlnová funkce je konstruována jiným způsobem (jiný druh orbitalů) ⇒ není to HF metoda!
Semiempirické metody
Kvantitativní MO teorie Schrodinger ⇓ Born-Oppenheimer ⇓ orbitalová approximace ⇓ LCAO approximace
Semiempirická aproximace zanedbání některých překryvových integrálů nahrazení některých integrálů funkcemi s empirickými parametry CNDO = Complete Neglect of Differential Overlap INDO = Intermediate Neglect of Differential Overlap NDDO = Neglect of Diatomic Differential Overlap MINDO/3 = Modified INDO MNDO = Modified NDO AM1 = Austin Model 1 PM3 = Parametric Model 3
fitování empirických parametrů
zanedbání překryvových integrálů
Možnosti semiempirických metod HF vlnové funkce molekulové orbitaly, energetické hladiny, náboje minima, tranzitní stavy, reakční koordináty silové konstanty, vibrační frekvence termodynamické vlastnosti řády vazeb, analýza σ a π orbitalů polarizabilita efekty solvatace
Omezení semiempirických metod zanedbání nebo parametrizace překryvových integrálů vnitřní orbitaly jsou zanedbány elektronová korelace je zahrnuta v parametrizaci
výpočetně méně náročné než ab-initio výpočetně náročnější než empirické metody
může vést k chybám
Software ab-initio: program GAUSSIAN program GAMESS program Spartan
semiempirické metody: MOPAC = Molecular Orbital Package program Spartan
Multiplicita 2S+1 • NO: 5+6=11=5.2+1. • CO: 4+6=10=5.2+0. • OO: 6+6=12=5.2+2.
2x1/2+1=2 2x0+1=1 2x1+1=3
Molekulová mechanika
Molekulová mechanika Empirical force-field Theoretical conformational analysis Strain energy calculation Westheimer-Hendrickson-Wiberg
Molekulová mechanika studium biologických makromolekul (tisíce atomů) umožňuje sledovat vliv okolí (solvent) atomy = sférické částice s bodovým nábojem interakce jsou založeny na klasické mechanice: tělesa (atomy) spojená pružinami (vazbami)
Silové pole (force field): rovnice silového pole soubor parametrů
Anatomie silového pole Energie(potenciální) : energie pnutí vazeb energie deformace vazebných úhlů energie rotace vazeb (torze) deformace planárních skupin (oop) křižné členy vdW energie elektrostatická energie vodíkové vazby
Rovnice silového pole energie pnutí vazeb energie rotace vazeb (torze)
vdW energie
energie deformace vazebných úhlů
elektrostatická energie
Parametry silového pole AMBER - parm94.dat (W. Cornell et all. 1995) definice atomových typů
vazby
úhly torze vdW poloměry
náboje
Původ parametrů Vazby a úhly: délky vazeb, velikosti vazebných úhlů - z RTG krystalografie silové konstanty - z empirických dat (vibrační frekvence)
Torze: přizpůsobení torzní potenciálové funkce tak, aby co nejpřesněji reprodukovala experimentální (nebo ab initio QM) energet. bariéry 9 možností definice torze H-C-C-H
Energie
Ethan
torze
Původ parametrů Nevazebné interakce: reálné plyny (van der Waals) sublimační tepla dipólové momenty
Termodynamické funkce: tvorná a spalná tepla disociační energie vazeb
Energie pnutí vazeb harmonická funkce (Hookův zákon)
Energie pnutí vazeb kubická korekce V(r)=a(r-ro)3 anharmonická funkce (Morse) V(r)=D[1-exp(-a(r-ro))]2-D H-H: D=108.9 kcal.mol-1 re=0.74 Å a=1.86 Å-1
Deformace vazebných úhlů
Vliv silových konstant
vazby/úhly
Energie rotace vazeb (torzní)
Vliv torzních parametrů
Dvojné vazby Deformace planárních skupin E(r)=0.5 K (φ - φ0)2 E(r)=0.5 V (1 - cos 2φ) Pseudo torzní uhel E(r)=0.5δ2
δ
Nevazebné interakce Van der Waals
elektrostatika
Nevazebné interakce 12-6 Lenard-Jones exp-6 Buckingham Repulze (Pauliho princip) Atrakce (Londonovy síly) 12-10 H-vazby
r*
Elektrostatická energie E(r) =
1 qi qj 4πε0 r ε =ε0 εr
vakuum
prostředí o relativní permitivitě εr
dielektrická konstanta (ne-)závislá na vzdálenosti
Celková energie Deformace vazebných délek V(b) Deformace vazebných uhlú V(θ) Deformace planárních skupin V(δ) Křižné členy V(b,θ,φ) Van der Waalsovské interakce V(r) Vodíkové vazby V(r) Vtotal=V(b)+V(θ)+V(δ)+V(b,θ,φ)+Vnb(r)+VHB(r)
Methan
Methan Parametry silového pole
Energie (single point)
Deformovaný methan Vazba C-H
K 340.0
Requil R 1.090 2.518
Energy 693.22
Ethan
Hyperplocha potenciální energie (PES) 1D plocha: E = f(φ)
Energie
17 16 15 -180
-60
φ
60
180
Restraints X Constraints (restraint = omezení, nátlak; constraint = donucení, přesné vymezení)
constraints = tvrdé vazné podmínky
silová konstanta K
restraints = měkké vazné podmínky
r0
r0
Aplikace metod molekulového modelování
Minimalizace energie E = f (struktury)
funkce n proměnných
Energie
= f (délky vazeb, pnutí vaz. úhlů, rotace vazeb, nevazebné interakce...) = hyperplocha
strukturní parametr
Minimum funkce z
Funkce 2 proměnných: 1. derivace v bodě xi = 0 2. derivace v bodě xi > 0
Funkce více proměnných: xi y
x
gradient v bodě xi ... vektor 1. derivací Hessian v bodě xi ... matice 2. parc. derivací
Algoritmus minimalizace Spádové metody: metoda největšího spádu (steepest descent)
směr negativního gradientu v bodě xi
iterativní proces: Energie
• startovní bod • směr pohybu? • prohledání regionu
konvergence?
• konvergenční kritéria? • nový bod
minimum
strukturní parametr
Konvergenční kritéria Energie
Rozdíl mezi 2 po sobě následujícími kroky:
rozdíl energie
rozdíl ve struktuře
strukturní parametr
Metody minimalizace derivační: metoda největšího spádu
daleko od minima
konjugované gradienty
blízko minima
Newton-Raphsonova metoda
pro malé molekuly vibrační frekvence
nederivační: simplexová metoda
Problém lokálních minim sedlový bod
Energie
Globální minimum
Lokální minima Reakční koordináta
Hledání minim na PES lokální metody = minimalizace (hledání lokálních minim) globální metody = hledání všech lokálních minim + globální minimum systematické prohledávání (grid, mapování) MD a odvozené techniky (simulované žíhání) Monte Carlo
Konformační prohledávání PES
Grid
Konformační prohledávání Molekulová dynamika
CICADA
Single coordinate driving
Konformační chování Dimethylfosfát
ε
ζ
ε
C-O-P-O
ζ
O-P-O-C
Konformační chování DMP
Konformační chování
O-P-O-C C-O-P-O
Solvent E(r) =
1 qi qj 4πε0 r ε =ε0 εr
vakuum
prostředí o relativní permitivitě εr
dielektrická konstanta (ne-)závislá na vzdálenosti
Explicitní solvent, H2O ☺ sledování specifických interakcí s molekulami H2O značná výpočetní náročnost
peptid (330 atomů) + 2400 molekul H2O
Dlouhodosahové elektrostatické interakce Systém N-částic: nejnáročnější část minimalizace = výpočet dlouhodosahových interakcí (Coulomb síly) N2 výpočtů elektrostatická energie
Hraniční podmínky Periodic Boundary Conditions (PBC)
Výpočet elektrostatických interakcí v PBC „odříznutí“ interakcí o vzdálenosti > cutoff
původní funkce
Energie
funkce při použití cutoffu
problém vyšších cutoffů
vzdálenost
Ewaldova sumace Periodický box
E(r) =
1 qi qj 4πε0 r
∞
E(r) = ∑
n=0
1
qi qj
4πε0 r + n
n ... násobek délky hrany boxu L
Ewald (1921): reálný prostor
E = f (L)
reciproký prostor
E = f (1/L)
Software • Hyperchem (Hypercube) • Chem-X (Chemical design) • AMBER (Kollman, U. California) • CHARMM, CERIUS • CHEMOFFICE (Cambridgesoft) • MM, MM4 (Allinger, Univ. Georgia) • ALCHEMY, SYBYL (Tripos) • BOBY (Zefirov, MGU Moskva, Springer) • MOLBLD (Boyd 1968, Pavelcik 1975)
Molekulová dynamika
Historie MD 1976 - 1985 „dark ages“ 1985 - 1994 free energy perturbation (FEP) ⇒ ΔG před. r. 1995: MD simulace nukleových kyselin (vysoce nabité systémy, nutno správně pracovat s dlouhodosahovými elektrostatickými interakcemi) - nestabilní (zkroucení duplexu, zlomení párů bazí).
Nedostatek výpočetní síly: Pro výpočet dlouhodosahových elektrostatických interakcí je užíván cutoff - useknutí interakcí - nestabilita. Neužívá se PME (známo od r. 1993) Pomalý vývoj silových polí - nemožnost kvalitních QM výpočtů Krátké simulace < 200ps
Historie MD 1995 - zvýšená počítačová síla (superpočítače) korektní zacházení s dlouhodosahovými nevazebnými interakcemi (elektrostatika, vdW) kvalitní QM výpočty na vyšší úrovni ⇒ silové pole 2. generace - např. AMBER parm94.dat ⇒ stabilní trajektorie
1998 - přesná interpretace struktur + ΔG MD simulace nukleových kyselin - DNA (duplex,..., kvadruplex), RNA (hairpin, ribozymy), modifikovaná DNA, RNA, protein/NA interakce
Molekulová dynamika Termodynamika
stav 1
stav 2
Termodynamika:
jaké stavy jsou možné?
Kinetika:
jak a jak rychle stavy konvertují?
stav n
Molekulová dynamika integrace Newtonových pohybových rovnic
Fi = mi ai
dE
dvi dt
dri
dri dE dri
d2ri = mi dt2
dt
Algoritmus MD Síla působící na částici = vektorový součet interakcí s ostatními částicemi F=m.a zrychlení
rychlost
pozice v čase t + Δt
Algoritmus MD Síla působící na částici = vektorový součet interakcí s ostatními částicemi
Aplikace termostatu: přeškálování rychlostí vzhledem k tepelné lázni
nová pozice v čase t + Δt
Počáteční konfigurace pro MD z experimentu přiřazení počátečních rychlostí - náhodně podle Maxwell - Boltzman distribuce (Gaussova distribuce) Generátor (pseudo)náhodných čísel: produkuje čísla ∈ < 0 , 1 > rovnoměrně přepočet na čísla odpovídající gaussovské distribuci
Časový krok v MD příliš malý - MD pokrývá jen omezenou část konformačního prostoru příliš velký - nestabilita trajektorie, vysoká E
Správný časový krok Δt : 1/10 nejkratšího vibračního pohybu C-H vazba vibruje s periodou 10fs ⇒ integrační krok Δt
= 1 fs
Constraints na vazby obsahující H (SHAKE algoritmus) ⇒ integrační krok Δt = 2 fs
Termostat napodobuje působení teploty na systém - působení makroskopických vibračních modů na dynamiku jednotlivých atomů = násobení každé rychlosti atomu faktorem, který tvoří požadovanou teplotu (každý krok nebo periodicky každých X kroků)
Simulované žíhání
Energie
zahřátí na vysokou teplotu
pomalé ochlazování
Postup MD+solvent volba vstupní struktury (např. z RTG), úprava (doplnění chybějících residuí, vodíků), optimalizace doplněných částí. přidání iontů (neutralizace systému), solvatace - přidání molekul vody (solvatační box) ekvilibrace systému: minimalizace vody a iontů krátká MD vody a iontů
Ekvilibrace systému Minimalizace a dynamika solventu a iontů constraints na molekulu peptidu
Několikastupňová minimalizace s postupným uvolňováním peptidu restraints na molekulu peptidu (25 kcal/mol/A, 20, 15, 10, 5, 0)
Zahřívání systému
0 ps
2 ps
Zahřívání 100 K → 300 K během 20 ps
Čas
20 ps
start produkční fáze MD
Teplota
Kontrola ekvilibrace
Produkční fáze MD Frekvence ukládání snímků = 500 kroků
1 krok MD = 0.002 ps = 2 fs
1 ps
0 ps
500 kroků
2 ps
500 kroků
1 snímek
1 snímek
... ... ...
MD trajektorie = soubor snímků
Čas
2 ns
Trajektorie - snímky
0 ps
200 ps
400 ps
600 ps
800 ps
1000 ps
1200 ps
1400 ps
1600 ps
1800 ps
Průměrná struktura z poslední fáze MD (např. 0.5 ns)
„Vázaná“ MD = MD s užitím restraints Restraints (vazné podmínky) ... získány z NMR
NMR
přiřazení signálů
vzdálenosti
MD + restraints
3D struktura
Problémy simulací Energie
?
časová škála: cíl:
μs - ms
realita: ps - ns
kvalita reprezentace: cíl:
komplexní systém
realita: molekulová mechanika
