Lesson 11 of 20

OpenMM — 사진 한 장에서 100나노초짜리 영상으로

제7장 「37도씨, 생리적 조건에서 100나노초. DFG-in 94.3%」

여기까지의 결론은 전부 사진 한 장에서 나왔습니다. AlphaFold 가 내놓은 좌표는 움직이지 않습니다.

그런데 단백질은 체온 37 °C 에서 끊임없이 떱니다. 사진에서 열려 있던 루프가, 실제로는 1 나노초 뒤에 도로 닫힐지도 모릅니다.

그걸 확인하는 도구가 분자동역학(Molecular Dynamics, MD) 이고, 이야기에서 쓴 프로그램이 OpenMM 입니다.

High school사진 한 장과 동영상의 차이

친구가 달리기하는 사진 한 장을 봤습니다. 다리가 벌어져 있네요. 이 친구는 빠른 걸까요?

알 수 없습니다. 딱 그 순간의 모습이니까요. 동영상을 봐야 압니다.

AlphaFold 가 준 건 사진 한 장입니다

	주는 것	알 수 있는 것
AlphaFold	정지 사진	어떤 모양인가
분자동역학 (OpenMM)	동영상	계속 그 모양인가, 얼마나 떠는가

단백질은 체온 37 ℃ 에서 끊임없이 떱니다. 사진에서 열려 있던 문이 1 나노초 뒤에 닫힐지도 모릅니다. 그걸 확인하려고 OpenMM 을 돌린 겁니다.

MD 가 하는 일은 딱 한 줄이다

F = m a

고등학교 물리 그대로입니다. 원자마다 받는 힘을 계산하고, 그 힘으로 가속도를 얻고, 아주 짧은 시간만큼 움직입니다. 그리고 반복합니다.

단계	무엇을 하나
1	지금 원자 위치에서 모든 원자쌍의 힘을 계산한다
2	F = ma 로 가속도를 얻는다
3	Δt 만큼 위치와 속도를 갱신한다
4	1로 돌아간다

이야기의 코드에서 Δt = 0.004 ps = 4 펨토초 입니다. 100 ns 를 채우려면 25,000,000 걸음이 필요합니다. 코드에 그대로 적혀 있습니다.

simulation.step(25_000_000) # 4fs x 25M steps = 100ns

High schoolF = ma 를 2500만 번 반복한다

분자동역학의 계산은 놀랍게 단순합니다. 중학교에서 배운 F = ma 입니다.

순서	하는 일
①	지금 위치에서 원자마다 받는 힘을 구한다
②	F = ma 로 가속도를 구한다
③	아주 짧은 시간 동안 움직인다
④	①로 돌아간다

'아주 짧은 시간' 이 얼마나 짧은가

4 펨토초입니다. 1초의 1000조분의 4.

100 나노초를 채우려면 이걸 2,500만 번 반복해야 합니다. 코드에 simulation.step(25_000_000) 이라고 그대로 적혀 있어요.

어려운 건 계산식이 아니라 반복 횟수입니다 — 그래서 슈퍼컴퓨터가 필요합니다

이야기의 코드를 한 줄씩

코드	무슨 뜻
`ForceField("amber14-all.xml", ...)`	역장(force field). 어떤 원자가 어떤 힘을 주고받는지 적어 놓은 규칙표. 실험과 양자계산으로 맞춰 놓은 근사다
`tip3pfb.xml`	물 모형. 단백질만 띄워 놓으면 실제와 완전히 다르게 움직인다. 물을 넣어야 한다
`nonbondedMethod=app.PME`	먼 거리 전기력을 다루는 방법. 없으면 전하가 있는 계에서 결과가 망가진다
`nonbondedCutoff=1.0 nm`	1 nm 보다 먼 원자끼리는 직접 계산을 생략한다. 안 그러면 계산이 안 끝난다
`LangevinMiddleIntegrator(310*kelvin, ...)`	310 K = 37 °C. 온도를 유지하면서 무작위 충돌까지 넣는 방식

310 K 라는 숫자가 이 시뮬레이션의 요점입니다. "체온에서도 그러한가" 를 묻는 것이니까요. 0 K 에서라면 아무것도 안 움직입니다.

이야기의 3.2 ± 0.4 Å 도 같은 방식입니다 — 프레임마다 잰 값의 평균과 표준편차입니다.

Python Runs in your browser

# MD 궤적이 어떻게 생겼는지, 그리고 '유지율 94.3%' 가 무엇을 센 값인지 봅니다.
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(943)
n_frames = 2000                      # 100 ns 를 50 ps 마다 저장했다고 치면 2000 장
t = np.linspace(0, 100, n_frames)    # ns

# 변이형: DFG-in(켜짐) 우물에 앉아 가끔만 넘본다 (Ornstein-Uhlenbeck 흉내)
def traj(center, pull, noise, seed):
    r = np.random.default_rng(seed); x = np.empty(n_frames); x[0] = center
    for i in range(1, n_frames):
        x[i] = x[i-1] + pull*(center - x[i-1]) + r.normal(0, noise)
    return x

ang_mut = traj(67.8,  0.03, 3.2, 1)
ang_wt  = traj(142.3, 0.05, 2.4, 2)

DFG_IN = 105.0                       # 이 각도보다 작으면 DFG-in(켜짐) 이라고 부르기로 한다
hold_mut = (ang_mut < DFG_IN).mean()
hold_wt  = (ang_wt  < DFG_IN).mean()

fig, ax = plt.subplots(figsize=(9, 3))
ax.plot(t, ang_mut, lw=.7, color="crimson", label="D816V")
ax.plot(t, ang_wt,  lw=.7, color="steelblue", label="wild type")
ax.axhline(DFG_IN, ls="--", c="gray"); ax.set_xlabel("time (ns)"); ax.set_ylabel("DFG angle (deg)")
ax.legend(fontsize=8); plt.tight_layout(); plt.show()

print(f"D816V   DFG-in 유지율 : {hold_mut*100:5.1f}%   (이야기의 값 94.3%)")
print(f"WT      DFG-in 유지율 : {hold_wt*100:5.1f}%")
print(f"D816V   루프 각도 요동 : {ang_mut.std():.1f} deg")
print()
print("'94.3%' 는 2000장의 스냅샷 중 몇 장이 문턱 아래였는지를 센 것입니다.")
print("그러니 문턱을 어디에 긋느냐가 그 숫자를 만듭니다.")

High school위 코드가 방금 무슨 일을 했나요

100 나노초 동안 DFG 각도가 어떻게 흔들리는지 그래프로 그렸고, '유지율' 이라는 숫자를 계산했습니다.

유지율 94.3% 는 무엇을 센 숫자인가

궤적을 2,000장의 스냅샷으로 저장했습니다. 그중 문턱(105°) 아래인 장면이 몇 장인지 센 겁니다.

세는 것	값
전체 스냅샷	2,000장
DFG-in 이었던 장면	약 1,886장
비율	94.3%

그래서 뭘 조심해야 하나

문턱을 105° 대신 90° 로 그으면 이 숫자가 달라집니다. 94.3% 는 자연이 정한 값이 아니라 사람이 문턱을 그어서 만든 값이에요. 논문에서 이런 비율을 보면 "문턱을 어디에 뒀나" 를 꼭 찾아보세요.

MD 를 읽을 때 반드시 물어야 할 넷

얼마나 오래 돌렸나. 100 ns 는 루프 움직임에는 충분하지만 큰 구조 변화에는 턱없이 짧습니다. "100 ns 동안 안 닫혔다" 는 "영원히 안 닫힌다" 가 아닙니다.
몇 번 돌렸나. 한 번만 돌리면 우연히 그렇게 나왔을 수 있습니다. 보통 시작 속도를 바꿔 서너 번(replicate) 돌립니다. 이야기에는 그 언급이 없습니다.
평형화(equilibration)를 얼마나 버렸나. 처음 몇 나노초는 구조가 자리를 잡는 시간이라 세면 안 됩니다.
문턱을 어디에 그었나. 위 코드에서 봤듯이 '유지율' 은 문턱이 정합니다. 105° 대신 90° 로 그으면 94.3% 가 다른 숫자가 됩니다.

이 넷은 /dock 강의 13~18강에서 훨씬 자세히 다룹니다.

그래서 이 시뮬레이션이 무엇을 확인했나

활성화 루프 RMSD 변동: 3.2 +- 0.4 A   (WT 대비 지속적 차이 유지)
DFG-in 상태 유지율:    94.3%
D816V 펩타이드 표면 노출: 평균 RSA 0.79 +- 0.08

세 줄이 각각 앞의 세 강을 다시 확인합니다.

줄	확인한 것	원래 어디서 나왔나
RMSD 3.2 ± 0.4 Å	루프 차이가 정지 구조의 착시가 아니었다	8강 (3.87 Å)
DFG-in 94.3%	켜진 자세가 계속 유지된다	9강 (67.8°)
RSA 0.79 ± 0.08	표면 노출이 흔들려도 유지된다	10강 (0.83)

여기서 ± 기호가 새로 생긴 것에 주목하세요. 7~10강의 값은 숫자 하나였습니다. MD 를 거치자 숫자와 폭이 되었습니다. 그게 정지 사진과 영상의 차이입니다.

High school± 기호가 새로 생긴 것에 주목하세요

7~10강의 결과와 11강의 결과를 나란히 놓아 보세요.

7~10강 (사진)	11강 (동영상)
RMSD 3.87 Å	RMSD 3.2 ± 0.4 Å
DFG 각도 67.8°	DFG-in 94.3% 유지
RSA 0.83	RSA 0.79 ± 0.08

± 가 붙었다는 게 무슨 뜻인가

"이 값 주변에서 이만큼 흔들린다" 는 뜻입니다. 숫자 하나였던 것이 숫자와 폭이 됐어요.

키를 잴 때 "170 cm" 라고 하는 것과 "아침엔 171, 저녁엔 169, 평균 170 ± 1" 이라고 하는 것의 차이입니다. 뒤가 훨씬 정직합니다.

과학에서 좋은 결과는 대개 ± 를 달고 나옵니다

🤔 Check yourself

100 ns 시뮬레이션에서 DFG-in 이 94.3% 유지되었습니다. 이 숫자를 만든 것 중 하나가 아닌 것은?

🤔 Check yourself

MD 시뮬레이션 코드에서 constraints=app.HBonds 를 켜는 실용적 이유는?