در کار مطالعاتی دیگری که توسط پاند^[۴۶] و شرتس^[۴۷] و همکارانش انجام شده است، انرژی­های آزاد پیوند مطلق^[۴۸] (استاندارد) هشت لیگاند مربوط به FK506 با FKBP12 (پروتئین پیوندی FK506) با بهره گرفتن از شبیه­سازی­های دینامیک مولکولی اختلال انرژی آزاد^[۴۹] با پتانسیل مرزی حلال عمومی^[۵۰] محاسبه می­ شود تا معلوم شود کارآمدی چنین راهکار محاسباتی در عمل تا چه حد است.
در اصل، شبیه­سازی­های دینامیک مولکولی FEPMD مبتنی بر مدل­های اتمی، قدرتمند­ترین و مطمئن­ترین روش جهت تخمین انرژی­های آزاد پیوند لیگاند­ها به ماکرومولکول­ها محسوب می­ شود.
بررسی محاسباتی نشان داده است که برای محاسبه­ی تمایلات نسبی اتصال در سیستم­های بیولوژیکی مهم، شبیه­سازی­های FEPMD مطمئن­تر از سایر روش­ها جواب می­دهد،[۳۲،۳۳] و می ­تواند به تأثیر حلال و انعطاف­پذیری بپردازد [۳۴]. این امیدواری وجود دارد که محاسبات مبتنی بر شبیه­سازی­های FEPMDدر برهمکنش­های پروتئین-لیگاند بتواند به ابزار مفید­ی در کشف دارو و بهینه­سازی آن تبدیل شود[۴۳-۳۵].

برای شبیه­سازی دقیق رفتار مولکول­ها، باید بتوانیم نوسانات دمایی و بر­همکنش­های محیطی بر­خاسته از سیستم­های پیچیده و متنوع مثل جایگاه اتصال پروتئین یا توده­ی محلول را توضیح دهیم. در شبیه­سازی­های FEPMD، پردازش مولکول­های حلال عموماً هزینه­ محاسباتی را تحت­الشعاع قرار می­دهد. یک روش حد­واسط که هم شامل پردازش صریح و هم غیر­صریح حلال باشد، شامل شبیه­سازی تعداد کمی مولکول حلال در مجاورت ناحیه­ی مورد نظر و همزمان بیان تأثیر حلال احاطه­کننده با یک پتانسیل مؤثر مرز حلال است[۴۹-۴۴]. چنین تقریبی یک راهکار جذاب برای کاهش هزینه­ محاسباتی محاسبات FEPMD به شمار می­رود، زیرا اختصاصی بودن اتصال، اغلب اوقات تحت تأثیر برهمکنش­های محلی در مجاورت لیگاند قرار می­گیرد، در حالی که نواحی دورتر از مولکول گیرنده، تنها به میزان کافی نقش دارد. روش به­کار­رفته در این بررسی، GSBP نام دارد. روش GSBP هم شامل میدان استاتیک ناشی از اتم­های دورتر ماکرومولکول که با حلال محافظت می­ شود و هم شامل میدان واکنش ناشی از پاسخ دی­الکتریک حلال که روی اتم­های ناحیه­ی مورد شبیه­سازی عمل می­ کند، می­باشد[۵۳].
این روش تعمیمی از پتانسیل مرزی حلال کروی است که برای شبیه­سازی یک حل­شونده در حلال آب طراحی شده است. در این روش، تمامی اتم­های موجود در ناحیه­ی درونی لیگاند، ماکرو­مولکول­ها یا حلال تحت تأثیر تشبیه­سازی قرار می­گیرد، در حالی که تأثیر اتم­های ماکرومولکول و حلالی که خارج از ناحیه­ی درونی قرار دارند، به طور غیر­صریح لحاظ می­ شود. همچنین می­توان هزینه­ محاسباتی شبیه­سازی­های FEPMD را کاهش داد و حتی دقت آن را هم با بهره گرفتن از چند ویژگی بهبود بخشید. اول، انرژی آزاد پیوند مطلق به چند مرحله­ متوالی تفکیک شده که در طول این مراحل، بر­همکنش­های احاطه­کننده­ لیگاند و پتانسیل­های مختلفی که انتقال، جهت­گیری و پیکربندی آن لیگاند را محدود می­ کنند، لحاظ یا نادیده گرفته می­ شود[۵۳].
دوم، نمونه گیری از پیکربندی لیگاند به کمک یک پتانسیل مهارکننده­ی مبتنی بر انحراف جذر میانگین مربعات نسبت به پیکربندی حالت مقید تقویت می­ شود.
سوم این که، سهم انرژی آزاد دفعی و جذب ناشی از بر­همکنش­های لنارد-جونز لیگاند با محیط اطرافش یعنی پروتئین و حلال با بهره گرفتن از جدا­سازی ویک- چاندلر- اندرسون^[۵۱] محاسبه می­ شود. این جدا­سازی همچنین همگرایی محاسبات FEPMD را بهبود می­بخشد.
چهارم این که، برای کاهش هزینه­ های محاسباتی، تنها تعداد کمی اتم در حوالی جایگاه اتصال به­ طور صریح شبیه­سازی می­شوند. در حالی که تمام تأثیرات اتم­های باقیمانده به­ طور غیر­صریح و با بهره گرفتن از روش GSBP لحاظ می­ شود.
هشت لیگاند مربوط به FK506 (لیگاندهای ۲، ۳، ۵، ۶، ۸، ۹، ۱۲، ۲۰) در شکل (۳-۳) نشان داده شده است. این لیگاندها بر اساس کار تجربی [۵۱-۵۰] شماره­گذاری شده است. لیگاند ۲۰، مولکول FK506 است [۵۲]. سه نوع مجموعه­ اولیه برای این محاسبات در نظر گرفته شده است. مجموعه اول x-rayشامل ساختار­های بلورین با لیگاند­های ۸، ۹ و ۲۰ است. مجموعه­ دوم، Mod است که به مدل­هایی برای لیگاند­های ۳ و ۵ اختصاص دارد که از طریق ساختار بلورین FKBP12 در کمپلکس با لیگاند ۹ ایجاد شده است. جایگزین­کردن گروه سیکلو­هگزین لیگاند ۹ با یک هیدروژن، لیگاند ۵ را می­دهد. در حالی که جایگزین کردن گروه فنیل متیل با یک هیدروژن، لیگاند ۳ را می­دهد. لیگاند­های ۳ و ۵، بسیار شبیه به لیگاند ۹ هستند، و مدل­سازی مستقیم قابل توجیه است. مجموعه­ سوم که توسط شرتس و پاند ارائه شده است منطبق با مختصات مدل­های داکینگ لیگاند­های ۲، ۳، ۵، ۶، ۱۲ و ساختار­های بلور لیگاند­های ۸، ۹، ۲۰ همراه با ۲۰۰ پیکو­ثانیه شبیه­سازی MD با حلال می­باشد. به این مجموعه MD گفته می­ شود.
شکل (۳- ۱ ) فرمول ساختاری هشت لیگاندی که در محاسبات مورد استفاده قرار گرفت[۵۳]
با توجه به جدول (۳-۲)، نتایج محاسبات برای هشت لیگاند با مقادیر تجربی و همچنین با نتایج محاسبات FEPMD گسترده توسط پاند، شرتس و همکارانش مقایسه شده است. این نتایج، توافق خوبی با نتایج تجربی دارد، خصوصاً در مورد لیگاند­هایی که برای آن­ها، این امکان وجود دارد که مجموعه­ اولیه کمپلکس خوبی یا از بلور­شناسی اشعه­ی X (لیگاند­های ۸، ۹، ۲۰) و یا از مدل­سازی مستقیم ساده (لیگاند­های ۳ و ۵) داشته باشد. تفاوت برای بیشتر لیگاند­ها در حد kcal.mol^-11 است که تقریباً در حد خطا­های آماری محاسباتی می­باشد.
انرژی­های آزاد پیوندی محاسبه­شده برای لیگاند خاص که از ساختار­های اولیه مختلف به­دست آمده، مثلاً MD و اشعه­ی ایکس لیگاند ۸ [۵۳] بسیار مشابه هستند. به نظر می­رسد خطا­ها برای ساختار­های بلوری در کل کوچک­تر باشد، انرژی آزاد حلال­پوشی لیگاند­ها، به جز مورد ۳ و ۵ که بر­اساس معادله­ (۳-۳۱) محاسبه شده، با آن نتایج سازگاری خوبی دارد،. در کل نتیجه می­گیریم که این شبیه­سازی بر­اساس مدل اتمی GSBP کاهش ­یافته، یک روش دقیق و کار­آمد است.
(۳-۳۱)
در معادله­(۳-۳۱)، انرژی آزاد پیکربندی لیگاندی است که از محیط اطراف خود جدا شده است و انرژی آزاد ناشی از بر­همکنش لیگاند با حلال است.
جدول (۳- ۲) نتایج محاسبات برای هشت لیگاند با مقادیر تجربی و همچنین با نتایج محاسبات FEPMD گسترده توسط پاند
جدول(۳-۲) انرژی آزاد پیوند و انرژی آزاد حلال پوشی در مقایسه با نتایج دیگر[۵۳]
۳-۲-۲- روش تدریجی
در این روش، در هر مرحله هامیلتونی، به مقدار بسیار کوچک و ثابتی تغییر داده می­ شود. این بدان مفهوم است که در هر مرحله اختلاف هامیلتونی H_λi+1 وH_λi بسیار کوچک است.
برای به­دست آوردن یک عبارت ریاضی برای روش تدریجی^[۵۲]، از معادله­ مربوط به روش اختلال ترمودینامیکی استفاده کرده، آن را به شکل یک سری تیلور می­نویسیم:[۱۷]
(۳-۳۲)
(۳-۳۳)
(۳-۳۴)
(۳-۳۵)
(۳-۳۶)
۳-۲-۳-خط سیر چند مرحله­ ای
خط سیر چند ­مرحله­ ای^[۵۳] روش جدیدی است که برای محاسبه­ی انرژی آزاد اتصال دی­ساکاریدها به­کار می­رود. این روش مبتنی بر ترکیباتی است که از چند خط سیر چند مرحله­ ای پیروی می­ کنند.
در این پژوهش، به بررسی اتصال هشت دی­ساکارید مختلف با گالکتین-۱ طحال گاو با بهره گرفتن از فرایند محاسباتی چند مرحله­ ای که توسط اشوریا^[۵۴] و آمتسل^[۵۵] معرفی شده، پرداخته شده است. در این فرایند، برای جدا کردن لیگاند از پاکت اتصال^[۵۶] و تخمین انرژی آزاد اتصال، این سیستم از یک حالت اتصالی^[۵۷] A به یک حالت غیر اتصالی B طی یک سری مراحل در طول مختصه­ی واکنش به جلو برده می­ شود. در طول این مراحل، این سیستم می ­تواند پیکربند­های غیر­تعادلی را نمونه گیری کند، اما در وقفه­­های مجزا سیستم اجازه دارد تا به تعادل برسد. از تبادل خط سیر بین این حالت­­های تعادلی به شیوه­­ای ترکیبی برای ایجاد یک توزیع کار- احتمال استفاده می­ شود. این توزیع، همراه با معادله­ زینسکی^[۵۸] برای محاسبه­ی انرژی آزاد فرایند به کار می­رود. انرژی­­های آزاد محاسبه­شده، انطباق خیلی خوبی با داده­­های آزمایشگاهی و تجربی برای اتصال دی­ساکاریدها دارد. خط سیر­های محاسبه­شده برای به­دست آوردن انرژی­های آزاد، جزئیات مولکولی فرایند غیر­اتصال را آشکار می­ کند که به درک ما از مکانیسم تمایل و اختصاصی بودن اتصال کربو­هیدرات کمک می­نماید[۵۴].
با داشتن یک مدل مینیمم­شده و تعادلی از قبل، خط سیر­های جداگانه با اختصاص سرعت­های تصادفی با توزیع ماکسول^[۵۹] در دمای K 298 به همه اتم­ها آغاز شد. این سیستم­ها قبل از آن که شبیه­سازی غیر­اتصالی^[۶۰] که در آن دی­ساکارید به­ طور مکانیکی از پروتئین به بیرون کشیده می­ شود را شروع کنند، برای مدت ps40 بیشتر در تعادل قرار گرفتند. برای اعمال نیروی کشنده، یک اتم مصنوعی^[۶۱] در مرکز جرم لیگاند تعریف شده است[۵۴]. سه­اتم دی­ساکارید به­ طور همزمان به اتم مصنوعی متصل شد، دو تا کربنی که در پیوند گلیکوزیدی مشارکت داشت و اکسیژن گلیکوزیدی. ثابت نیرو برای هر مهار، kcal.Å^-۲mol^-110 قرار داده شد (فاصله­ی اتم مصنوعی به اتم کربن: Å ۵/۱، فاصله­ی اتم مصنوعی به اتم اکسیژن Å ۲/۱). از اتم مصنوعی در طول شبیه­سازی برای کشیدن دی­ساکارید از جایگاه اتصال استفاده شد. برای جلوگیری از چرخش­ها و همسو­شدن دی­ساکارید در طول جدا شدن، محور کشش به عنوان محور دوم اینرسی مولکول دی­ساکارید تعریف شد. استفاده از بزرگترین محور اینرسی باعث ایجاد تصادف و برخورد بین پروتئین و لیگاند در طول فرایند کشیدن می­شد[۵۴].
هر کدام از حدود بین دی­ساکارید و اتم مصنوعی، به­عنوان یک محدوده­ هماهنگ در نظر گرفته شده است به طوری که نیرو به صورت معادله­ (۳-۳۷) نشان داده می­ شود:
(۳-۳۷)
که در آن موقعیت اتم متصل­شده به اتم مصنوعی و متناظر با موقعیت هدف در نقطه­ی تعادلی هر محدوده­ هماهنگ است (یعنی Å ۵/۱ از اتم قلابی برای کربن­ها و Å ۲/۱ برای اکسیژن). نیروی وارد بر دی­ساکارید با در نظر گرفتن تغییر شکل هر حدود محاسبه شد. ثابت نیروی مؤثر وقتی سهم هر فنر در راستای محور کشش^[۶۲] تصویر می­ شود، حدود kcal. ^-۲mol^-17 است. سهم کار هر کدام از این اتم­های متصل­شده به اتم مصنوعی متحرک در مرحله­ N به­ صورت زیر است:
(۳-۳۸)
که در آن، نیروی تصویر شده بر روی محور کشش و دیفرانسیل در جهت کشیدن است [۵۴].
قند از جایگاه اتصال پروتئین با بهره گرفتن از ترکیبی از انتقالات اتم مصنوعی همراه با تعادلات سیستم به نحوی که قبلاً شرح داده شد، جدا می­ شود. اتم مصنوعی در مراحل Å ۱ کشیده می­ شود. بعد از هر مرحله، سیستم اجازه می­یابد تا برای مدت ps5/2 به تعادل برسد. از این جا به بعد، به این مراحل، مراحل کوچک گفته می­ شود. وقتی ۱۰ مرحله­، کوچک که معادل پیش­روی به اندازه­ Å۱۰ است کامل شد سیستم برای مدت ps 50 به تعادل می­رسد. به این مرحله، «مرحله­ طولانی» گفته می­ شود. هر مجموعه از هشت مرحله­ طولانی، که معادل انتقال دی­ساکارید به میزان Å ۸ از جایگاه اتصال است، تشکیل یک مسیر را می­دهد.[۵۵].
برای هر کمپلکس گالکتین/دی­ساکارید، ۲۸ خط سیر محاسبه­شده تعداد خط سیر­ها بایستی به نحوی انتخاب می­شد تا تخمین­های انرژی آزاد برای هر یک از دی­ساکاریدها همگرا به خطای کمتر از kcal.mol^-1 ۳/۰ باشد. انرژی­های آزاد با روش ترکیب خط سیر چندمرحله­ای برآورد شد. این اولین باری است که از روش MSTC برای محاسبه­ی انرژی­های آزاد اتصال استفاده می­ شود و با توجه به نتایج جدول (۳-۳)، انرژی­های آزاد اتصال محاسبه شده توافق خیلی خوبی با اندازه ­گیری آزمایشگاهی و تجربی نشان می­دهد و این روش به اندازه­ کافی حساس است تا بین دی­ساکارید­های مختلف تفاوت قائل شود[۵۵].
جدول (۳- ۳) انرژی آزاد اتصال برای کمپلکس ­های گالکتین-۱/دی­ساکارید مختلف[۵۵]
۳-۲-۴- انتگرال­گیری ترمودینامیکی
به دلیل اینکه تمام کارهای این پایان نامه بر اساس روش انتگرال­گیری ترمودینامیکی انجام شده است درفصل بعد به طور مفصل در مورد این روش بحث شده است
۳-۳- کاربرد روش­های محاسبه­ی اختلاف انرژی آزاد
۳-۳-۱- چرخه­های ترمودینامیکی
بسیاری از پدیده ­های مورد توجه در مدل­سازی مولکولی شامل تعادل میان مولکول­هایی است که با هم برهمکنش کووالانسی دارند، معادله­ انرژی آزاد با ثابت تعادل به­ صورت است. به­عنوان نمونه، اتصال دو لیگاند (L₁ و L₂) را به یک مولکول گیرنده ® در نظر می­گیریم. این دو لیگاند می ­تواند دو بازدارنده­ی یک آنزیم باشد. چرخه­ی ترمودینامیکی برای دو فرایند اتصال در شکل (۳-۲) نشان داده شده است.

شکل (۳- ۲)چرخه ترمودینامیکی برای اتصال لیگاند­های L₁ و L₂ به گیرنده R. [17]
تمایل نسبی اتصال دو لیگاند برابر است که معمولاً به­ صورت ΔΔG نوشته می­ شود. به هر حال، باید بتوانیم به­وسیله­ شبیه­سازی فرایند واقعی اتصال، ΔG₁ و G₂Δ را محاسبه کنیم. برای این منظور باید لیگاند و گیرنده را از یک فاصله اولیه دور به­تدریج به هم نزدیک کنیم تا در نهایت کمپلکسی میان آن­ها تشکیل شود. اما معمولاً طی چنین فرآیندی تغییرات زیادی در ساختار گیرنده، لیگاند و حلال اطراف آن­ها صورت می­گیرد که نمونه­برداری کامل از فضای فاز با مشکل مواجه می­سازد[۱۷].
انرژی آزاد یک تابع حالت است و در نتیجه، مجموع تغییرات آن در یک چرخه­ی ترمودینامیکی، برابر با صفر است، بنابراین ، و ΔG₃ برابر با اختلاف انرژی آزاد و لیگاند در محلول، و ΔG₄ نیز برابر اختلاف انرژی آزاد دو کمپلکس بین­مولکولی است. اختلاف انرژی­های ΔG₃ و ΔG₄ با هیچ فرایند آزمایشگاهی قابل محاسبه نیستند. اما محاسبه­ی آن­ها با بهره گرفتن از رایانه نسبتاً ساده است. اختلاف انرژی آزاد فقط به نقاط انتهایی بستگی دارد و اگرچه این فرایند­ها غیر­فیزیکی­اند، محاسبه­ی اختلاف انرژی آزاد آن­ها به­وسیله­ شبیه­سازی معتبر­تر از دو فرایند فیزیکی متناظر با ΔG₁ و ΔG₂ است، مخصوصاً اگر لیگاند­های L₁ و L₂ ساختار مشابهی داشته باشند[۱۷].
۳-۳-۲- محاسبه­ی انرژی آزاد مطلق
در برخی موارد، می­توانیم چرخه­های ترمودینامیکی طراحی کنیم که ما را قادر به محاسبه­ی انرژی آزاد مطلق یک فرایند با بهره گرفتن از شبیه­سازی دینامیک مولکولی می­سازند. شکل (۳-۳) یک چرخه­ی ترمودینامیکی را برای اجتماع L و R و تشکیل یک کمپلکس LR در دو فاز گازی و محلول نشان می­دهد[۱۷].