💊
分子動力学シミュレーション解析であるwaterswappingを論文に沿ってやってみた。

2025/01/06に公開
gromacs
本記事はリガンド-タンパク質 MDsimulation（分子動力学シミュレーション）を論文に沿って行なった結果に対し、

論文で使われている解析方法であるwaterswappingを行なっています。
動作検証済み環境

Windows 11 Home, 13th Gen Intel(R) Core(TM) i7-13700, 64 ビット オペレーティング システム、x64 ベース プロセッサ, メモリ：32GB
本記事は以下の論文のフォローになります。

※本論文のMD simulationはAMBERを用いて行なっておりますが、本記事ではGROMACSでフォローしています。
https://www.nature.com/articles/s41598-022-10364-z
本論文は標的の選定方法からスクリーニング、スクリーニングした薬物候補化合物の毒性、物性予測、MDシミュレーションまでしており、ほぼ無料のツールですることができます。in silico創薬の基本が詰まっているので、一緒に勉強していきましょう。

 宣伝本記事を見てくださり、ありがとうございます。
インシリコ創薬についてより学びたい方は
拙著

https://zenn.dev/labcode/books/xvsw0lszpqir77

で学び、さらに色々な方法で新薬探索を楽しんでいただければと思います！
また化合物の評価を行いたい場合は

https://zenn.dev/labcode/books/ar7sn7jlhhafgm

を見ていただければと大変嬉しいです。

 In silico創薬とその流れ
 In silico 創薬、in silicoスクリーニングとはインシリコ創薬（in silico drug discovery）は、コンピュータを駆使し、シミュレーションやデータ解析を用いて新薬を設計・発見する手法です。
これにより、従来の実験的な方法に比べてコストと時間を大幅に削減できるとされています。
その中でもin silicoスクリーニングは、膨大なデータベースから有望な候補物質を迅速に特定し、その効果や副作用を予測することで、創薬の初期段階での効率化が図られます。
一見難しそうなin silico創薬ですが、現在では様々なアプリケーションやwebサイトがあり、それらを駆使すれば、誰でも簡単に創薬をすることもできます。本記事では、それらのアプリケーションを駆使し、in silico創薬を行った論文をもとに、手法をわかりやすく説明していきます。

 In silico創薬の流れ
ターゲット選定、準備：疾患の原因となる分子（ターゲット）を特定し、その3D構造を準備します。

化合物ライブラリの構築：in silicoスクリーニングに使用する化合物集団（＝ライブラリ）を作成します。

in silicoスクリーニング：コンピュータ上で数百万の化合物を対象に、ターゲット分子との結合親和性を評価します。

分子動力学シミュレーション（MD simulation)：候補化合物とターゲットの相互作用を動的に解析し、安定性や効果を予測します。

4.1.環境構築:Windowsでの環境構築を示しています。

4.2.MD simulation:MD simulationを行います。結果の可視化を行います。

4.3.RMSD, RMSF, RoG, hydrogen bond, RDF 解析:構造変化や分子間相互作用を評価するための解析手法を行います。

4.4.MM-GBSA解析とMM-PBSA解析、decomposition解析、アラニンスキャン:分子間相互作用や結合エネルギーへの各残基の寄与を解析します。

4.5.Water swapping解析（本記事):リガンド-タンパク質間の結合自由エネルギー計算を行います。npx zenn preview

物性、毒性評価：過去のデータを用いて、新たな候補物質の特性や副作用を予測するモデルを構築します。

（実験的検証）：in silico解析で得られた候補物質を実験室で合成し、実際の効果や安全性を検証します。
すでに公開したものについては、リンクを貼っています。
本記事はFigure 9のWaterswappingにおけるBAR, FER, TIの値を測定します。（詳細な説明は後述）

 Waterswapping解析とは
 1. MM-PBSAおよびMM-GBSA法の限界
以前ご紹介したMM-PBSAやMM-GBSA法は、シミュレーションスナップショットから自由エネルギーを計算する手法です。しかし、これらにはいくつかの課題があります。まず、シミュレーション中のリガンド構造の変化により、どの部分が結合エネルギーに寄与しているかを明確に特定するのが難しいという点です。また、暗黙的溶媒モデルを使用するため、リガンドとタンパク質間で水分子が橋渡しをする効果を考慮できないという限界があります。

 2. WaterSwap法の仕組み
WaterSwap法は、リガンドと水分子を入れ替え（スワップ）し、その過程でのエネルギー変化を基に自由エネルギーを計算する手法です。この方法の利点は、水分子の役割を直接モデル化できる点にあります。また、計算の収束性が高く、安定した結果を得ることが可能です。さらに、リガンドのどの部分が結合エネルギーに寄与しているかを精密に評価することができます。

 3. WaterSwap法のメリット
WaterSwap法の最大のメリットは、収束性と安定性に優れている点です。他の手法に比べて、結合エネルギーが安定して収束するため、信頼性の高い結果が得られます。また、水分子の橋渡し効果を含めた詳細な結合解析が可能で、複雑なリガンドや水分子の関与が重要な系にも適用しやすい手法です。

 OpenBioSimの環境構築まずWaterswappingが実行できるツールであるOpenBioSimについて、インストラクション（2. Easy installation - Run in a conda environment）に従って環境構築していきます。
全コードはこちら
# miniconda インストール
curl -L -O "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh"
bash Miniforge3-$(uname)-$(uname -m).sh

# 実行環境作成
conda create -n openbiosim "python<3.10"
conda install -n base -c conda-forge mamba
mamba install -n openbiosim -c openbiosim biosimspace
conda activate openbiosim
python -m ipykernel install --user --name=openbiosim --display-name "Python (openbiosim)"

#必要に応じて以下も実行
sudo apt update
sudo apt install libgsl-dev
conda install -n openbiosim -c conda-forge gsl

#実行確認
import BioSimSpace as BSS
!waterswap --help

以下に一つずつ解説していきます。
まず以下でminigforgeをインストールしてください。
curl -L -O "https://github.com/.../Miniforge3-$(uname)-$(uname -m).sh"

curl:
Webからデータを取得するためのコマンドラインツールです。
URLを指定すると、そのデータをダウンロードできます。


L:
リダイレクトを追跡するオプションです。
指定した URL がリダイレクトされている場合でも、自動的に正しい場所にアクセスします。


O:
ダウンロードしたデータをファイルとして保存するオプションです。
ファイル名はサーバー側で指定された名前になります。

URL部分 ("https://github.com/.../Miniforge3-$(uname)-$(uname -m).sh"):

uname:
システムの情報（オペレーティングシステムの名前）を取得するコマンド。

$(uname) は現在のシステム（例: Linux, Darwin など）を動的に取得します。


uname -m:
マシンのアーキテクチャ情報を取得します（例: x86_64, arm64）。

$(uname -m) は現在のシステムのアーキテクチャ（CPUの種類）を動的に取得します。


この URL は、現在使用しているシステムに対応する Miniforge インストーラスクリプトを自動的に選択します。
bash Miniforge3-$(uname)-$(uname -m).sh

bash:
シェルスクリプトを実行するためのコマンド。
ダウンロードした Miniforge3-$(uname)-$(uname -m).sh を実行します。


Miniforge3-$(uname)-$(uname -m).sh:
上でダウンロードしたインストーラファイル名。

$(uname) と $(uname -m) が展開され、具体的なファイル名になります（例: Miniforge3-Linux-x86_64.sh）。

このコマンドを実行すると、インストーラが起動し、Miniforge のインストールが進みます。
以下をターミナルで実行し、実行環境を作ります。
（mambaというcondaよりも高速なものを使っています。）
conda create -n openbiosim "python<3.10"
conda install -n base -c conda-forge mamba
mamba install -n openbiosim -c openbiosim biosimspace
conda activate openbiosim
python -m ipykernel install --user --name=openbiosim --display-name "Python (openbiosim)"
conda create -n openbiosim "python<3.10”

conda create:
新しい仮想環境を作成するコマンド。


n openbiosim:
作成する仮想環境の名前を openbiosim に指定します。


"python<3.10":
仮想環境内で使用する Python のバージョンを指定します。
この場合、Python 3.10 未満（例えば、3.9.13 など）がインストールされます。
BioSimSpace は特定の Python バージョンを要求するため、この指定が重要です。

conda install -n base -c conda-forge mamba

conda install:
パッケージをインストールするコマンド。


n base:
Conda の「ベース環境」にインストールすることを指定します。


c conda-forge:
Conda-Forge チャネル（コミュニティで管理される信頼性の高いリポジトリ）を指定。
Mamba は Conda-Forge で配布されています。


mamba:
Conda の代替ツールで、高速なパッケージ管理を提供します

mamba install -n openbiosim -c openbiosim biosimspace

mamba install:
パッケージをインストールするコマンド（Mamba 版）。


n openbiosim:
仮想環境 openbiosim にインストールすることを指定。


c openbiosim:
BioSimSpace を配布している OpenBioSim チャネルを指定します。


biosimspace:
分子シミュレーション用のライブラリです。
化学計算や分子シミュレーションを Python で簡単に扱えるようにするパッケージ。

conda activate openbiosim

conda activate:
指定した仮想環境をアクティブ化するコマンド。
アクティブ化後、その環境内の Python やパッケージが使用されます。


openbiosim:
先ほど作成した仮想環境の名前。

python -m ipykernel install --user --name=openbiosim --display-name "Python (openbiosim)”

python -m ipykernel install:
Jupyter Notebook 用に仮想環境をカーネルとして登録するコマンド。


-user:
ユーザーレベルで設定を適用します（システム全体に影響を与えません）。


-name=openbiosim:
仮想環境の内部名を指定します。
この名前は Jupyter が認識するカーネルの名前になります。


-display-name "Python (openbiosim)":
Jupyter Notebook 上で表示される名前を指定します。
ここでは "Python (openbiosim)" という名前で選択可能になります。

vscodeで環境を設定します。
右側の赤文字と記載している箇所をクリックすると、真ん中で仮想環境が設定できるようになるので、設定してください。
（openbiosimがなければ、Select Another Kernel...→Python environments…選択できるようになります。）
以下がキチンと実行できれば、環境構築は完了です。
import BioSimSpace as BSS
（私の環境の場合、追加で以下も入れました）
sudo apt update
sudo apt install libgsl-dev
conda install -n openbiosim -c conda-forge gsl

sudo apt update:システムのパッケージリストを最新の状態に更新します。

sudo apt install libgsl-dev:GNU Scientific Libraryの開発パッケージ（libgsl-dev）をインストールします。

conda install -n openbiosim -c conda-forge gsl:Anacondaの仮想環境「openbiosim」にGSLをインストールします。
続いて、waterswappingが使えるか確認します。
!waterswap --help
以下のように出力されれば、waterswapは使えます。
Starting waterswap: number of threads equals 24

usage: waterswap [-h] [--description] [-H] [--author] [--version]
                 [-l [LIGAND]] [-t [TOPOLOGY_FILE]] [-c [COORDINATE_FILE]]
                 [-C [CONFIG]]
                 [--lambda_values LAMBDA_VALUES [LAMBDA_VALUES ...]]
                 [-n [NUM_ITERATIONS]] [-f]

Calculate absolute binding free energies using waterswap

optional arguments:
  -h, --help            show this help message and exit
  --description         Print a complete description of this program.
  -H, --help-config     Get additional help regarding all of the parameters
                        (and their default values) that can be set in the
                        optionally-supplied CONFIG file
  --author              Get information about the authors of this script.
  --version             Get version information about this script.
  -l [LIGAND], --ligand [LIGAND]
                        Supply the name of one of the residues in the ligand
                        whose binding free energy is to be calculated. By
                        default, the ligand will be the first non-protein,
                        non-solvent molecule in the input topology file.
  -t [TOPOLOGY_FILE], --topology_file [TOPOLOGY_FILE]
                        The Amber topology file containing the solvated
...

waterswap is built using Sire and is distributed under the GPL. For more
information please visit https://sire.openbiosim.org/waterswap.html, or type
'waterswap --description'
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...

 GROMACSファイルをAMBERファイルに変換全コードはこちら。以下に解説していきます。
 conda install -c conda-forge ambertools=23
 pdb4amber -i protein_A_chain.pdb -o fixed_protein.pdb --reduce
 antechamber -i ligand.mol2 -fi mol2 -o ligand.mol2 -fo mol2 -c bcc -at gaff
 parmchk2 -i ligand.mol2 -f mol2 -o ligand.frcmod
 ~tleap.inファイルの作成~
 tleap -f tleap.in
以下に解説していきます。
 conda install -c conda-forge ambertools=23
まずAMBERTOOLSをインストールします。

conda:
Pythonやそのライブラリを管理するためのパッケージ管理ツールです。
仮想環境の作成、依存関係の解決、パッケージのインストール・更新を行います。


install:
指定したパッケージをインストールするためのサブコマンドです。
ここではambertools=23をインストールします。


c conda-forge:

conda-forgeというチャンネル（リポジトリ）からパッケージを取得します。


ambertools=23:
AmberToolsのバージョン23をインストールします。

まずProtein.pdbをpymolなどを用いて、A鎖だけにして、Protein_A_chain.pdbにしてください。
waterswappingはAMBER対応なので、以下のコードでAMBER形式に変更します。
pdb4amber -i protein_A_chain.pdb -o fixed_protein.pdb --reduce

pdb4amber:
AmberToolsの一部で、PDB形式の分子構造データをAmberフォーマットに適合させるためのツール。
入力PDBファイルの検証、修正、水素原子の追加、ジスルフィド結合の識別などを行います。


i protein_A_chain.pdb:
入力ファイル（protein_A_chain.pdb）を指定します。
ここでは、protein_A_chain.pdbが処理対象のPDBファイルです。


o fixed_protein.pdb:
出力ファイル名を指定します。
修正後のPDBファイルがfixed_protein.pdbとして保存されます。


-reduce:

ジスルフィド結合の自動検出と修正、および欠落している水素原子の追加を行うオプション。

リガンドのファイルもAMBER用に合わせていきます。
antechamber -i ligand.mol2 -fi mol2 -o ligand.mol2 -fo mol2 -c bcc -at gaff
parmchk2 -i ligand.mol2 -f mol2 -o ligand.frcmod
antechamber -i ligand.mol2 -fi mol2 -o ligand.mol2 -fo mol2 -c bcc -at gaff
i ligand.mol2
入力ファイルligand.mol2を指定します。
fi mol2
入力ファイルの形式をmol2形式として指定します。
o ligand.mol2
修正・処理された出力ファイルをligand.mol2という名前で保存します。
fo mol2
出力ファイルの形式をmol2形式に指定します。
c bcc
BCC（Bond Charge Correction）法を用いて分子の部分電荷を計算します。
at gaff
General Amber Force Field（GAFF）の原子タイプを各原子に割り当てます。
parmchk2 -i ligand.mol2 -f mol2 -o ligand.frcmod
i ligand.mol2
入力ファイルligand.mol2を指定します。
f mol2
入力ファイルの形式をmol2形式として指定します。
o ligand.frcmod
出力ファイルをligand.frcmodという名前で保存します。このファイルにはAmberシミュレーション用の不足パラメータが補完されます。
以下のtleap.inファイルを作成してください。
tleap.in
source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p

# タンパク質をロード
protein = loadPdb "fixed_protein.pdb"

# リガンドをロード
loadamberparams ligand.frcmod
ligand = loadMol2 "ligand.mol2"

# タンパク質とリガンドを結合
complex = combine {protein ligand}

# 水ボックスの追加
solvateBox complex TIP3PBOX 10.0

# イオンを追加して中性化
addIons complex Cl- 0
addIons complex Na+ 0

# パラメータを保存
saveAmberParm complex complex.top complex.inpcrd
savePdb complex complex.pdb

quit

このtleapスクリプトは、Amber分子動力学シミュレーション用のシステムセットアップを行うためのものです。以下に各行の役割を詳しく説明します。
力場ファイルのロード

source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p

leaprc.protein.ff14SB: タンパク質用のff14SB力場をロードします。ff14SBは、Amberで広く使用される力場の1つです。

leaprc.gaff: GAFF (General Amber Force Field)をロードします。これはリガンドや非標準分子に対して使用されます。

leaprc.water.tip3p: TIP3P水モデルをロードします。分子動力学シミュレーションにおける標準的な水分子モデルです。
タンパク質をロード
protein = loadPdb "fixed_protein.pdb

loadPdbコマンドで、fixed_protein.pdbファイルからタンパク質構造を読み込みます。
読み込んだ構造をproteinという変数に格納します。
リガンドのパラメータと構造をロード
loadamberparams ligand.frcmod
ligand = loadMol2 "ligand.mol2"

ligand.frcmod: リガンド用の力場パラメータファイル (frcmod形式) を読み込みます。
GAFFで生成されたリガンドのパラメータを含んでいます。


ligand.mol2: リガンドの構造を含むmol2形式のファイルを読み込み、ligandという変数に格納します。
タンパク質とリガンドの結合
complex = combine {protein ligand}
*combine*コマンドで、読み込んだタンパク質とリガンドを1つのシステム (複合体) に結合します。
結合したシステムをcomplexという変数に格納します。
水ボックスの追加
solvateBox complex TIP3PBOX 10.0
**solvateBox**コマンドで、complexに水分子を追加して溶媒環境を作成します。

TIP3PBOX: TIP3P水モデルを用いた直方体の水ボックスを指定します。

10.0: システムの各方向に10.0 Åの水分子を追加することを意味します。
イオンを追加して中性化
addIons complex Cl- 0
addIons complex Na+ 0

addIonsコマンドで、システムにイオンを追加します。

Cl-およびNa+を指定して、システムの電荷を中性化します。

0は、ここでは特定の数値ではなく、「中性化するために必要な数のイオン」を自動計算するよう指定しています。
パラメータと構造の保存
saveAmberParm complex complex.top complex.inpcrd
savePdb complex complex.pdb

saveAmberParm:

complexシステムのパラメータおよび初期構造を保存します。

complex.top: Amberのトポロジーファイル (力場パラメータや接続情報)。

complex.inpcrd: 初期座標ファイル。


savePdb:

complex.pdb: 結合したシステムをPDB形式で保存します。

tleapの終了
quit
*quit*コマンドでtleapを終了します。
以下で作成したtleapを実行していきます。
tleap -f tleap.in

tleap: AmberToolsの一部であるシステムセットアップツール。

f: ファイルモードを指定し、tleapに入力ファイルを渡します。

tleap.in: tleapスクリプトファイル (テキスト形式) で、システム構築のための命令が記述されています。

 Waterswapping解析実行以下でwaterswappingを実行していきます。終わるのに約2日かかりました。
waterswap -t complex.top -c complex.inpcrd -l UNL

waterswap: 結合自由エネルギーを計算するための手法。リガンドと水分子を「スワップ」することで、リガンドの結合親和性を評価します。

t complex.top: トポロジーファイルを指定します。
complex.top: Amberで生成された複合体 (タンパク質 + リガンド + 水分子) の力場パラメータおよび結合情報が含まれたファイルです。


c complex.inpcrd: 初期構造 (座標) ファイルを指定します。
complex.inpcrd: Amberで出力された座標情報を持つファイル。最初のシミュレーション条件 (エネルギー最小化や平衡化後) を示します。


l UNL: リガンドの残基名を指定します。
UNL: リガンドの識別名 (Amberフォーマットに従い指定されます)。
リガンドがcomplex.topおよびcomplex.inpcrdに含まれていることが前提です。


 結果終了すると、output ディレクトリが生成しています。
この中に各種ログや結果が格納されています。
次のコードを実行してください。
analyse_freenrg -i output/freenrgs.s3
このコマンドは、WaterSwapやその他の自由エネルギー計算プログラムで生成された自由エネルギー解析結果を処理・分析するために使用されます。

analyse_freenrg:
自由エネルギー解析を行うプログラム。
WaterSwapや関連ツールで実行された計算結果を解析する際に使われます。


output/freenrgs.s3:
入力ファイルのパスとファイル名。

output/ディレクトリにあるfreenrgs.s3というファイルを解析対象に指定しています。
このファイルには、WaterSwap計算で得られた自由エネルギーのデータが格納されています。

以下のように結果が出てきます。
# Free energies 
# Bennetts = -16.339873459301437 +/- 0.019185393764403726 kcal mol-1#
# FEP = -15.185634029041388 +/- 9.27384986669758 kcal mol-1#
# TI = -14.766822185190964 +/- 1.4769497989874303 kcal mol-1 (quadrature = -15.205205079031563 kcal mol-1)#
1. Bennett’s Acceptance Ratio (BAR)
BAR法は、状態間の自由エネルギー差を計算するための高精度な方法です。この手法は、状態間のエネルギー重なり（例: リガンドが結合している状態と結合していない状態の重なり）を利用することで効率的に計算を行います。具体的には、逆方向（例:A→B とB→A）のエネルギー情報を統合することで、少ないサンプリングでも非常に収束性の高い結果を得ることができます。そのため、BAR法は最も信頼性の高い自由エネルギー計算手法とされています。

Bennetts = -16.34 ± 0.019 kcal/mol
−16.34kcal/molは、リガンドがタンパク質に結合する際に放出されるエネルギーを表しています。負の値であるため、結合が熱力学的に有利である（安定した結合を形成する）ことを示します。
誤差は±0.019 kcal/molと非常に小さいため、この値は高い精度で収束しており、信頼性が非常に高いです。
2. Free Energy Perturbation (FEP)
FEP法は、2つの状態間の自由エネルギー差を統計力学的に直接計算する方法です。エネルギー差
ΔUを基に、統計的平均をとることで自由エネルギー差を求めます。この方法は、計算理論がシンプルで実装が容易であり、計算コストも比較的低いという利点があります。しかし、状態間のエネルギー重なりが少ない場合や、サンプリングが不足している場合には結果が収束しにくく、誤差が大きくなるという欠点があります。そのため、FEP法は単純なシステムや初期の概算計算に適していると言えます。
FEP = -15.19 ± 9.27 kcal/mol
−15.19kcal/molは、リガンドの結合エネルギーを示しますが、誤差が非常に大きいため、この値の信頼性は低いです。

±9.27 kcal/molは、サンプリングが不十分であるか、状態間の重なりが不足している可能性を示唆しています。
3. Thermodynamic Integration (TI)
TI法は、状態間の自由エネルギー差をエネルギー傾き（∂U/∂λ）を積分することで計算する方法です。ここでλは、状態間の相互作用をスケールするためのパラメータで、0から1の範囲で変化します。この手法は理論的に明確で柔軟性が高いことが特徴です。

TI = -14.77 ± 1.48 kcal/mol (quadrature = -15.21 kcal/mol)
−14.77kcal/mol（または−15.21kcal/mol）は、リガンドの結合エネルギーを示します。他の手法（特にBAR法）と近い値を示しており、妥当な結果と考えられます。

±1.48 kcal/molは中程度の収束性を示しており、信頼性はBAR法より劣るものの、FEP法より高いです。
まとめ
全般的に-15kcal/mol付近の結合エネルギーを持っていることがわかります。これによりリガンドはタンパク質に対して高い結合親和性を持つと予測できます。

 参考文献https://github.com/CCPBioSim/xswaps-workshop/tree/master

https://sire.openbiosim.org/index.html#
宣伝

In silico創薬とその流れ

In silico 創薬、in silicoスクリーニングとは

In silico創薬の流れ

Waterswapping解析とは

1. MM-PBSAおよびMM-GBSA法の限界

2. WaterSwap法の仕組み

3. WaterSwap法のメリット

OpenBioSimの環境構築

GROMACSファイルをAMBERファイルに変換

Waterswapping解析実行

結果

参考文献

Discussion
