Journal of Computer Chemistry, Japan
Online ISSN : 1347-3824
Print ISSN : 1347-1767
ISSN-L : 1347-1767
総説
プロトンの水和自由エネルギー: 酸解離定数および標準水素電極電位の高精度計算
松井 亨喜屋武 茜庄司 光男重田 育照
著者情報
ジャーナル フリー HTML

2016 年 15 巻 5 号 p. 184-191

詳細
Abstract

本総説では,分極誘電体モデルと量子化学計算に基づく酸解離定数(pKa)の計算手法について述べる.この手法では,参照分子に対して計算により得られる自由エネルギー差と実験により得られるpKa値から官能基毎の線形関係を導くことで,アミノ酸の側鎖のpKaの半定量的な計算が可能になる.ペプチド3量体の計算においては,周囲のアミノ酸の水素結合の効果によってpKaが単量体とくらべ大きく(3 pKa単位)異なることがわかる.また本手法は標準水素電極電位の計算にも適用可能であり,いくつかの酸化還元反応の誤差はCCSD (T)/aug-cc-pVDZでは0.1V以内であり,B3LYPでも同程度の精度で酸化還元電位が計算される.

1 はじめに

酵素の機能は周囲の環境,すなわち,溶媒の種類,イオン,および溶液の温度およびpHなどの状態によって支配されている.タンパク質の表面には極性をもつアミノ酸残基や荷電性アミノ酸残基が数多く存在することから,溶液のpH値に応じてアミノ酸側鎖がプロトンを受け取ったり(プロトン化),放出したり(脱プロトン)することで,タンパク質全体の構造安定性が変化する.このようなアミノ酸残基は,触媒中心やプロトン移動経路などのようなタンパク質の活性部位に位置しているので,アミノ酸がプロトン化しているかどうかを知ることは,そのタンパク質の機能を知る上で極めて重要な因子となる.

アミノ酸残基がプロトン化してるかどうかを判断するための基準は,実験的には酸解離定数(pKa)の情報に基づいている.pKaは,酸や塩基による滴定実験により直接的に決定されるが,多くのプロトンを供与•受容するアミノ酸残基がある場合,その位置を特定することは困難となる.また,NMRはアミノ酸残基のpKa値の変化を測定するための最も有効な方法の一つであり,実際の小タンパク質のpKa値は化学シフトのpH依存性により決定されている [1,2,3,4,5].個々のアミノ酸残基のpKa値は精密に測定されているが,タンパク質になった時には周囲のタンパク質環境(例えば,水素結合,分極,疎水場など様々な効果の競合の結果)に大きく依存することが明らかとなっている.

理論的な側面では,経験的なのpKaの予測手法(PROPKAなど [6])が普及しており,分子動力学計算を実行する前のプロトン化状態の決定にしばしば用いられている.しかしながら,これらの方法はあくまで結晶構造に対する経験的なパラメータに基づいているため,非物理的なpKa値を算出することもあり,その値に対して信頼性が低いこともある.他にもPROPKAでは,補因子や小分子なのpKaの算出ができないことから,それらの系に関してはより現実的なモデルの適用が望まれる.例えば,定性的なpKa値は自由エネルギー摂動法を用いた分子動力学(MD)シミュレーションを実行することにより得ることができる.また,波動関数理論や密度汎関数法(DFT)のような量子化学計算を用いることで,PROPKAやMDなどにおける計算上のパラメータ依存性を回避することができる.小化合物のpKa値に関しては,量子化学計算を用いた研究例が数多くするを報告されており,その精度の詳細はHoとCooteの総説で議論されている [7].

pKa値を正しく得るためには,酸の解離に伴う溶液中の分子の自由エネルギー差を高精度に求める必要がある.Sprikらは,DFTに基づく第一原理MDから酸の解離に伴う自由エネルギープロファイルを得るための方法を提案し,水溶液中でのアミノ酸を含む多くの化合物のpKa値を得ることに成功している [8,9].この方法は,第一原理スキームに最も近いが,自由エネルギーを得るために多くの計算時間がかかるため,タンパク質のpKa値の計算には向かない.したがってその高い計算コストを避けるためには,静的な量子化学的方法に基づいた自由エネルギー計算手法の適用が不可欠となる.近年,溶媒効果に関しては,様々な量子化学計算手法と分極連続体モデル(PCM)との混合手法が開発され,高精度溶媒和ギブズエネルギー計算により水溶液中の脱プロトン化の傾向を議論することが可能となっている [10,11,12,13,14,15].その精度は鷹野とHoukによって報告されているように,実験値と極めてよく一致する [16].

しかしながら第一原理理論を用いたpKa計算においても,従来法ではいくつかの問題があることも分かっている.そのうちの1つは,方法論の選択(量子化学計算手法,基底関数,PCMの空洞モデルとパラメータ)であり,その選択によって誤差は大きく変わる.また,熱力学サイクルを用いて水溶液中での反応ギブズエネルギーを計算する際には気相中の計算を介するが,気相中ではアミノ酸の双性イオンの最適化構造は不安定であることから,計算することができくなるという問題もある.最も取り除くのが困難な問題点は,プロトンの水和自由エネルギー値である.従来法では,いくつかの基準値(例えば,高度計算の値や実験値)を採用し,pKaの計算において一定であると仮定する.プロトンの水和自由エネルギーの値そのものは,プロトンが電子を持たないことから計算できない.その代わりに,H3O+やさらに大きなプロトン溶媒和クラスターの解離反応やその他のプロトン交換反応からプロトンの水和自由エネルギーを算出するが,クラスターのサイズや交換反応の選択には恣意性がある.また系の増大に伴い構造の多様性も生ずることから,プロトンの水和自由エネルギーを高精度に求めることは,極めて困難となる.そのため,これまで多くの場合,定量的にpKa値を再現することができなかった.

この問題を克服するために我々は,水溶液中での酸解離ギブズエネルギー差の計算値とpKaの実験値に関して線形フィッティングを適用し,その検量線から未知化合物のpKa値を計算する新しい方法を提案した.その際,化学的な性質の違いに着目し,官能基ごとの検量線を作成することで,より高精度な計算を実現している.例えば,5-置換ウラシルのpKaの誤差が0.2のpKa単位に収まるなど [17],従来法の誤差(1∼2のpKa単位)を大きく凌駕する結果を得た [17,18,19,20].

また,本手法を用いることで,同じくその算出にプロトンの水和自由エネルギーを用いる必要がある標準水素電極電位(プロトンの解離エネルギーを基準とする酸化還元電位)の高精度計算にも成功している [21,22,23].

本総説においては,近年我々の提案しているpKaおよび標準水素電極電位の統一的計算手法に関して解説する.

2 pKaおよびEredoxの新規計算手法

pKaは通常水溶液中での酸解離反応   

HA H + + A
において,ギブズエネルギー差   
Δ G 0 = G aq ( A ) G aq ( HA ) + G aq ( H + ) (1)
を用いて以下のように書き下せる.   
p K a = Δ G 0 2.303 R T (2)

ここで,Rは気体定数, Tは絶対温度である.本研究ではプロトンの溶媒和自由エネルギーを最初から分からないものとして定数Gaq(H+)=GHとおく.また水溶液中の化合物のギブズエネルギーは 連続誘電体モデル(PCM)と振動計算により導出可能とし,ギブズエネルギーに誤差を補正するスケーリングファクターsをかけて,   

k = s 2.303 R T , C 0 = s 2.303 R T G H
とおくと(2)式は   
p K a = k ( G aq ( A ) G aq ( HA ) ) + C 0 (3)
と近似できる.実験で得られたpKaと計算によるAとHAのギブズエネルギー差は線形関係にあるので,線形フィッティングによってk, C0を決定できる.他の化合物に関しては,この式に計算によるギブズエネルギー差を代入することにより,pKaの予測値を算出する.ここで,プロトンの溶媒和自由エネルギーGHは,2つのパラメータの比    GH = C0/k    より算出できる [19,21].

一方,酸化還元電位Eredoxは   

E redox = G aq ( ox ) G aq ( red ) F E SHE
で与えられる.ここで,第1項は酸化体(ox)と還元体(red)の水溶液中でのギブズエネルギー差,Fはファラデー定数,ESHEは標準水素電極電位であり,GHおよびギブスエネルギーGgas(H2)気相中での水素分子のgを用いて   
E SHE = G H G gas ( H 2 ) / 2 F
と表すことができる.したがって,様々な計算レベルや官能基に対してGHを求めることにより,統一的にpKaEredoxが計算可能である [21].

この方法を様々な系に適用するために,特定の官能基を持つ化合物の実験値を利用した.本研究では,特に断らない限りCPCM/UFFを用い,Gaussian09を用いてB3LYP/6–31++G (d,p)レベルの計算を行っている [24].

3 結果

この理論を用いることで,特定の官能基を持つ化合物に対して半定量的にpKa値を導出することができる.今回はその例をいくつか示す.

3.1 検量線の準備

この方法を様々な酸•塩基に対して実行するため,まずはpKaの実験値がある小分子(24種類)に関して自由エネルギー差を計算し,pKa値との相関を取った [25,26,27,28,29,30].Figure 1に示されるように,全ての官能基に対して直線でフィッティングできることがわかる(R2因子 = 0.99).しかしながら,官能基によって傾きや切片の大きさが異なる.つまり,従来1つであるはずのGHがどの官能基から放出されるかによって値が異なることを意味する.しかも,これらの値は採用する量子化学研究手法やPCMの設定,基底関数によっても異なる.よって,本手法は官能基や手法ごとに実行する必要がある.

Figure 1.

 Correlation between experimental pKa values and calculated Gibbs energy for 24 compounds.

3.2 サリチル酸における水素結合の効果

複数のプロトンサイトを持つサリチル酸において,分子内水素結合がpKaに与える影響を考えるため,Figure 2にあるような2つの構造を考えてCOOH, フェノールのプロトンが順にとれる二段階の脱プロトン反応において,プロトンがある状態とない状態のギブズエネルギーの差からpKaを見積もる.

Figure 2.

 Conformations of salicylic acid (gray, red, and while indicate carbon, oxygen, and hydrogen atom, respectively. In structure (a), OH group forms a hydrogen bond to COOH group, while in structure (b) there exists no hydrogen bond.

計算の結果,分子内の水素結合の影響によりaの方が6 kcal/mol程度安定である.aの構造ではpKaは2.69と13.29, bの構造ではpKaは4.10, 10.76と算出された.実験値は2.81と13.4であることから,aの構造を支持している.abで大きく値が異なる理由はCOOHが抜けることでCOOとなる際に,aではOHがCOOと水素結合を強化しやすい形になり,脱プロトン化した状態も安定になる.この形になった場合フェノールのプロトンが脱離しにくくなり,第2段階目でのpKaが大幅に上昇する原因となる.一方で,bの構造では各々が独立しているために,安息香酸(pKa = 4.20)とフェノール(pKa = 10.02)と似た値になる.分子内の水素結合により,pKa値が2.5程度変化しうることがこれらの結果から示唆される.

3.3 アミノ酸のpKa

Table 1は本計算手法をアミノ酸側鎖に適用した結果である.Asp, Glu, Tyrにおいては文献値を0.1pKa単位以内の誤差で定量的に再現している.その一方で,Cys, His, Lysでは文献値と0.4–0.5pKa単位の誤差と多少大きくなってはいるが,従来法の誤差(1 pKa単位以上)に比べ,極めて正確にpKaを予測することが可能である.

Table 1.  Calculated pKa values of side chains in amino acid monomers and the corresponding experimental values [31-33]. Calculated pKa values of side chains of X in amino acid trimer GXG (G=Gly, X=Cys, Ser Tyr).
Group pKaCalc. pKaExp.
Asp -COOH 3.78 3.86
Glu -COOH 4.23 4.25
Cys -SH 8.85 8.33
His Imidazole 6.49 6.04
Lys -NH3+ 10.92 10.53
Ser Alcohol –OH 15.92 -
Tyr Phenol –OH 10.07 10.07
GCG -SH 8.47
GSG Alcohol –OH 12.80
GYG Phenol –OH 10.01

次に,タンパク質環境そのものを検証する前に,小さなポリペプチドの系におけるpKaを算出した.それぞれの構造は直鎖上の3量体をモデリングし, 構造最適化,および,振動解析によって自由エネルギーを算出している.Gly-Cys-Glyなどのトリペプチドにおいては,CysのpKaが8.44に下がることから,3量体を形成することにより,アミノ酸の構造からのズレや電荷や分極方向の変化がpKaに影響することがわかる.さらに,Gly-Ser-Glyにおいては,単量体よりも3 pKa単位減少している.これは,SerのOH基がアミノ酸骨格と水素結合しているためである.このように,SH基とOH基の水素結合能もpKaに大きく影響する.一方で,Gly-Tyr-Glyは側鎖が大きく他と相互作用しにくいことから,大きな影響は見られない.

3.4 Trpケージの側鎖のpKa

この手法を小タンパク質であるTrpケージ (20アミノ酸残基,構造をFigure 3に示す [34])に適用した.このサイズのタンパク質では,振動解析が事実上困難であることから,ギブズエネルギー差の代わりに電子エネルギー差を用いて,パラメータを再フィッティングした(R2因子 = 0.97).また,計算においてはpKaを算出する際に着目するアミノ酸残基に関しては6–31++G (d,p)を用い,その他のアミノ酸残基については6-31G (d)を用いた.

Figure 3.

 Structure of Trp-cage by NMR measurement (PDBID:1L2Y). Dotted circles indicates the protons in Tyr and Ser of interest.

Trpケージにはチロシンが1残基,セリンが2残基含まれている.そこで,B3LYPレベルの計算で最適化しTyr3のpKa値を導出したところ,10.02 となり単体のチロシンの側鎖と同様の結果となった.これはチロシンがタンパク質の外側を向いていて,タンパク質を構成する他のアミノ酸残基と 水素結合を形成していないためである.一方で,シニョリンには環境の異なる2つのセリン(Ser13とSer14)が存在する.Ser13のpKa値は20.80,Ser14は14.69であった.前者はタンパク質の内側にあり,単体のセリンの値(15.92)よりも極めて大きくなることがわかった.これは,Ser13の側鎖がタンパク質の他の部位と水素結合しているためである.一方で後者は,単体のものより若干小さいが,こちらの場合はSer14の側鎖がタンパク質の外側を向いており水素結合はしていないものの,タンパク質全体による分極によってpKa値が変化しているものと考えられる.

本解説では小タンパク質に対して全電子状態計算を用いた解析を行っている.さらに大きなタンパク質のpKaを予測するためには,QM/MM法 [35]やフラグメント分子軌道法 [36]など,大規模系を取り扱う手法が必須となる.また,水チャネルやプロトン移動経路など,タンパク質の内部では誘電環境も水溶液中と異なることから,溶媒和モデルにおいても改良が必要である.またその際には,パラメータの決定などにも工夫が必要である.したがって,さらに大きなサイズのタンパク質の生体内機能に重要な部位のpKaの定量的な評価は,この分野における大きなチャレンジであり,是非,克服されるべき課題であるとといえる.

3.5 pKaと標準水素電極電位(SHE)に対する計算手法依存性

これまでの計算はタンパク質などの大規模系への適用を想定して,全てB3LYPを用いたものであった.ここでは,同様の計算に関してその手法依存性を検証すると共に,もう一つの物理量である酸化還元電位,すなわち標準水素電極電位(SHEポテンシャル)の計算精度に関して議論する.

典型的な第1級アルコール分子に関して,pKa値(実験値)と計算の脱プロトンエネルギーの計算手法依存性をFigure 4に示す [21].Figure 4にあるように, アルコール分子においては全ての計算手法で線形関係にあることが分かる.特に,摂動論の次数を上げることで系統的に精度が向上する.また,近似直線の切片から算出されたプロトンのギブスエネルギーとSHEポテンシャルをTable 2に示す.Gaussian-3 (G3)法で−265.5 kcal/molと,実験値(−265.9 kcal/mol)に極めて近い結果が得られた.一方で,SHEポテンシャルは,計算手法によって大きく変化するが,Coupled cluster法(CCSD (T))やG3法などの高精度計算ではIUPACが推奨する実験値(4.44V)に極めて近い値を算出している.一方で,Hartree-Fock (HF)法では3.83Vと実験値よりも0.6V程度低い値を出しているが,これはGgas(H2)を計算する際,電子相関が欠如しているために,過剰に低く見積もられているためである.

Figure 4.

 Methodology dependence of pKa values and Gibss energies. Hartree-Fock (HF), Møller-Presset perturbation theory (MPn), and Coupled cluster (CC) results are shown.

Table 2.  Methodology dependence of calculated and experimental standard hydrogen electrode potential (V) and Gibbs energy (kcal/mol).
Method SHE (V) GH (kcal/mol)
B3LYP/6–31++G** 4.8 −259.75
HF/CBS-Limit 3.83 −268.35
CCSD (T)/aug-cc-pVDZ 4.52 −265.17
G3B3 4.53 −265.46
Exp. (IUPAC) 4.44 −265.9

これを利用して,小分子における8種類の酸化還元反応電位の実験値 [37]との絶対誤差はFigure 5のようになった.SHEポテンシャルを調整しなければ,B3LYPやHFにおいては大きな誤差(それぞれ0.26Vおよび0.48V)が生じている.一方で,CCSD (T)では平均誤差が0.1V程度とこの時点でも正確な値が出ている.ここで,Table 2で得られたSHEポテンシャルを適用するとB3LYPでの平均誤差は0.08Vとなり,CCSD (T)などの高精度な計算(0.08V)と同程度の誤差に収まる.また,HFにおいても誤差は大幅に減少している(0.17V)ことから,HFやDFTを用いて酸化還元電位を算出する場合,SHEポテンシャルの補正が不可欠といえる.

Figure 5.

 Absolute error between calculated and experimental standard hydrogen electrode potential (V) for 8 redox reactions using experimental SHE=0.44V and methodology dependent calculated SHE (see Table 2).

4 まとめ

本研究では,実験値のpKa値と量子化学計算から得られるプロトン解離に伴う溶液中の自由エネルギー差から官能基(や計算手法,基底関数,溶媒和モデル)に依存する検量線を作成し,未知の化合物のpKa値を予測する手法に関して解説した.従来法では,計算が困難なプロトンの水和自由エネルギーに対して実験値を採用するが,本手法ではそれらを官能基ごとフィッティングによって求めている.ハイブリッド密度汎関数法を用いた本手法の精度は,従来法での1∼2 pKa単位に比べて大きく改善し,最大で0.5 pKa単位となる.本手法をサリチル酸,アミノ酸,タンパク質な度に適用し,その局所構造依存性に関して議論した.

酸化還元反応で用いられる標準水素電極電位に関しても,同様の手法を適用することで,高精度に酸化還元電位を求めることが可能となる.その精度は,CCSD (T)などの高精度計算ではIUPACが推奨する実験値を用いても,本手法を適用しても同程度であるが,ハイブリッド密度汎関数法に関しては大きく改善し,CCSD (T)と同程度の誤差となった.これらの手法は,遷移金属錯体やさらに大きな生体内分子の半反応の酸化還元電位の計算にも用いられており,非常に汎用性が高い.

Acknowledgment

本研究は著者ら(主に松井,重田)が兵庫県立大学に在籍していた時に始め,所属がバラバラになりながらも,大阪大学,筑波大学と所属を変えながら続けた成果をまとめたものです.これらの研究に関わりました,全ての共同研究者に感謝申し上げます.

参考文献
 
© 2016 日本コンピュータ化学会
feedback
Top