Lab 5: Proteinfolding#
I moderne biokjemi er beregningsbaserte verktøy blitt veldig viktige. Denne datalaben gir deg erfaring med to svært ulike tilnærminger til proteinstruktur: atomær simulering (her: Monte Carlo-simulering) og KI-drevet prediksjon. Du vil se hva en Monte Carlo-simulering faktisk gjør – og hva den ikke klarer – og bruke dette som utgangspunkt for å forstå hvorfor metoder basert på kunstig intelligens, som AlphaFold, representerer et gjennombrudd. Gjennom praktiske analyser i PyMOL og AlphaFold ska du i denne laben arbeide med proteinstrukturer, legemiddelbinding, mutasjonseffekter og prediksjonsusikkerhet.
Del 1: Monte-Carlo-simuleringer#
I denne simuleringen bruker vi en Monte Carlo-metode (MC) for å utforske hvordan et protein kan folde seg. I stedet for å følge proteinets faktiske bevegelse over tid (slik molekylærdynamikk gjør), foreslår vi tilfeldige endringer i proteinets form — her ved å rotere en del av kjeden rundt et tilfeldig punkt. For hvert forslag beregner vi om den nye konformasjonen har lavere eller høyere energi enn den forrige. Hvis energien går ned, aksepterer vi alltid det nye forslaget. Hvis energien går opp, aksepterer vi det likevel med en viss sannsynlighet gitt av Metropolis-kriteriet:
ΔE er energiøkningen og T er temperaturen. Dette betyr at små forverringer aksepteres oftere enn store, og at høy temperatur gjør systemet mer «eventyrlystent». Denne mekanismen hindrer simuleringen i å sette seg fast i en dårlig konformasjon tidlig, og lar den utforske et bredt landskap av mulige strukturer. Etter hvert senkes temperaturen gradvis — dette kalles simulert annealing — slik at systemet til slutt konvergerer mot en stabil, lavenergikonformasjon.
Kjør koden nedenfor. Denne implementerer en enkel foldealgoritme basert på MC med hydrofobe residuer (H) og hydrofile/polare residuer (P). Se animasjonen og beskriv hva du observerer:
Observer foldingen. Hvordan endrer proteinets konformasjon seg over tid?
Hvordan varierer energien under simuleringen? Legg merke til både den nåværende energien og den beste energien som oppnås underveis.
Hvilken rolle spiller tilfeldige bevegelser i prosessen?
Hvilke faktorer tar modellen hensyn til?
Hva mangler for å gjøre simuleringen mer realistisk?
Endre parameteren “sequence” systematisk og kjør simuleringen. Prøv ulike sekvenser: En med flere/bare hydrofobe (H) aminosyrer, en med flere/bare polare (P) aminosyrer, og blandede sekvenser.
Diskuter hvordan aminosyresekvensen påvirker den endelige foldingen. Hvilke faktorer er viktige i denne modellen, og hvilke effekter er viktige i virkeligheten? Hvordan påvirker antall H og P den endelige energien?
Tips: Det kan være lurt å bruke plot_best_configuration-funksjonen (som plottes til slutt) for å visualisere den beste konformasjonen for hver sekvens du tester.
import random
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# ─────────────────────────────────────────────────────────────
# HP-MODELL FOR PROTEINFOLDING (Monte Carlo / simulert annealing)
#
# Energifunksjon: Gaussisk brønnpotensial mellom H–H-par.
# E_ij = –exp(–(r – 1)² / 0.5) for r < 3.0
# Minimum ved r = 1 (E = –1), avtar jevnt mot 0 med avstand.
# Energien er et kontinuerlig (ikke-heltalls) tall.
#
# Trekk: KUN pivot-rotasjoner (90°/180°/270° om én akse).
# Pivot-trekk garanterer at alle bindingslengder forblir = 1,
# slik at kjeden aldri brytes.
#
# Algoritme: simulert annealing med eksponentiell kjøleplan.
# ─────────────────────────────────────────────────────────────
### ENERGIFUNKSJON ###
def beregn_energi(posisjoner, sekvens):
"""
Beregner total energi ved hjelp av et Gaussisk brønnpotensial
mellom alle ikke-bundne H–H-par med avstand < 3.0.
Potensial: E_ij = –exp(–(r – 1)² / 0.5)
Minimum (–1) når r = 1, avtar jevnt mot 0 ved større avstand.
"""
energi = 0.0
n = len(sekvens)
cutoff = 3.0
for i in range(n):
for j in range(i + 2, n):
if sekvens[i] == 'H' and sekvens[j] == 'H':
dx = posisjoner[i][0] - posisjoner[j][0]
dy = posisjoner[i][1] - posisjoner[j][1]
dz = posisjoner[i][2] - posisjoner[j][2]
r = math.sqrt(dx*dx + dy*dy + dz*dz)
if r < cutoff:
energi += -math.exp(-((r - 1)**2) / 0.5)
return energi
### STARTKONFIGURASJON ###
def generer_startposisjoner(sekvens):
"""
Legger alle residuer på rad langs x-aksen (utstrakt kjede).
Bindingslengde mellom naboer er 1 overalt.
"""
return [(float(i), 0.0, 0.0) for i in range(len(sekvens))]
### PIVOT-TREKK ###
def roter_vektor(v, matrise):
x, y, z = v
return (
matrise[0][0]*x + matrise[0][1]*y + matrise[0][2]*z,
matrise[1][0]*x + matrise[1][1]*y + matrise[1][2]*z,
matrise[2][0]*x + matrise[2][1]*y + matrise[2][2]*z,
)
ROTASJONER = [
[[1,0,0],[0,0,-1],[0,1,0]], # x 90°
[[1,0,0],[0,-1,0],[0,0,-1]], # x 180°
[[1,0,0],[0,0,1],[0,-1,0]], # x 270°
[[0,0,1],[0,1,0],[-1,0,0]], # y 90°
[[-1,0,0],[0,1,0],[0,0,-1]], # y 180°
[[0,0,-1],[0,1,0],[1,0,0]], # y 270°
[[0,-1,0],[1,0,0],[0,0,1]], # z 90°
[[-1,0,0],[0,-1,0],[0,0,1]], # z 180°
[[0,1,0],[-1,0,0],[0,0,1]], # z 270°
]
def forsøk_pivot(posisjoner, sekvens, temperatur, grenser):
"""
Pivot-trekk:
1. Velg et tilfeldig indre residuum som pivot.
2. Roter alle residuer på den ene siden av pivoten med en
tilfeldig 90°/180°/270°-rotasjon om én akse.
3. Avvis trekket dersom det oppstår overlapp eller grensebrudd.
4. Aksepter/forkast etter Metropolis-kriteriet:
– Alltid aksepter dersom ny energi ≤ gammel energi.
– Aksepter med sannsynlighet exp(–ΔE / T) ellers.
Pivot-trekk bevarer alltid bindingslengde = 1 for alle nabopar.
"""
n = len(posisjoner)
if n < 3:
return posisjoner, False
pivot = random.randint(1, n - 2)
pivot_pos = posisjoner[pivot]
rotasjon = random.choice(ROTASJONER)
nye_pos = list(posisjoner)
for i in range(pivot + 1, n):
rel = (
posisjoner[i][0] - pivot_pos[0],
posisjoner[i][1] - pivot_pos[1],
posisjoner[i][2] - pivot_pos[2],
)
rotert = roter_vektor(rel, rotasjon)
nye_pos[i] = (
pivot_pos[0] + rotert[0],
pivot_pos[1] + rotert[1],
pivot_pos[2] + rotert[2],
)
# Avvis ved overlapp (to residuer på samme posisjon)
if len(set(nye_pos)) < n:
return posisjoner, False
# Avvis ved grensebrudd
for pos in nye_pos:
if not (grenser['x'][0] <= pos[0] <= grenser['x'][1] and
grenser['y'][0] <= pos[1] <= grenser['y'][1] and
grenser['z'][0] <= pos[2] <= grenser['z'][1]):
return posisjoner, False
gammel_E = beregn_energi(posisjoner, sekvens)
ny_E = beregn_energi(nye_pos, sekvens)
delta_E = ny_E - gammel_E
if delta_E <= 0 or random.random() < math.exp(-delta_E / temperatur):
return nye_pos, True
return posisjoner, False
### SIMULERT ANNEALING ###
def kjør_simulering(
sekvens = "HPPPHHHHPPHHHH",
antall_steg = 20000,
starttemperatur = 2.0,
slutttemperatur = 0.05,
opptaksintervall = 50,
grenser = None,
):
"""
Kjører MC-simuleringen med simulert annealing.
Temperaturen kjøles ned eksponentielt fra starttemperatur til
slutttemperatur: T(steg) = T_start · (T_slutt/T_start)^(steg/N)
Dette gir god utforskning tidlig og fin konvergens mot slutten.
Parametre:
----------
sekvens : Aminosyresekvens med 'H' og 'P'
antall_steg : Totalt antall pivot-forsøk
starttemperatur : Høy T → mange trekk aksepteres (utforskning)
slutttemperatur : Lav T → nesten bare forbedringer aksepteres
opptaksintervall : Lagre konformasjon for animasjon hvert N-te steg
"""
if grenser is None:
grenser = {
'x': (-10, len(sekvens) + 10),
'y': (-10, 10),
'z': (-10, 10),
}
posisjoner = generer_startposisjoner(sekvens)
energi = beregn_energi(posisjoner, sekvens)
beste_pos = list(posisjoner)
beste_energi = energi
pos_historikk = [list(posisjoner)]
E_historikk = [energi]
kjølefaktor = (slutttemperatur / starttemperatur) ** (1.0 / antall_steg)
T = starttemperatur
for steg in range(1, antall_steg + 1):
T *= kjølefaktor
posisjoner, _ = forsøk_pivot(posisjoner, sekvens, T, grenser)
energi = beregn_energi(posisjoner, sekvens)
if energi < beste_energi:
beste_energi = energi
beste_pos = list(posisjoner)
if steg % opptaksintervall == 0:
pos_historikk.append(list(posisjoner))
E_historikk.append(energi)
return beste_pos, beste_energi, pos_historikk, E_historikk, grenser
### ANIMASJON ###
def animer_folding(
sekvens = "HPPPHHHHPPHHHH",
antall_steg = 20000,
starttemperatur = 2.0,
slutttemperatur = 0.05,
opptaksintervall = 50,
grenser = None,
intervall_ms = 100,
):
"""
Kjører simuleringen og animerer foldingsprosessen.
Røde kuler = H (hydrofobe), blå kuler = P (polare).
Returner animasjonsobjektet. Vis det med:
HTML(anim.to_jshtml())
som siste linje i Jupyter-cellen.
"""
beste_pos, beste_energi, pos_historikk, E_historikk, grenser = kjør_simulering(
sekvens, antall_steg, starttemperatur, slutttemperatur,
opptaksintervall, grenser,
)
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection='3d')
fargekart = {'H': 'crimson', 'P': 'steelblue'}
def oppdater(bilde):
ax.clear()
pos = pos_historikk[bilde]
xs = [p[0] for p in pos]
ys = [p[1] for p in pos]
zs = [p[2] for p in pos]
farger = [fargekart.get(r, 'gray') for r in sekvens]
ax.scatter(xs, ys, zs, c=farger, s=120)
for i in range(len(pos) - 1):
ax.plot(
[pos[i][0], pos[i+1][0]],
[pos[i][1], pos[i+1][1]],
[pos[i][2], pos[i+1][2]],
color='gray', linewidth=2,
)
# Tette aksegrenser rundt kjeden så den alltid fyller plottet
margin = 2
cx = (max(xs) + min(xs)) / 2
cy = (max(ys) + min(ys)) / 2
cz = (max(zs) + min(zs)) / 2
rekkevidde = max(max(xs)-min(xs), max(ys)-min(ys), max(zs)-min(zs), 4) / 2 + margin
ax.set_xlim(cx - rekkevidde, cx + rekkevidde)
ax.set_ylim(cy - rekkevidde, cy + rekkevidde)
ax.set_zlim(cz - rekkevidde, cz + rekkevidde)
beste_hittil = min(E_historikk[:bilde + 1])
ax.set_title(
f"Steg: {bilde * opptaksintervall}\n"
f"Nåværende energi: {E_historikk[bilde]:.3f} "
f"Beste energi: {beste_hittil:.3f}",
fontsize=9,
)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
anim = FuncAnimation(
fig, oppdater, frames=len(pos_historikk),
interval=intervall_ms, blit=False, repeat=False,
)
print(f"Beste energi funnet: {beste_energi:.3f}")
plt.close()
return anim, beste_pos, beste_energi
### PLOT AV BESTE KONFIGURASJON ###
def plot_beste_konfigurasjon(beste_pos, beste_energi, sekvens):
"""
Viser den beste konformasjonen som ble funnet under simuleringen.
Kjør denne i en egen celle etter animasjonscellen.
"""
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection='3d')
fargekart = {'H': 'crimson', 'P': 'steelblue'}
xs = [p[0] for p in beste_pos]
ys = [p[1] for p in beste_pos]
zs = [p[2] for p in beste_pos]
farger = [fargekart.get(r, 'gray') for r in sekvens]
ax.scatter(xs, ys, zs, c=farger, s=100)
for i in range(len(beste_pos) - 1):
ax.plot(
[beste_pos[i][0], beste_pos[i+1][0]],
[beste_pos[i][1], beste_pos[i+1][1]],
[beste_pos[i][2], beste_pos[i+1][2]],
color='gray', linewidth=2,
)
ax.set_xlim((-10, len(sekvens) + 10))
ax.set_ylim((-10, 10))
ax.set_zlim((-10, 10))
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title(f"Beste konfigurasjon\nEnergi: {beste_energi:.3f}")
plt.show()
# ═══════════════════════════════════════════════════════════════
# KJØR SIMULERINGEN HER – endre sekvens og parametere etter ønske
# ═══════════════════════════════════════════════════════════════
sekvens = "HPPPHHHHPPHHHH"
anim, beste_pos, beste_energi = animer_folding(
sekvens = sekvens,
antall_steg = 20000,
starttemperatur = 2.0,
slutttemperatur = 0.05,
opptaksintervall = 50,
intervall_ms = 100,
)
# HTML(anim.to_jshtml()) må være siste linje i cellen for å vise animasjonen
HTML(anim.to_jshtml())
Beste energi funnet: -18.890
Animation size has reached 20994873 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.