시뮬레이션으로 해석하는 나노 세계의 역학

박여은․서지호․송채정․이수현

전북과학고

분자 시뮬레이션이 적용된 실제 연구 사례를 분석하고, 향후 응용 분야를 고찰해보자.

분자동역학 시뮬레이션(MD)는 분자의 위치와 속도를 시간에 따라 계산하며 분자의 움직임을 계산적으로 추적하는 방법이다. MD는 Force Field를 통해 얻은 힘을 입력받아 운동을 계산한다. 이때 Force Field는 크게 고전적인 Force Field와 반응성 Force Field로 나눌 수 있다.

고전적인 Force Field(Classical FF)는 비반응성 환경에서 물리적 상호작용을 다루는 데에 집중되어 있다. 이는 고전적인 Force Field가 Harmonic potential만을 따르기 때문에 나타나는 특징이다. 이러한 특징을 반영하여 고전적인 Force Field는 비반응성 환경에서 특정 물질들의 상호작용을 해석하는 연구에 주로 사용된다. 고전적 Force Field가 적용되는 분야는 대표적으로 단백질 접힘, 리간드 결합, 지질 이중층, DNA 염기 서열 구조 등이 있다. 본 R&E 보고서 6쪽에 명시된 실제 연구 사례가 이와 같은 Force Field가 사용된 경우이다. 이는 OPLS-AA 기반 사례인데, OPLS-AA(Optimized Potentials for Liqid Simulations – All Atom)는 MD에서 자주 쓰이는 Classical Force Field 중 하나이다. 유기 용매, 생체 분자, 고분자 등의 성질과 구조를 예측하기 위한 모델로, 주로 액체 상태의 분자를 설명하는 Force Field이다. 연구에서는 이를 이용하여 그래핀 나노기공을 통한 기체 분리 메커니즘 분석을 진행하였다. 그래핀은 2D 나노소재로, 본 연구에서 사용한 그래핀 모델은 1nm보다 작은 크기인 sub-nanometer이다. 이 그래핀의 탄소 고리 사이를 기공이라 하였을 때, 이산화탄소와 메테인이 그래핀 표면에서 어떻게 흡착되고 기공을 통과하며 탈착되는지에 대해 알아보기 위해서는 먼저 이산화탄소와 메테인의 성질에 대해 알아볼 필요가 있다. 이산화탄소와 메테인은 모드 비극성이고, 분자의 크기는 메테인이 더 크다. 두 기체 모두 Van der Waals 힘을 이용하여 그래핀과 상호작용이 가능한데, 이는 기체 분자의 구조 변형 없이 유도 쌍극자만을 이용하여 결합하거나 상호작용 하는 것이다. 메테인의 경우 분산력만 작용하여 흡착 에너지도 작고 방향성도 없다고 볼 수 있다. 이산화탄소의 경우, 만약 그래핀에 질소 도핑을 한다면 그래핀 기공과의 흡착력이 높아질 것이다. 이러한 고전적인 Force Field가 향후 응용될 수 있는 분야에는 기체와 고체의 상호작용-반례가 존재할 수 있으나-이나 수용액 속에서 구조만 유지하며 반응하지 않는 생체 고분자 등의 연구가 있을 것이다.

그림 . Graphene 기공을 통과하는 Gas molecule

반응성 Force Field(ReaxFF)는 결합 생상과 해리를 동적으로 반영할 수 있어 화학 반응을 나타내기 적합하다. ReaxFF는 결합 차수(bond order) 개념이 도입되어 원자 간 결합이 동적으로 형성되거나 끊어지는 반응을 처리할 수 있다. 본 R&E 보고서 6~7쪽에 언급되는 사례가 이러한 Force Fireld가 사용된 예시이다. 에탄올의 알루미늄 나노입자 표면 반응 시뮬레이션에서는 반응성 요소들이 중심을 이루도록 하여 고전적 Force Field를 사용할 수 없다. 고전적인 FF를 사용할 경우 에탄올 분자 내부의 결합과 알루미늄 분자 내부의 결합은 끊어지지 않고, 화합물 간 새로운 결합 또한 생성할 수 없다. 에탄올은 수소 결합, hydroxyl group이 존재하고 산화-환원 반응이 가능한 특징이 있고, 약한 Van der Waals, 수소 결합, 표면 극성 상호작용을 통해 초기에는 흡착된 상태일 것이다. 이후 에탄올의 수산기가 알루미늄과 결합을 형성하면, 알루미늄 표면에서 에탄올의 탈수소화가 발생한다. 생성된 수소는 알루미늄과 반응하거나 주변으로 방출된다. 이를 위해 ReaxFF는 실시간으로 결합 차수를 계산하여 원자 간 전자 밀도 분포를 조정한다. 또한 특정 결합이 끊어지거나 생성되는 조건을 자동으로 추적하여 시뮬레이션 한다. 이러한 것과 같이 반응성을 연구하는 곳에 사용될 수 있고, 탈수소 반응 등을 촉매하는 것들에 대한 에너지 방출과 보관과 관련된 연구를 진행할 수 있다. 또한 화학적 촉매의 반응을 해석하고 경로를 파악할 수 있을 것이다.

그림 . 에탄올, 알루미늄 분자의 반응

위와 같은 분석을 통해 알 수 있는 것은, MD가 특정 학문 분야에 치우치지 않고, 생물, 화학, 에너지, 입자 물리 등 다양한 부분에서 시뮬레이션을 통한 해석에 사용될 수 있다는 점이다. 또한 인간이 실험적으로, 자연적으로 관측하기 힘든 현상이나 인공지능을 통한 예측 또한 시간이 너무 오래 걸리거나 불가능 한 경우, MD를 유용하게 사용할 수 있을 것으로 추정된다. 예를 들면, 분자 단위의 실험이나 전자의 이동을 관측해야 하는 경우, 진핵생물의 체내에서 관측되는 많은 신호 전달 경로 등이 이에 해당한다.

평소 단백질-단백질 상호작용에서 진행되는 연쇄적인 신호 전달 경로에 대해 관심이 있다. 세포는 신호를 전달할 때 외부에서 신호 분자를 통해 신호룰 수용한 후 신호를 전달하는데, 이때 중계 분자를 통해서 신호를 전달하기 때문에 신호를 받은 단백질(효소) 분자들이 변화하는 것을 필요로 한다. 따라서 MD를 사용한다면 초기 상호작용은 시뮬레이션 할 수 있을 것이다. 다만, 세포 수준에서의 신호 전달 경로를 시뮬레이션 한다면 가능한 상태 조합이 기하급수적으로 증가하게 된다. 따라서 시뮬레이션으로 이를 구현하기엔 다소 무리가 있어 보이기도 한다.

분자 시뮬레이션으로 계산 가능한 물리량과 실험 결과와의 차이를 이해해보자

시뮬레이션 결과를 해석하기 위해 알아야 하는 이론적 배경과 모델은 아래와 같다.

Velocity Verlet 알고리즘

: 현재 시점의 위치와 속도, 가속도를 바탕으로 다음 시점의 값을 계산하는 방식이다. 시간에 따른 원자의 위치와 속도 계산에 사용된다.

Force Field(힘장)

: 분자 내 또는 분자 간 작용하는 다양한 상호작용을 수학적으로 모델링한 일련의 함수들의 집합이다. Force Field의 종류로 고전적 Force Field, ReaxFF가 있다.

PES(Potential Energy Surface)

: 분자 시스템의 원자 위치에 따른 포텐셜 에너지의 변화를 시각화한 함수이다. 분자 구조가 변화할 때 어떤 경로로 움직이며 어떤 에너지 장벽을 넘게 되는지 예측할 수 있다.

결합 차수(bond order)

: 원자 사이의 결합 강도와 존재 여부를 나타낸 수치이다. ReaxFF에서 원자 간 결합의 세기가 연속적으로 모델링할 수 있도록 한다.

앙상블

: 실제 가능한 수많은 물리적 성질의 상태들의 통계적 평균을 개념화한 것이다. 앙상블을 사용하여 실험 조건을 모사할 수 있다. 앙상블의 종류는 아래와 같다.

- NVE: 입자수(N), 부피(V), 에너지(E)를 고정하고 온도와 압력을 변화시킨다. 외부와 에너지 교환이 없어 고립계를 모사한다.

- NVT: 입자수(N), 부피(V), 온도(T)를 고정하고 에너지와 압력을 변화시킨다. 실온 조건의 실험을 모사할 수 있다.

- NPT: 입자수(N), 압력(P), 온도(T)를 고정하고 에너지와 부피를 변화시킨다. 실제 실험 조건과 가장 유사하다.

6. ergodicity

: 하나의 시스템을 충분히 오랜 시간 시뮬레이션 하면 그 시간 평균이 전체 상태 공간에 대한 평균과 같아진다는 성질이다. 시간 평균과 앙상블 평균이 같아지는 것으로 볼 수 있다.

이러한 것들을 이용하여 계산된 시뮬레이션 값을 실험값과 비교하기 위해서는 물리량이 통계적 평균값으로 정의된다는 것을 명확히 하고 가야한다. 실험값은 수많은 입자와 상호작용한 평균적인 결과를 의미하고 MD 시뮬레이션 결과는 충분한 시간 동안의 평균을 의미한다. 이는 ergodicity에 의한 결과이다. 또한 실험 조건과 맞는 앙상블을 골라서 시뮬레이션 해야 더 정확하고 괴리가 적은 결과를 도출할 수 있다. 시뮬레이션이 충분한 시간, 샘플링을 제공하지 못하면 에너지 장벽이 높거나 자유도가 높은 시스템에서는 전체 상태 공간을 제대로 탐색하지 못해 실험값과 시뮬레이션 값이 통계적으로 수렴하지 않아서 차이가 발생할 수 있다. 이러한 것들을 고려하여 비교하면 보다 정확한 비교가 가능할 것이다.

시뮬레이션으로 계산된 값을 측정한 값과 어떻게 비교하는지에 대해서는 아래의 표와 같은 예시를 들 수 있다.

물리량	시뮬레이션	실험
에너지	평균 퍼텐셜 에너지	열용량 등 실험값 등
투과 속도	시간에 따른 분자 수 변화율	기체 투과 실험값 등
반응 경로	bond order 변화, PES 분석	활성화 에너지

이렇게 얻은 값들은 온도, 압력, 앙상블 등의 조건을 맞추고 단위도 통일하면 비교가 가능하다. 또한 값들의 평균값, 표준편차를 계산하고 시각화 한다면 한 눈에 비교할 수 있다.

시스템을 직접 설계하고, 적절한 Force Field를 선택하여 시뮬레이션으로 물리적 특성을 예측해보자

Ubiquitin은 재활용되는 작은 단백질로, 다른 단백질과 공유결합하여 다른 단백질의 기능이나 활성 상태를 변화시키는 특성이 있다. 이를 유비퀴틴화라고 하고, 이때 기질은 세포 내에서 매우 다양하다. 세포에는 유비퀴틴화를 통한 활동 기작이 거의 대부분을 차지한다. 단백질 중 Lys(리신)을 통해 생성된 유비퀴틴 사슬을 달고 있는 단백질은 Proteasome에 의해 인식되어 분해된다. 여기서 proteasome은 세포 내에서 불필요하거나 손상된 단백질을 분해하는 항상성 유지와 관련된 분자이다. Ubiquitin은 단백질이기 때문에 아미노산으로 이루어져 있어 그 구조가 명확하게 밝혀져 있다. 따라서 MD 시뮬레이션이 가능할 것이다.

예측하고자 하는 시스템은 ubiquitin 주위의 물 분자 정렬이다. 단백질 표면에서 물 분자가 어떻게 정렬되고 분포하는지 분석하여 단백질과 용매 사이 상호작용의 원리를 이해하고자 한다. 이때 작용하는 요인으로는 정전기적 상호작용, 수소 결합, 소수성 효과 등 다양한 요인이 있을 것이다. 사용할 물 모델은 TIP3P 물 모델로, MD에서 물을 모델링하기 위해 사용되는 간단한 모델이다. 수소 원자 2개와 산소 원자 1개로 구성되어 있다. 앙상블은 생리조건에 따른 실험을 진행하므로 NPT가 적합해 보인다. 온도는 300K에 1 atm으로 설정한다. 이를 분석하는 방법으로는 단백질 표면의 산소, 질소 등의 원자와 물 분자의 산소 사이 거리 분포를 계산하는 것, 단백질과 물 사이 생성된 수소 결합이 개수, ubiguitin의 친수성, 소수성 잔기 주변의 물 밀도를 비교하는 것이 있다. 그리고 물 분자의 위치를 시각화 한다. 시스템 특성상, 펩타이드가 구조를 형성하여 기능을 갖는 단백질이 등장하고, 소분자 등과 같은 생체 분자가 용액 속에 노출된 것에 초점을 두고 있으므로 적합한 Force Field는 CHARMM이다.

그림 . ubiguitin의 작용 과정

그림 . TIP3P 물 분자 모델

만약 시뮬레이션 결과와 실험값이 차이를 보이는 경우, 그 원인으로는 크게 3개로 추정할 수 있다. 첫 번째로, 고전적인 Force Field는 실제 분자의 양자역학적 전자 구조를 단순화 시키는데, 이때 수소 결합 세기와 정전기적 인력 계산이 완벽하게 되지 않아 오차가 생겼을 수 있다. 두 번째로, ergodicity를 명확히 하지 않았을 수 있다. 오랜 시간 동안 시뮬레이션 하여 하나의 값으로 수렴시켜야 하는데 단백질의 수화 구조의 평형 완전하지 않아 완벽히 수렴되지 않았을 수 있고, 이 때문에 시간 평균이 실험에서 측정된 장기적인 평균과 다를 수 있다. 세 번째로, 실험 조건과 차이가 존재할 수 있다. 시뮬레이션과 실험이 온도, 압력, 단백질의 농도 등이 완전히 같지 않은 경우 오차가 발생할 수 있다. 이러한 요소가 일치하지 않는다면 단백질에 물 분자가 정렬되는 양상 또한 달라지게 되어 오차가 생길 수 있다.

이렇게 차이가 보일 경우 이를 보완할 수 있는 방법에는 먼저 Force Field를 다시 선택하는 것이 있다. 물 모델을 더 정밀한 전자 3개의 TIP4P 물 모델로 바꾸거나 단백질의 Force Field를 바꿀 수 있다. 또한 시뮬레이션 시간을 증가시켜 완벽하게 ergodicity를 충족시키는 방법도 있다. 오차를 명확히 하기 위해 시간 평균값과 표준편차 그래프를 그려 확인하고 순간값과 평균 구조를 비교하는 방법도 있다. 중요한 것은 시뮬레이션과 실험값이 정확하게 일치할 수 없다는 사실을 염두에 둬야한다는 것이다.