실크 펩타이드 모델은 bombyx mori silkworm 실크 서열 기반으로 구성하였다. Bombyx mori silkworm 실크는 비결정질(TGSSGFGPYVANGGYSGYEYAWSSESDFGT)구조와 결정질(GAGAGSGAGAGYGAGAGSGAAS)구조로 이루어져 있다. 특히, 결정질은 베타 시트 및 턴(turn)과 코일(coil)형태의 2차구조로 이루어져 있고, 특히 베타 시트는 실크의 구조 및 물성에 중요한 영향을 미치는 것으로 알려져 있으며 베타 시트의 비율에 따라 기계적 물성이 달라지는 것으로 보고된다⁽¹⁴⁾. 그리고 베타 시트를 구성하는 실크 가닥(strand)은 서로 역 평형(anti-parallel)로 구성되며 크기는 21 nm × 6 nm × 2 nm인 것으로 보고된다⁽¹³⁾. 이 연구에서는 베타 시트 결정질을 UCSF Chimera를 통해 실크 펩타이드 모델을 모사하였다.

모든 시뮬레이션에 앞서 에너지 최소화를 진행하였다. 이때 조건은 허용 오차 값이 1.0 kJ/mol이고 단계가 0.002 nm인 최대 1 000 000단계가 사용되었다. 전체 분자동역학 시뮬레이션은 Amber99SB-ildn force-field⁽¹⁵⁾와 함께 GROMACS 5.0.7을 사용하여 수행되었으며, xyz 방향의 주기적 경계 조건(PBC)에서 0.002 ps의 시간 간격으로 leap-frog 적분 알고리즘을 적용하였다. PBC box 유형은 계산 효율성을 최대화하기 위해 직사각형 모양으로 설정하였고, 시뮬레이션은 200 ns 동안 지속되었으며 열적 커플링(thermal coupling)의 온도 및 시간은 300 K 및 0.1 ps로 설정하였다. 압력 조건은 Parrinello-Rahman barostat으로 설정하였으며, 기준 압력은 1 bar 및 시간 상수는 2ps로 설정하였다^(16,17).

실크의 기계적 물성을 평가하기 위해 ENM을 적용하였다. ENM이론은 원자가 가지고 있는 고유의 위치에너지에 근거하고, 원자간 위치 에너지가 평형하다고 가정하며, 가정된 모델은 매우 안정적이라고 가정한다. 이 때의 ENM의 전체 위치 에너지는 다음과 같이 표현할 수 있다.

여기서, N은 C_α 원자의 총 개수, h는 헤비사이드(heaviside) 함수를 나타내며 ΔR_i 및 ΔR_j 각각 위치 R_i 및 R_j에서 i^th번째 및 j^th번째 측쇄의 변위이다. 기호 r_C는 원자 간 결합 여부를 판단하는 차단 거리(cut-off distance)로, 1.2 nm의 값일 때 그에 해당하는 힘 상수 γ는 7.113 kJ/mol이다. 추가적으로, 실크 기반 나노 소재 내 원자 간 발생하는 힘 장의 종류를 세분화하여 원자 간 비공유 결합(수소 결합)일 경우 10배⁽¹⁸⁾와 공유 결합일 경우 100배⁽¹⁹⁾의 힘 상수 값을 부여하였다. 공유결합은 실크 단백질 서열 상 i번째 원자에 대해 i-1^th와 i+1^th 원자 사이에 적용하였고, 수소결합의 경우 분자동역학 시뮬레이션 궤적 클러스터 중 대표 클러스터 평균 구조의 수소 결합 형태를 기준으로 선정하였다.

ENM의 위치 에너지는 행렬 표기법으로 식 (1)과 같이 표현할 수 있다.

여기서, ∆R은 ENM의 3N 차원 벡터이고 H는 3N × 3N 크기의 Hessian 행렬로 ENM의 강성을 표현한다. ENM에서 모드 해석(normal model analysis, NMA)을 수행하면 6가지 강체 모드를 제외한 저주파 범위 내에서 비강체 운동과 고유 진동수를 얻을 수 있다. 이러한 모드 내에서 굽힘(bending), 인장(axial), 비틀림(torsional)과 변형(deformation) 모드는 등가 연속체 탄성 구조의 인장 및 비틀림 탄성 계수를 계산하기 위해 선택되었다.

연속체 역학에서 오일러-베르누이 빔 모델(euler-bernoulli beam model)에 의해서 굽힘, 비틀림 및 축강성에 대한 운동방정식을 정리하면, 연속체 탄성 구조의 인장 및 비틀림 탄성 계수는 고유 진동수를 사용하여 다음과 같이 표현할 수 있다^(20~21).

여기서 E_B, G_B, 그리고 Y_T는 각각 실크 단백질에 대한 굽힘 강성, 비틀림 전단 계수 및 축 탄성 계수이고, ω_B, ω_A, 그리고 ω_T는 각각 굽힘, 비틀림, 축 탄성 변형에 대한 해당 고유 진동수이다. 기호 ρ는 구조의 질량 밀도를 나타내고, L은 특성 길이이며, A는 단면적, β_n는 모드 형상에 대한 종속 상수이다. 굽힘, 비틀림 및 축강성 모드에 해당하는 고유 진동수와 각각의 굽힘, 비틀림 및 축 강성 계수와 같은 기계적 특성은 변형 모드에 대한 고유 진동수와 관련하여 식 (3~5)와 같이 나타낼 수 있다. 특히 굽힘 모드에 대해서는 실크의 방향인 X축을 제외한 Y축 및 Z축에서의 굽힘을 고려하여 각각의 단면 계수를 고려하여 해당 모드 및 영률를 계산하였다. 여기서 π는 원의 둘레에 대한 지름의 비율이다.

대표적인 구조를 추출하기 위해 GROMACS의 gmx_cluster plugin을 사용하여 Jarvis-Patrick 클러스터 분석을 수행하였고, 가장 큰 클러스터를 대표 모델로 사용하였다. 클러스터는 분자동역학 시뮬레이션에 대한 결과 프레임이 도출되는데, 이때 프레임에 따라 실크 단백질의 구조적 유사성을 기준으로 같은 프레임으로 간주한 것을 의미한다. 데이터는 RMSD(root-mean-square deviation), 2차 구조 분포 분석을 통해 베타 시트의 구조를 추출하고 이를 이용하여 구조 및 안정성을 분석하였다. 이때 2차 구조 분포는 VMD(visual moelcular dynamcis) 1.9.3. 타임라인 콘솔을 사용하여 궤적 구조에서 계산하였다.

Bombyx mori 실크의 결정질에 대해서 소수성 상호작용 및 반데르발스(Van der Waals) 상호작용에 대한 에너지를 계산하기 위해 GROMACS의 gmx_mmpbsa plugin을 사용하여 에너지 계산을 수행하였다⁽²²⁾. 분석에 사용된 모델을 2.5절의 대표 클러스터 모델을 선택하였고, 이에 대해 분자 기계적(molecular mechanics, MM)에너지, polar/non-polar 에너지, 그리고 반데르발스 및 정전기적(electrostatics) 에너지를 계산하였다. 위의 모든 에너지의 총합은 결합(binding) 에너지이며, 결합 에너지가 낮아질수록 더 강하게 결합한다는 것을 의미한다. 분자 기계적 에너지는 분자 모델에 대한 본드(bond), 각도(angle), 이면체(dihedral)에 대한 결합 에너지를 의미하며, polar 및 non-polar에너지는 Poisson-Boltzmann 방정식에 의거하여 물 분자와의 상호작용을 고려하여 에너지를 계산한다⁽²³⁾.

실크 단백질의 결정질에 대하여 분자 동역학 기반 평형 시뮬레이션을 수행하였다. 결정질 모델 중에서 베타 시트는 실크 단백질 기반 소재의 기계적 물성에 크게 영향을 미치는 요인으로 알려져 있다⁽²⁷⁾. 베타 시트는 (GAGAGS)라는 실크 아미노산의 조합과 (GAGAGY)조합으로 이루어져 있다. 구체적으로, 베타 시트를 이루는 아미노산인 글리신(G), 알라닌(A), 세린(S), 그리고 티로신(Y)의 아민 그룹과 카르복실 그룹 사이의 서로 수소 결합을 하여 강하게 결합되고, 베타 시트를 구성하는 아미노산은 아미노산끼리 강하게 결합하는 소수성 작용 및 반데르발스 작용으로 강하게 결합되어 있는 것으로 알려져 있다^(24,25). 이와 반대로, (GAAS) 아미노산 조합로 되어있는 부분은 친수성 작용이 강하기 때문에 턴 및 코일 형태의 2차 구조를 형성하고 접힘 현상이 있는 것을 발견할 수 있다. 이는 이전 연구에서 보고된 경향과 일치한다⁽¹³⁾. (Fig. 1 참조) 이에 대해 정량적인 결과 분석을 통해서도 알 수 있는데, Fig. 2(a)는 분자 동역학 시뮬레이션 동안의 구조의 유사성을 RMSD 값으로 분석한 결과이다. 그 결과, 실크 단백질의 RMSD 값은 매우 안정되어 있는 것을 확인할 수 있다. 또한, Fig. 2(b)는 실크 단백질에서 결정질의 2차 구조의 비율을 분석한 결과이다. 이를 보면 실크 단백질의 기계적 물성과 매우 관련이 높은 베타 시트 비율이 약 60 %로 매우 높은 것을 확인할 수 있다. 또한 (GAAS) 아미노산 조합로 인한 코일 및 턴의 2차 구조도 형성되어 있는 것을 확인할 수 있다. 이러한 이유는 (GAGAGS) 아미노산 조합 및 (GAGAGY) 아미노산 조합이 수소 결합을 통해 연결되어 있기 때문이다. 이 때의 수소 결합의 수는 평균적으로 1082.2개인 것으로 분석되었다(Fig. 2(c) 참조). 따라서, 실크 단백질의 베타 시트는 (GAGAGS) 아미노산 조합 및 (GAGAGY) 아미노산 조합의 수소 결합에 의해서 강하게 결합하여 형성될 수 있음을 알 수 있으며, 수소 결합의 수가 많을수록 실크 단백질의 기계적 물성에 영향을 미치는 중요한 요소인 것을 확인할 수 있다.

Fig. 1 
Equilibration MD simulation of silk peptide beta-sheet

Fig. 2 
MD results of silk peptide beta-sheet

실크 단백질 기반 나노 소재의 성능을 고찰하기 위해서 탄성 네트워크 모델 및 고유 모드 해석을 이용하여 기계적 물성을 분석하였다. 이를 기반으로 실크 기반 나노 소재의 기계적 특성을 측정하였다. 여기서 강체 운동에 대한 진동 모드는 제외하고, 최저 주파수(7-12 th)는 ~ 0.1 THz이고 해당 진동 모드는 실크 기반 나노 소재의 진동 거동과 유사한 굽힘 모드이다. 단면적의 이방성으로 인해 두 가지 기본적인 굽힘 모드를 각각 soft. 굽힘 모드와 stiff. 굽힘 모드를 고려하였습니다. 각 고유 진동수에 따른 모드는 Fig. 3(a) 및 Table 1으로 정리하였다. 이에 대해 각 모드에 대해 영률은 다음과 같다. Soft. 굽힘(16.419 GPa), stiff. 굽힘(20.904 GPa), 비틀림(3.461 GPa), 그리고 인장(14.376 GPa)이다. 특히, soft. 굽힘 및 stiff. 굽힘에 대한 탄성계수가 비틀림 및 인장의 탄성계수가 높은 것을 알 수 있는데, 이는 실크 기반 나노 소재에 대한 우세한 변형이 굽힘 변형임을 의미한다. 기존 연구에 따르면 누에 실크에 대한 영률은 약 5 GPa ~ 25 GPa로 알려져 있는데⁽²⁶⁾, 이러한 경향과 매우 잘 맞는 것을 확인할 수 있다. 이렇게 우수한 기계적 물성의 원인은 실크 단백질의 결정질 중에서 베타 시트 형성 비율과 수소결합이 높아짐에 따라 더 강화되는 것을 알 수 있다⁽²⁷⁾. 이에 대한 이유는 베타 시트 형성의 중요한 원인은 수소 결합, 소수성 효과, 반데르발스 결합으로 알려져 있는데 베타 시트 구조는 이를 통해서 구조적으로 강하게 결합한 상태로 존재하기 때문이다⁽¹³⁾. 소수성 효과 및 반데르발스 결합을 추가적으로 분석하였는데, Fig. 3(c)에서 반데르발스 에너지는 약 -11.412 kJ/mol, 소수성 효과와 연관되는 분자 기계적 에너지(-13.512 kJ/mol) 및 Polar 에너지(-3.662 kJ/mol)인 것을 확인할 수 있다. 이러한 정량적 결과는 베타 시트가 반데르발스 결합 및 소수성 효과로 결합되어 있는 것을 시사한다. 최종적으로 결합 에너지는 -22.337 kJ/mol로써, 베타 시트의 우수한 기계적 물성의 원인은 반데르발스 결합 및 polar, 분자 기계적 에너지를 통해 확인하였다. 따라서, 실크 단백질 기반 나노 소재의 기계적 물성에 미치는 중요한 요인은 반데르발스 결합과 소수성 작용인 것을 추가적으로 알 수 있다.

Fig. 3 
Mechanical properties of silk peptide beta-sheet