[LNP] Step5_combine_coating

## 이전 글 ##

2026.11.04 - [분류 전체보기] - [LNP] Step4_create_top_hexagon

[LNP] Step4_create_top_hexagon

 

 

TS2CG로 지질 코팅을 만들고, 그 코팅을 육각형 코어(hexagon_fake_whole.pdb)와 합쳐서 코팅된 구조 coated_test.pdb를 만든다. 아래에서 더 자세히 뜯어볼 예정이다.

 

우선 다시 전체 흐름을 복기해 보겠다.

• Step 1~4: 4개 채널 생성 → 채널 추출 → 육각형 코어로 복제 → 코어 topology 작성
• VMD와 Blender로 코어 표면을 삼각형 메쉬(.tsi)로 만듦
• Step 5: 코팅을 생성하고 코어와 합침
• Step 6 이후: 코팅까지 포함한 topology 작성 → 진공 minimization → solvation과 중화

 

 

---

# 셀 라이브러리 import

import numpy as np
import MDAnalysis as md

import warnings
warnings.filterwarnings('ignore')

줄	하는 일	쓰이는 곳
numpy as np	수치 계산	np.max, np.min 로 코어의 z 범위를 구하고, np.unique 로 residue 번호의 중복을 없앰
MDAnalysis as md	구조 파일(gro/pdb)을 읽고, 원자를 선택하고, 합치고, 쓰기	md.Universe, select_atoms, md.Merge, .write 등 핵심 작업
warnings.filterwarnings('ignore')	모든 경고 메시지를 끔	MDAnalysis가 pdb/gro를 읽을 때 원소 정보나 질량을 추정하지 못하여 나오는 경고를 숨기는 용도. (Martini CG bead 이름은 실제 원소가 아니라서 경고가 많이 나옴)

 

 

---

그 다음으로는 아래와 같은 경고가 나와있다.

## Important note: Please do not proceed if you are not familiar with TS2CG.
### You can learn more about TS2CG here: https://cgmartini.nl/docs/tutorials/Martini3/TS2CG/

 

TS2CG에 대해서 익숙하지 않으면 이후 과정을 진행하지 말라는 안내문이다.

TS2CG에 대해서 공부가 필요하면 아래 사이트를 참고하면 된다.

 

https://cgmartini.nl/docs/tutorials/Martini3/TS2CG/

TS2CG 2.0 as a membrane builder – Martini Force Field Initiative

 

 

---

# Step 1: Create lipid coating for hexagon core

Step 1은 TS2CG 지질 조성 파일 input.str 을 작성하는 과정이다.

IL = 'MC3H'

build_file = open('input.str' , 'w')
build_file.write(f'[Lipids List]\n')
build_file.write(f'Domain 0\n')
build_file.write(f'{IL}  0.25  0.25 0.45\n')
build_file.write(f'DSPC  0.1  0.1 0.45\n')
build_file.write(f'CHOL  0.25  0.25 0.45\n')
build_file.write(f'End\n')
build_file.close()

 

TS2CG는 다음의 두 단계로 막을 만든다.

1. plm 이 표면 매쉬를 받아 지질을 놓을 자리를 계산한다.
2. pcg 가 그 자리에 실제 지질 분자를 배치한다.

input.str 은 pcg 에게 어떤 지질을, 어떤 비율로, 얼마나 촘촘히 놓을지 알려주는 설정 파일이다.

 

IL = 'MC3M' 은 이온화 지질의 이름을 변수로 둔다.

이는 다른 이온화 지질로 바꿔 보기 쉽게 하려는 것이며, MC3H는 DLin-MC3-DNA의 양성자화된(+1 전하) 형태이다.

산성 조건의 코어를 가정하며, 이 이름은 Martini3.LIB (분자의 bead 좌표 템플릿)에 정의되어 있어야 한다.

 

[Lipids List] ... End 사이에 지질 목록을 작성한다.

TS2CG는 표면을 여러 도메인으로 나누어 도메인마다 다른 조성을 줄 수 있다.

어느 도메인인지는 .tsi 메쉬의 꼭짓점마다 지정된다. 여기서는 도메인 0 하나만 써서 표면 전체를 같은 조성으로 덮는다.

 

지질 목록은 ( 이름    위 leaflet 비율    아래 leaflet 비율    APL )의 형식으로 작성한다.

• 비율은 같은 leaflet 안에서의 상대 비율을 의미한다. 위아래 값을 같게 둬서 양쪽 leaflet의 조성이 대칭이다. monolayer모드에서는 한쪽 값만 쓰인다.
• APL(area per lipid) : 지질 하나가 차지하는 면적으로, 이 값을 통해 같은 표면에 지질을 몇 개 놓을지 정한다. 작게 잡을수록 촘촘해진다.

 

---

# 메쉬 out_test.tsi 준비하기

이 단계는 A, B 두 부분으로 나뉜다.

 

 

A. hexagon_fake_whole.pdb 만들기

 

이전 노트북 파일에서 만든 Hexagon.pdb를 열면 분자 206개가 z 방향 주기 경계에서 둘로 잘려있다.

한쪽 끝은 z ≒ 0 에, 다른 쪽 끝은 z ≒ 103에 있다.

코어 채널은 z 방향으로 끝없이 이어진 구조를 잘라낸 것이라 이 상태로는 표면을 만들 수 없다.

 

GROMACS로 잘린 분자를 다시 이어 붙여준다.

이 부분은 노트북 파일에 설명이 나와있지 않아 유추해야 하지만, 채널 자체는 그대로 두고 분자 하나하나만 온전하게 만들어 이름이 "fake whole"일 것으로 추정해본다.

 

 

1. 작업 폴더로 이동

 

먼저 작업 파일로 이동한다.

$GMX 는 GROMACS 실행 파일 경로를 줄여 쓰려는 변수이다.

cd 작업파일위치
GMX=절대경로/gromacs-2026.3-install/bin/gmx

 

 

2. grompp로 .tpr 만들기

$GMX grompp -f MDPs/min_vac.mdp -c Hexagon.pdb -p hexagon.top -o hex_whole.tpr -maxwarn 5

 

잘린 분자를 이어 붙이려면 GROMACS가 "어느 원자끼리 한 분자이고 서로 결합되어 있는지"를 알아야 한다.

이 정보는 좌표 파일(pdb)에는 없고 topology에만 있다. .tpr은 좌표와 topology를 하나로 묶은 파일이다.

• -f MDPs/min_vac.mdp : 시뮬레이션 설정 파일로, 여기서는 시뮬레이션을 돌리지 않으니 형식상 아무 mdp나 넣는다.
• -c Hexagon.pdb : 코어 좌표 (Step 3 노트북 결과)
• -p hexagon.top : 코어 topology (Step 4 노트북 결과)
• -o hex_whole.tpr : 출력 파일명
• -maxwarn 5 : 사소한 경고로 인해 멈추지 않게 함

저 코드를 실행하면 hex_whole.tpr과 mdout.mdp (설정 기록)가 생긴다.

 

 

3. trjconv 로 분자 이어 붙이기

echo 0 | $GMX trjconv -s hex_whole.tpr -f Hexagon.pdb -pbc mol -o hexagon_fake_whole.pdb

• -pbc mol : 분자마다 중심을 상자 안에 두고, 경계 너머로 넘어간 원자를 같은 쪽으로 붙여서 분자를 온전하게 만듦
• echo 0 : 어떤 원자들을 출력할지 묻는 질문에 0 (System, 전체) 라고 자동으로 답한다.

결과는 hexagon_fake_whole.pdb 가 나온다.

파일	원자 수	가장 긴 분자
Hexagon.pdb	22995	102.8 Å (잘려서 상자 양 끝에 걸쳐 있음)
hexagon_fake_whole.pdb	22995	22.8 Å (정상)

 

 

---

 

B. out_test.tsi 만들기

 

이제 다음을 진행한다.

1. VMD에서 코어를 QuickSurf 표면으로 그려 vmdscene.stl 로 저장한다.
2. Blender에서 표면을 고른 삼각형으로 다시 짜고(voxel remesh) TS2CG 형식으로 저장한다.

 

먼저 프로그램 VMD와 Blender를 아래 링크에서 각각 다운로드 받아준다.

VMD의 경우에는 이메일 등록이 필요하다.

 

- VMD

https://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD

Software Downloads

 

- Blender

https://www.blender.org/download/

Download — Blender

 

 

---

 

다운로드 받은 vmd 파일을 압축 해제하고 인스톨 해준다.

tar -xzf vmd-2.0.1a1.bin.LINUXAMD64.tar.gz
cd vmd-2.0.1a1

sudo VMDINSTALLBINDIR=$HOME/.local/bin
sudo VMDINSTALLLIBRARYDIR=$HOME/.local/lib/vmd

./configure
cd src
make install

 

설치가 끝나면 확인하고 파일 생성을 진행한다.

which vmd
cd 실행폴더위치
vmd -dispdev text -e render_stl.tcl -args hexagon_fake_whole.pdb vmdscene.stl
ls -la vmdscene.stl

 

 

 

이제 이어서 Blender도 이어가면 된다.

sed -i 's/bpy.ops.vm.stl_import(/bpy.ops.wm.stl_import(/' remesh.py
grep -n "stl_import" remesh.py     # bpy.ops.wm.stl_import(filepath=in_file) 로 나와야 함
ls -la vmdscene.stl                # VMD 단계 결과물이 있는지 확인
./blender-5.2.3-linux-x64/blender --background --python remesh.py
ls -la out_test.tsi

 

 

이렇게 하면 out_test.tsi를 만들 수 있다.... 아 힘들어;;

 

 

 

다음으로 아래 코드를 실행해주면 된다.

 

처음에 -rescalefactor를 13으로 진행하고 코팅과 코어 사이의 틈과 겹침을 재보았다.

13은 겹침은 없지만 코팅이 코어에서 1 nm 넘게 떠있었고, 이 틈이 나중에 vaccum minimization으로는 닫히지 않거나 나중에 물을 채울 때 이 틈으로 물이 들어갈 위험이 있어서 수치를 12로 바꾸어서 돌려주었다.

 

 

#mono 2
$TS2CG/PLM -TSfile out_test.tsi -rescalefactor 12 12 12 -smooth -monolayer 1 -o mono_layer2
$TS2CG/PCG -str input.str -LLIB Martini3.LIB -defout mono2 -dts mono_layer2

#Bilayer
$TS2CG/PLM -TSfile out_test.tsi -rescalefactor 12 12 12 -smooth -o Bi-layer2 -bilayerThickness 2
$TS2CG/PCG -str input.str -LLIB Martini3.LIB -defout Bi-layer2 -dts Bi-layer2

 

나는 Bilayer 먼저 돌리고 mono layer를 돌려주었다.

이 코드는 코어를 감쌀 지질 껍질을 만드는 역할을 한다.

 

TS2CG는 "모양 틀을 주면 그 표면에 지질 분자를 빼곡히 꽂아 주는" 프로그램이다.

틀은 out_test.tsi가 되며, VMD와 Blender로 생성해준 코어 겉모양의 표면 그물망이다.

이 틀에 PLM이 표면 위에 "여기에 지질 하나" 를 의미하는 점을  촘촘히 찍고, 점마다 지질이 향할 방향도 정해준다.

 

- rescalefactor는 틀이 실제 크기보다 작게 저장되어 있어서 12배로 키워주는 것이다. 

13배로는 껍질이 너무 커서 코어와 틈이 생기기 때문에 12로 바꿔주었다.

 

-monolayer는 점을 1겹만 찍는다는 의미이고, -bilayerThickness 2 를 주고 -monolayer를 빼면, 표면 안팎으로 2 nm 간격을 두고 2겹을 찍는다.

 

 

이제 PCG는 그 자리에 실제 지질 분자를 꽂아준다.

input.str 은 어떤 지질을 어떤 비율로 놓을지 정한 파일이다. 

 

 

사실 터미널에 돌리는 순서는 중요하지 않다. 둘 다 코어 전체를 감싸는 껍질이 각각 형성되기 때문에 다음으로 할 과정은 필요한 부분만 오려 붙이는 것이다.

 

 

 

---

# Step 2: Join hexagon core with the coating

지금까지 형성된 코어는 지름이 약 20 nm 이고, 높이가 약 11 nm 인 육각기둥의 모양이다.

채널 7개가 육각형으로 모여 있고, 채널 하나는 지질로 된 가느다란 관이다.

관 안쪽에는 물과 siRNA가 들어 있으며, 지질 머리는 안쪽 물을 향하고 꼬리는 바깥을 향하는 inverse hexagonal 배열을 가진다.

 

이 채널들은 원래 z 방향으로 끝없이 이어진 구조(구조 경계)였는데, 그걸 약 11 nm 높이에서 잘라낸 것이 지금의 코어이다.

 

그럼 지금 윗면과 옆면을 생각해 보았을 때 지질을 다르게 덮어주어야 한다.

옆면은 채널 바깥쪽 벽으로 지질 꼬리가 바깥을 향해 드러나 있어서 지질 1겹만 덮으면 된다.

윗면과 아랫면은 관을 칼로 자른 단면 모양으로 물과 RNA, 지질이 뒤섞여 드러나 있어서 지질 2겹으로 두껍게 덮어야 한다.

 

 

....라고 생각했지만 아니라고 한다..!!!

(아직 이 노트북 파일들을 제대로 이해하지 못한듯...;;;;;; 다시 공부해야겠다.)

 

두 껍질은 모두 같은 틀인 out_test.tsi 를 따라 만들어졌는데, 그 틀이 코어에 딱 붙은 모양이 아니라 코어보다 한 겹 크게 잡힌 표면이므로 괜찮다고 한다. VMD QuickSurf가 원자 바깥으로 여유를 두고 표면을 그렸고, 그걸 12배로 키웠기 때문이라고 한다.

 

 

---

- 구조 파일 세 개 불러오기

m2 = md.Universe('mono2.gro')
m1 = md.Universe('Bi-layer2.gro')
c  = md.Universe('hexagon_fake_whole.pdb')

univers = [m1, m2, c]

• md.Universe(...)는 MDAnalysis에서 구조 파일 하나를 여는 기본 명령이다.
• 열면 파일 안의 모든 원자(bead) 정보가 담긴 객체가 만들어진다. 원자 이름, residue 이름과 번호, 좌표 등이 들어있다.
• .gro 는 좌표를 nm 단위로 저장하지만, MDAnalysis는 읽을 때 자동으로 Å 으로 바꾼다. 그래서 이후 계산은 모두 Å 단위다.

 

 

- 세 구조의 중심을 원점으로 옮기기

 

for u in univers:
    u.atoms.positions = u.atoms.positions - u.atoms.center_of_geometry()

 

세 파일은 각각 다른 프로그램에서 만들어져 좌표계가 다르다.

이대로 합치면 코어와 어긋나기 때문에 셋 다 중심을 (0, 0, 0)으로 옮겨서 코어가 껍질 한가운데 오게 한다.

 

 

 

- 자를 높이 기준선 두 개 정하기

c_max = np.max(c.atoms.positions[:,2]) - 5
c_min = np.min(c.atoms.positions[:,2]) + 5

• 먼저 코어 모든 원자의 z 좌표만 뽑는다.
• 코어에서 가장 높고 낮은 원자의 z를 뽑아서 각각 5 Å 씩 안쪽으로 당긴다.

5 Å 를 안쪽으로 당겨주는 이유는, 코어 끝 높이에서 바로 자르면 윗면 조각과 옆면 조각이 딱 맞닿기만 하기 때문이다. 그러면 모서리에 틈이 생기기 쉽다. 반면에 5 Å 안쪽에서 자르면 두 조각이 위아래로 조금 겹쳐 틈이 막힌다.

 

 

 

- 값 확인

c_min

 

실행 결과 np.float32(-46.102932) 가 출력되었다.

 

 

 

- Bilayer 에서 윗면과 아랫면 조각 잘라내기

l1 = ' '.join(f'{i}' for i in np.unique(m1.select_atoms(f'prop z >= {c_max}').resids))
M1_upper = m1.select_atoms(f'resid {l1}')

l2 = ' '.join(f'{i}' for i in np.unique(m1.select_atoms(f'prop z <= {c_min}').resids))
M1_lower = m1.select_atoms(f'resid {l2}')

 

Bilayer 껍질에서 윗면 조각과 아랫면 조각을 잘라내는데, 지질이 반토막 나지 않게 하려고 두 번에 걸쳐서 고른다.

 

먼저, 선을 넘은 원자를 찾는다. 그 다음에 그 원자가 속한 지질을 통째로 가져온다.

윗면에서는 높이 51.6 Å 보다 위에 걸친 지질을 가져온다. 935개다.

아랫면에서는 높이 -46.1 Å 보다 아래에 걸친 지질을 가져온다. 1010개이다.

 

 

 

- 아랫면 조각 확인

M1_lower

 

<AtomGroup with 10830 atoms> 이 출력되면 정상이다.

 

 

 

- Monolayer에서 옆면 띠 잘라내기

l = ' '.join(f'{i}' for i in np.unique(m2.select_atoms(f'prop z < {c_max} and prop z > {c_min}').resids))
sides =  m2.select_atoms(f'resid {l}')

 

Bilayer와 같은 방식인데, 대상이 monolayer이고 조건이 -46.1 < z < 51.6 으로 차이가 있다.

원자 하나라도 이 구간에 있으면 그 지질을 통째로 가져온다. 결과는 1006개 지질, 10833개 원자이다.

 

 

 

- 네 조각 합치기

combined = md.Merge(M1_upper.atoms, M1_lower.atoms, sides.atoms, c.atoms)

• md.Merge 는 여러 원자 묶음을 이어 붙여 새 Universe 하나를 만든다.
• 좌표는 각 조직이 가진 값을 그대로 쓴다. 이전에 중심을 맞춰 두었기 때문에 제자리에 들어간다.
• 원자 순서는 윗면 → 아랫면 → 옆면 → 코어 이다. 

 

 

- 파일로 저장하기

combined.atoms.write('coated_test.pdb')

 

합친 구조를 coated_test.pdb 로 저장한다. 이 파일은 다음 노트북 파일의 입력이 된다.

 

이후에 이어질 셀들은 모두 합친 구조에 지질이 몇 개 들어가는지 세는 확인용 코드이다.

 

 

 

- 코팅만 따로 합치기

coat = md.Merge(M1_upper.atoms, M1_lower.atoms, sides.atoms)

 

코어를 배고 윗면, 아랫면, 옆면만 합쳐서 코팅 지질만 세어보려는 코드다.

np.unique(coat.atoms.resids).shape

• 코팅에 있는 지질 번호의 종류 수를 센다.
• (2394,) 가 나오는데, 실제 지질은 2951개이다. 이는 monolayer, bilayer 파일 둘 다 지질 번호를 1부터 매겨서, 번호가 같은 지질이 하나로 세어진 것이다.
• 정확히 세려면 coat.residues.n_residues 를 쓰면 된다.

 

 

- 저장한 파일 다시 열기

u = md.Universe('coated_test.pdb')

 

파일이 제대로 저장됐는지 확인하는 코드이다.

 

 

 

- 지질 종류별 개수 세기

L1P = u.select_atoms('resname MC3H and name NP').n_atoms
CHOL = u.select_atoms('resname CHOL and name ROH').n_atoms
DSPC = u.select_atoms('resname DSPC and name PO4').n_atoms

 

지질 하나에 딱 하나씩만 있는 bead를 골라 개수를 센다. 그러면 bead 수가 곧 분자 수가 된다.

MC3H는 NP bead (머리의 아민), ChOL은 ROH bead (머리의 OH), DSPC는 PO4 bead (인산기)가 된다.

이는 코어 안의 지질까지 포함한 전체 개수이다.

 

이어지는 셀에서 이 값으로 비율(%)를 계산한다.

 

 

 

- 지질 총 개수

total = L1P + CHOL + DSPC

 

세 지질 개수를 더한다. 2216 + 1707 + 491 = 4414 개이다.

 

 

 

- 지질별 비율(%)

(L1P/total)*100    # MC3H
(CHOL/total)*100   # 콜레스테롤
(DSPC/total)*100   # DSPC

 

각 지질 수를 총 개수로 나누고 100을 곱해 %로 바꾼다.

 

각각 MC3H 50.2%, CHOL 38.7%, DSPC 11.1% 로 나오면 된다.

코어 지질까지 포함한 LNP 전체 조성이라서 이전에 정한 코팅 조성 (약 42:42:17)과 다르게 나오는게 정상이다.

코어에는 MC3H가 많고 DSPC가 없기 때문이다.

 

 

 

- 들어 있는 분자 종류 확인

np.unique(u.select_atoms('all').resnames)

 

구조 안의 residue 이름을 중복 없이 보여준다.

A, U, MC3H, CHOL, DSPC, W, ION 이 7개만 나오면 정상이다.

 

 

 

 

ㅠㅠㅠ~~~~ 이제 끝이다!!!!

다음주 화요일에 있는 조인트 미팅에서 이거 어떻게 돌아가는지 발표해야 하기 때문에 빨리 다음 노트북으로 넘어가야 한다....