نوع مقاله : مقاله پژوهشی

نویسندگان

1 گروه بیوفیزیک، دانشکده علوم زیستی، داشگاه تربیت مدرس، تهران، ایران

2 دانشکده ریاضی، دانشگاه تربیت مدرس، تهران، ایران

3 ایران، تهران، دانشگاه تربیت مدرس، دانشکده علوم زیستی، گروه بیوفیزیک ایران، تهران، مرکز تحقیقات علوم بنیادی، پژوهشکده علوم زیستی

چکیده

در میان انواع مختلف فازهای لیپیدی، ساختار لایه‌ای از اهمیت ویژه‌ای در درک فرآیندهای مرتبط با غشای سلولی اعم از همجوشی غشای لیپیدی برخوردار است. ازاین‌رو تحقیق در مدل‌سازی فاز L_α لیپیدی، بخش قابل‌توجهی از شبیه‌سازی ماکرومولکول‌های زیستی را به خود اختصاص داده است. در این تحقیق از درشت‌دانه‌سازی سامانمند به‌منظور مدل‌سازی فاز لایه‌ای مولکول لیپید dioleoylphosphatidylethanolamine (DOPE)-1،2 استفاده کرده‌ایم. مدل درشت‌دانه بدون حلال بوده و مولکول لیپید از ۱۴ دانه تشکیل می‌شود به‌طوری‌که انرژی پتانسیل سیستم شامل جملات پیوندی، غیر پیوندی و الکتروستاتیک هست. مدل‌سازی از روی شبیه‌سازی اتمی مولکول‌های لیپید در محتوای آبی %۴۴٫۶ وزنی-وزنی با استفاده از روش معکوس بولتزمن تکرارشونده و مونت‌کارلو معکوس انجام گرفت. اگرچه همگرایی مناسب در توابع توزیع پیوندی و زاویه‌ای در پایان فرآیند معکوس بولتزمن تکرارشونده حاصل گردید، اما همگرایی دقیق در توابع توزیع پیچشی و توابع توزیع شعاعی در پایان مونت‌کارلو معکوس بدست آمد. با استفاده از پتاسیل‌های حاصله، فاز لایه‌ای L_α در کمتر از 10 نانوثانیه بدست آمد. علاوه بر این، پتانسیل‌های درشت‌دانه قادر به مدل کردن همجوشی چندلایه‌ی لیپیدی و ساختار حد واسط ذرات لیپیدی هستند.

کلیدواژه‌ها

موضوعات

عنوان مقاله [English]

Simulation of the L_α lipid phase behavior using systematic coarse-graining of dioleoylphosphatidylethanolamine molecule

نویسندگان [English]

  • saeed mortezazadeh 1
  • yousef jamali 2
  • Hossein Naderimanesh 3

1 Department of Biophysics, Faculty of Biological Sciences, Tarbiat Modares University, Tehran, Iran

2 Department of Mathematics, Tarbiat Modares University, Tehran, Iran

3 Dept. of Nanobiotechnology/Biophysics, Faculty of Biological Sciences, Tarbiat Modares University, Tehran, I.R. of Iran Research Institute of Biological Sciences, Research Centre of Basic Sciences, Tehran, I.R. of Iran

چکیده [English]

Among the various types of lipid phases, the bilayer structure is of particular importance in understanding the processes associated with cell membranes such as lipid membrane fusion. Hence, the research in modeling of L_α lipid phase has been a significant part of the simulation of bio-macromolecules. In this study, we have exploited the systematic coarse-graining approaches for modeling of the 1,2- dioleoylphosphatidylethanolamine (DOPE) lipid molecule. The implicit solvent coarse-grained model of the lipid molecule is composed of 14 beads so that the potential energy of the system consists terms of bonded, non-bonded, and electrostatic. Constructing the Coarse-grained model was accomplished based on the atomistic simulation of lipid molecules with the water content of 45 %wt using the methods of iterative Boltzmann inversion and inverse Monte Carlo. Although proper convergence in the bond and angle distribution functions was achieved at the end of iterative Boltzmann inversion procedure, the exact convergence in the torsion angle distribution functions and the radial distribution functions was obtained at the end of the Monte Carlo process. Using the obtained results, the lamellar L_α lipid phase was obtained in less than 10 nanoseconds. In addition, the coarse-grained potentials are capable of modeling the lipid multilayer fusion and the intermediate structure of lipidic particles.

کلیدواژه‌ها [English]

  • molecular dynamics
  • coarse-graining
  • iterative boltzmann inversion
  • inverse monte carlo

شبیه‌سازی رفتار فاز لیپیدی  با استفاده از درشت‌دانه‌سازی سامان مند مولکول dioleoylphosphatidylethanolamine

سعید مرتضی زاده1، یوسف جمالی2 و حسین نادری منش1،3*

1 ایران، تهران، دانشگاه تربیت مدرس، دانشکده علوم زیستی، گروه بیوفیزیک

2 ایران، تهران، دانشگاه تربیت مدرس، دانشکده علوم ریاضی

3 ایران، تهران، مرکز تحقیقات علوم بنیادی، پژوهشکده علوم زیستی

تاریخ دریافت: 12/9/97                تاریخ پذیرش: 2/10/97

چکیده

در میان انواع مختلف فازهای لیپیدی، ساختار لایه‌ای از اهمیت ویژه‌ای در درک فرآیندهای مرتبط با غشای سلولی اعم از همجوشی غشای لیپیدی برخوردار است. ازاین‌رو تحقیق در مدل‌سازی فاز لیپیدی  ، بخش قابل‌توجهی از شبیه‌سازی ماکرومولکولهای زیستی را به خود اختصاص داده است. در این تحقیق از درشت‌دانه­سازی سامان مند به‌منظور مدل‌سازی فاز لایه‌ای مولکول لیپید dioleoylphosphatidylethanolamine (DOPE)-1،2 استفاده شده است. مدل درشت‌دانه بدون حلال بوده و مولکول لیپید از ۱۴ دانه تشکیل می‌شود به‌طوری‌که انرژی پتانسیل سیستم شامل جملات پیوندی، غیر پیوندی و الکتروستاتیک هست. مدل‌سازی از روی شبیه‌سازی اتمی مولکولهای لیپید در محتوای آبی 6/44 درصد وزنی- وزنی با استفاده از روش معکوس بولتزمن تکرارشونده و مونت‌کارلو معکوس انجام گرفت. اگرچه همگرایی مناسب در توابع توزیع پیوندی و زاویه‌ای در پایان فرآیند معکوس بولتزمن تکرارشونده حاصل گردید، اما همگرایی دقیق در توابع توزیع پیچشی و توابع توزیع شعاعی در پایان مونت‌کارلو معکوس به دست آمد. با استفاده از پتاسیل‌های حاصله، فاز لایه‌ای  در کمتر از 10 نانوثانیه به دست آمد. علاوه بر این، پتانسیلهای درشت‌دانه قادر به مدل کردن همجوشی چندلایه لیپیدی و ساختار حد واسط ذرات لیپیدی هستند.

واژه‌هایکلیدی:دینامیک مولکولی، درشت‌دانه‌سازی، معکوس بولتزمن تکرارشونده، مونت‌کارلو معکوس

* نویسنده مسئول، تلفن ۸۲۸۸۴۴۱۰-۲۱ ، پست الکترونیکی: [email protected]

مقدمه

 

فسفولیپیدها مواد اصلی غشاهای سلولی هستند که سلولها را از محیط اطراف خود جدا می‌کنند. علاوه بر نقش کلیدی در زیست‌شناسی، فسفولیپیدها خواص زیست‌سازگاری بالایی دارند که به‌عنوان نانوحامل در تحویل دارو مورد استفاده قرار می‌گیرند(10). زیست‌سازگاری بالا و تحویل قابل‌ کنترل داروها، فسفولیپیدها را به‌عنوان یک سامانه تحویل دارو در کاربردهای پزشکی تبدیل کرده است(2 و 19). این قبیل توانایی فسفولیپیدها برگرفته از خاصیت چند­ریختی در ایجاد طیف گسترده‌ای از ساختارهای خودآرا ازجمله ساختارهای لایه‌ای، ساختارهای کروی مانند میسل و لیپوزوم، ساختارهای لوله‌ای مانند فاز شش ضلعی و برخی دیگر از ساختارهای پیچیده مانند فاز مکعبی لیپیدی می‌باشند. مدل مایع موزاییک به‌عنوان یک پیش‌فرض اصلی برای درک ساختار و عملکرد غشاء با ارجاع به ساختار دولایه مایع­بلور استفاده می‌شود.

دو عامل مهم در شکل‌گیری ساختارهای خودآرای لیپیدی نقش دارند. اولین عامل در ایجاد شکل کروی یا مخروطی لیپید، رابطه بین اندازه (عرض) سر و دم لیپید است که خود بستگی به ماهیت شیمیایی سر لیپید و طول زنجیره­ها و همچنین تعداد و محل قرارگیری هیدروکربنهای غیراشباع دارد که موجب القای یک عامل بسته‌بندی خاص برای لیپیدها می‌شود. ازاین‌رو، دوگانه‌دوستی فسفولیپیدها موجب تجمع لیپیدی ناشی از برهمکنشهای آب‌گریز در محلول آبی می‌شود. عامل دوم ناشی از پارامترهای ترمودینامیکی مانند دما و میزان محتوای آب نمونه است. در این راستا، مدل‌سازی مولکولی ساختارهای لیپیدی، اطلاعات ارزشمندی در مورد مکانیسم خودآرایی آنها ارائه می‌دهد که یکی از چالشهای اساسی در تحقیقات مواد نرم محسوب می‌شود. ازاین‌رو شبیه‌سازی دینامیک مولکولی به‌عنوان ابزاری مناسب در تحقیقات ماکرومولکولهای زیستی محسوب می‌گردد(1).

مقیاسهای زمانی و فضایی پدیده‌های خودآرایی بسیار بزرگ‌تر از آن است که بتوان از طریق دینامیک مولکولی اتمی ردیابی کرد. برای رفع این محدودیت، یک رهیافت، کاهش تعداد درجات آزادی سیستم با حذف جزییات غیر ضرور به‌منظور سرعت بخشیدن به محاسبات است. این روش، درشت‌دانه‌سازی نامیده می‌شود که به دو رویکرد متفاوت بالا به پایین و پایین به بالا تقسیم‌بندی می‌گردد(30). رویکردهای پایین به بالا که به درشت‌دانه‌سازی سامان مند مشهور است از نتایج شبیه‌سازی اتمی با به‌کارگیری روشهای مختلفی مانند معکوس بولتزمن تکرارشونده(29)، مونت‌کارلو معکوس(21)، تطبیق نیرو(15) و آنتروپی نسبی(33) برای ساخت مدل درشت‌دانه استفاده می‌کند. با توجه به بازتولید توابع توزیع برگرفته از شبیه‌سازی اتمی، روشهای فوق مبتنی بر ساختار می‌باشند به‌جز روش تطبیق نیرو که به جای توابع توزیع، برآیند نیرو‌های اتمی را بازتولید می‌کند. مرور خوبی بر روی روشهای مختلف درشت‌دانه‌سازی سامان مند در مقاله برینی و همکارانش قابل یافت است(4). در طرف دیگر، هدف اصلی رویکردهای بالا به پایین، تولید داده‌های کلیدی تجربی مانند انرژی آزاد مولکول از حلال قطبی به غیر قطبی است. مدل مارتینی(23) معروف‌ترین آنهاست که از انرژی پتانسیلهای لنارد-جونز و الکترواستاتیک برای برهمکنشهای غیر پیوندی و انرژی پتانسیل فنر برای برهمکنشهای پیوندی و زاویه‌ای استفاده می‌کند. هر دو این رویکردها دارای محدودیتهای خاص خود هستند به‌طوری‌که روشهای پایین به بالا می‌تواند جزئیات بیشتری از برهمکنشهای بین ذرات را به‌عنوان پتانسیل نیروی متوسط به دست ‌آورند، درحالی که رویکردهای بالا به پایین، یک چارچوب میدان‌نیرویی مناسب جهت توسعه به سیستمها و مولکولهای دیگر را فراهم می‌کند(14).

در میان ساختارهای مختلف فسفولیپیدها، فاز لایه‌ای  از اهمیت خاصی در بررسی عملکرد و ساختار غشای سلولی برخوردار است. هر دو رویکرد بالا به پایین و پایین به بالا برای مطالعه شکل‌گیری انواع مختلفی از ساختارهای خودآرای لیپیدی به‌ویژه ساختارهای دولایه استفاده‌شده است. از میان روشهای بالا به پایین می‌توان به مدل مارتینی توسط مارینک و مارک(22) و مدل اِلبا توسط اُرسی و اِسِکس(26) اشاره کرد. در هر دو مدل از مولکولهای آب در شبیه‌سازی درشت‌دانه استفاده می‌شود به‌طوری‌که در مدل مارتینی هر ذره حلال معادل ۴ مولکول آب اتمی است درحالی که در مدل اِلبا هر مولکول آب معادل یک مولکول آب اتمی به انضمام یک دوقطبی الکتریکی است. از روشهای مختلف درشت‌دانه‌سازی سامان مند نیز در بررسی شکل‌گیری ساختار لایه‌ای لیپیدها استفاده ‌شده است. شِلی و همکاران(34) از روش معکوس بولتزمن تکرارشونده برای به دست آوردن پتانسیلهای درشت‌دانه به همراه حلال آبی نگاشت شده از هر مولکول آب اتمی استفاده کردند. وُث و همکارانش(16 و 20) مدل درشت‌دانه بدون حلالی را با استفاده از روش تطبیق نیرو ارائه کرده‌اند. در تعدادی از مدلهای بدون حلال، چالش برهمکنشهای آب‌گریز در فرآیند خودآرایی را با معرفی یک برهمکنش اضافی حل کرده‌اند. در همین راستا، وَنگ و دِسِرنو(38)، سُدت و هید-گوردُن(36) از برهمکنشهای بلند بُرد جاذبه­ای برای تقلید اثر آب­گریزی دمهای لیپیدی به‌منظور ایجاد و تثبیت ساختار دولایه بهره بردند.

در تحقیق پیشرو، شرایط لازم برای تولید مدل درشت‌دانه بدون حلال با استفاده از روش درشت‌دانه‌سازی سامان مند با استفاده از روشهای معکوس بولتزمن تکرار­شونده و مونت­کارلو معکوس به­منظور حصول فاز  مولکول DOPE مورد بررسی قرار گرفته است.

مواد و روشها

مبانی نظری: هدف از درشت‌دانه‌سازی سامان مند، دستیابی به پتانسیلهای درشت‌دانه براساس طرح نگاشت مدل اتمی است. پتانسیلهای درشت‌دانه معمولاً شامل جملات پیوندی، غیر پیوندی و کولنی می­باشند. حدس اولیه برای جملات پیوندی و غیر پیوندی را می­توان از طریق معکوس بولتزمن مستقیم به دست آورد(29).

 

 

 

(1)                               

جایی که  و  و به ترتیب برابر ثابت بولتزمن و دما می‌باشند. میانگین تابع توزیع فاصله ، زاویه  و زاویه پیچشی  از روی مسیر شبیه­سازی اتمی نگاشت‌شده محاسبه می­شود. مخرج آرگومان توابع لگاریتم، تبدیل ژاکوبی مختصات کروی به‌منظور نرمال­سازی فضایی توابع توزیع است. بااین‌حال، پتانسیل اولیه معمولاً نمی­تواند توابع توزیع مرجع اتمی، مخصوصاً برهمکنشهای غیرپیوندی را بازتولید کند. ابتدایی­ترین روش برای بهبود پتانسیلهای درشت‌دانه، معکوس بولتزمن تکرارشونده است (31) که در فرآیندی از شبیه­سازیهای متوالی، اصلاح پتانسیل انجام می­گیرد. پتانسیلهای درشت‌دانه به صورت جدولی از اعداد می­باشند که هر یک از مقادیر این جداول در طی تکرار اُم بر اساس رابطه زیر اصلاح می­گردد:

(2)               

جایی که ۱  ۰، عامل تنظیم کننده همگرایی روش تکرارشونده است. هنگامی که مقدار اُم از تابع توزیع  از شبیه‌سازی اُم به سمت هیستوگرام مرجع  نزدیک شود،  اصلاح پتانسیل درشت‌دانه به سمت صفر میل می­کند.

فرضیه اصلی در روش معکوس بولتزمن تکرارشونده این است که تمامی مقادیر جداول پتانسیل درشت‌دانه مستقل از یکدیگرند به­طوری­که همبستگی متقابل میان برهمکنشها در اصلاح پتانسیل درنظر گرفته نمی‌شود. این محدودیت منجر به مشکل همگرایی در سیستمهای چندجزئی شامل انواع مختلفی از برهمکنشهای غیر پیوندی می‌گردد. برای حل این مشکل، از روش تکرارشونده دیگری موسوم به مونت­کارلو معکوس استفاده می­شود که از ماتریس همبستگی متقابل بین برهمکنشهای مختلف برای اصلاح پتانسیلهای درشت‌دانه استفاده می­کند. جزئیات روش مونت­کارلو معکوس در مقاله لوبارتسف (21 و 25) آمده است. این روش یک مسئله حل معکوس ماتریس است، به­طوری­که پتانسیلهای درشت‌دانه بر اساس رابطه زیر اصلاح می­شود:

 

(3)       

در مقاله اصلی روش مونت­کارلو معکوس، ماتریس همبستگی متقابل  از روی نمونه­برداری مونت­کارلو محاسبه شده است، با این ‌وجود هر الگوریتم نمونه­برداری دیگری مانند شبیه­سازی دینامیک­ مولکولی، مادامی‌که توزیع کانونی فضای فاز را تولید کند، قابل استفاده است. با این حال، برای جلوگیری از وقوع خطای آماری، شبیه­سازی طولانی مدت در هر تکرار برای محاسبه ماتریس همبستگی مورد نیاز است. ماتریس همبستگی را می­توان بین تمام جملات برهمکنشهای پیوندی و غیر پیوندی محاسبه کرد.

شبیه‌سازی اتمی: شبیه­سازی اتمی مولکولهای DOPE در آب با استفاده از نرم افزار Gromacs 5.0.7 (3) انجام شد. با توجه به مقاله (25)، ۶۰ مولکول DOPE در ۲۰۰۰ مولکول آب متناسب با محتوای آبی 6/44 درصد وزنی- وزنی به مدت ۵۰۰ نانوثانیه با شروع از ساختار تصادفی مولکولهای لیپید شبیه­سازی گردید. از آنجا که بار کل هر مولکول لیپید صفر است، نیازی به اضافه کردن یون برای خنثی کردن جعبه نیست. میدان نیروی Slipids (8 و 17) برای مولکول DOPE و مدل tip3p (18) برای مولکولهای آب استفاده گردید. شبیه‌سازی در جعبه مکعبی با گام زمانی ۲ فمتوثانیه در هنگرد NPT، در دمای ۳۰۳ کلوین و فشار ۱ اتمسفر انجام شد. روش V-rescale (7) با جفت­شدگی زمانی 5/0 پیکوثانیه و روش پارینلو-رهمان (27) با جفت­شدگی زمانی ۵ پیکوثانیه به ترتیب برای کنترل دما و فشار به‌کار گرفته ­شد. به‌منظور کنترل ارتعاشات پیوندهای متصل به اتم هیدروژن، از قیدگذاری Lincs (11) استفاده گردید. برای محاسبه پتانسیل الکترواستاتیک، روش ذره-شبکه­ای ایوالد (9) با شعاع قطع 2/1 نانومتر استفاده شد. همچنین، شعاع قطع درونی و بیرونی در محاسبه پتانسیل لنارد-جونز به ترتیب برابر با 8/0 و 2/1 نانومتر انتخاب گردید(37). قبل از شروع شبیه‌سازی، کمینه‌سازی انرژی با استفاده از روش تندترین­کاهش نیز انجام پذیرفت. مسیر شبیه­سازی اتمی بعد از ۱۰۰ نانوثانیه تعادل‌رسانی به‌منظور محاسبه توابع توزیع مرجع مورد بهره­برداری قرار گرفت. نرم افزار VMD.1.9.2 (13) و packmol (24) به­منظور نمایش و تولید ساختار تصادفی اولیه مورد استفاده قرارگرفتند.

طرح نگاشت: اولین گام در ساخت یک مدل درشت­دانه، نحوه بیان سیستم بر اساس یک طرح نگاشت است. در این مطالعه، مولکولهای آب از شبیه­سازی درشت­دانه حذف شدند در حالی که هر مولکول DOPE که شامل ۱۲۹ اتم است به یک مولکول درشت­دانه با ۱۴ دانه نگاشت گردید که طرح نگاشت آن در شکل ۱ نشان داده شده است. مدل درشت­دانه از شش نوع دانه شامل N، P، CO، C3، CdB و C4 تشکیل شده که هر کدام از یک گروه شیمیایی مجزا نگاشت شده است. بار جزئی هر دانه از جمع اتمهای درگیر در آن محاسبه می‌شود. در مدل درشت­دانه مولکول DOPE شش نوع برهمکنش پیوندی شامل N-P، P-CO، CO-C3، C3-C3، C3-CdBو C3-C4 تعریف گردید. علاوه بر این، ۷ نوع پتانسیل زاویه‌ای شامل NP-CO، P-CO-C3، CO-P-CO، CO-C3-C3، C3-C3-CdB، C3-CdB-C3 و CdB-C3-C4، و ۶ نوع پتانسیل زاویه پیچشی شامل NP-CO-C3، P-CO-C3-C3، CO-P-CO-C3، CO-C3-C3-CdB، C3-C3-CdB-C3 و C3-CdB-C3-C4 بر اساس انواع برهمکنشهای پیوندی استخراج شد. علاوه بر برهمکنشهای پیوندی، ۲۱ نوع برهمکنش غیر پیوندی بین ۶ نوع مختلف دانه وجود دارد. برهمکنشهای غیرپیوندی برای جفت‌دانه‌هایی که درگیر برهمکنشهای پیوندی، زاویه‌ای و یا پیچشی نباشند محاسبه می‌شود.

 

 

شکل ۱- طرح نگاشت مولکول DOPE. نگاشت گروههای مختلف شیمیایی با رنگ متمایز نشان داده شده است. N (گروه آمین)، P (گروه فسفات)، CO (گروه استر)، C3 (گروه پروپان)، CdB (گروه پروپِن)، و C4 (گروه بوتان) را نشان می‌دهد.

 

در ادامه، طرح نگاشت برای تبدیل مسیر اتمی به مسیر درشت­دانه به‌منظور محاسبه توابع توزیع مرجع مورد استفاده قرار گرفت.

شبیه‌سازی درشت‌دانه: شبیه‌سازی سیستمهای درشت‌دانه با استفاده از نرم افزار LAMMPS-2017 (28) در هنگرد NVT در دمای ۳۰۳ کلوین انجام شد. LAMMPS مجهز به بسته نرم افزاری GPU (5 و 6) هست که اجازه شبیه‌سازی موازی بر روی انواع سخت افزارهای مرکب از CPU و GPU را فراهم می‌آورد که موجب افزایش بازدهی شبیه‌سازی سیستمهای گوناگون مخصوصاً مواردی که از پتانسیلهای جدول‌بندی شده تشکیل شده‌اند، می‌گردد. ازآنجایی‌که شبیه‌سازی درشت‌دانه بدون حلال انجام گرفته است، تنظیم سیستم با محتوای آبی متفاوت از روی تغییر چگالی حجمی مولکولهای لیپید انجام گرفت. در همین راستا، از دینامیک لانژوین با جفت‌شدگی زمانی 5/0 پیکوثانیه و گام زمانی ۱۰ فمتوثانیه برای نمونه‌برداری فضای فاز استفاده شد. قبل از هر شبیه‌سازی از روش گرادیانهای مزدوج برای کمینه‌سازی انرژی پتانسیل سیستم استفاده شده است. شعاع قطع در پتانسیلهای غیر پیوندی برابر ۲ نانومتر تنظیم شد به‌طوری‌که برای صفر کردن انرژی پتانسیل و نیرو در شعاع قطع از یک تابع درجه ۳ در فاصله 8/1 تا ۲ نانومتر استفاده گردید.

محاسبه پتانسیل درشت‌دانه با استفاده از روشهای تکرارشونده: پتانسیلهای پیوندی و غیر پیوندی درشت‌دانه در طی فرآیندهای معکوس بولتزمن تکرارشونده و مونت‌کارلو معکوس به دست می‌آیند. روش معکوس بولتزمن تکرارشونده برای شروع مورد استفاده قرار گرفت، به‌طوری‌که پتانسیل اولیه برای راه اندازی شبیه‌سازی درشت‌دانه از معکوس بولتزمن مستقیم توابع توزیع مرجع تخمین زده شد. پس از تکمیل ۲۰ مرحله از فرآیند معکوس بولتزمن تکرارشونده، ۲۰ مرحله فرآیند مونت‌کارلو معکوس برای به دست آوردن بهترین پتانسیل درشت‌دانه انجام شد. زمان شبیه‌سازی هر تکرار در فرآیند معکوس بولتزمن تکرارشونده ۲۷ نانوثانیه با فاصله زمانی نمونه‌برداری 5/2 پیکوثانیه بود که ۲ نانوثانیه ابتدایی هر مسیر به عنوان تعادل‌رسانی سیستم حذف گردید. زمان شبیه‌سازی و تعادل‌رسانی در ۱۰ مرحله اول فرآیند مونت‌کارلو معکوس به ترتیب برابر ۵۴ و ۴ نانوثانیه تعیین شد، درحالی که برای ۱۰ مرحله انتهایی به ۱۰۸ و ۸ نانوثانیه رسید. دو برابر شدن زمان شبیه‌سازی و تعادل‌رسانی، این اطمینان را می‌دهد تا از محاسبه ماتریس همبستگی خطادار ناشی از نمونه‌برداری ناکافی فضای فاز جلوگیری شود. پتانسیل نهایی اصلاح شده درشت‌دانه در پایان فرآیند مونت‌کارلو معکوس برای مطالعات بیشتر مورد استفاده قرار گرفت. از آنجایی که نرم‌افزار‌های موجود برای به دست‌آوردن پتانسیلهای درشت‌دانه شامل نواقصی هستند، برنامه‌ای به زبان پایتون نوشته شد که دربرگیرنده تمام مراحل درشت‌دانه‌سازی سامان مند ازجمله: تبدیل مسیر اتمی به درشت‌دانه بر اساس ماتریس نگاشت، تولید فایل ساختاری درشت‌دانه، محاسبه توابع توزیع و اصلاح پتانسیلهای درشت‌دانه طی فرآیندهای تکرارشونده هست. نرم‌افزار‌های مشابه دیگر مثل ‌MagiC (25)، پتانسیل پیچشی را پوشش نمی‌دهد. از طرفی نرم‌افزار VOTCA (32)، از تابع درجه دوم و نمایی برای برون‌یابی پتانسیلهای پیوندی استفاده می‌کند. در هر حال استفاده از این توابع موجب می‌شود پیچش کامل حول پیوند رخ ندهد و زوایای پیچشی به محدوده 180- و 180+ درجه نزدیک نشوند. در کد نوشته شده، پتانسیلهای پیچشی به صورت دوره‌ای دورن و برون‌یابی می‌شوند به‌طوری‌که نقیصه فوق برطرف می‌گردد. از طرفی، پتانسیلهای غیر پیوندی در نرم‌افزار MagiC با استفاده از تابع نمایی به سرعت در شعاع قطع پتانسیل به صفر میل می‌کند. در نرم افزار VOTCA هر پتانسیل غیرپیوندی به‌طور مستقل به سمت صفر جا به جا می‌شود که اثرات دافعه و جاذبه‌ای در مرز شعاع قطع از بین می‌رود. در کد نوشته شده با تقلید از روش تعویض تابع (37) در محاسبه پتانسیلهای لنارد- جونز در شبیه‌سازی اتمی، از یک تابع درجه 5 برای صفر شدن پتانسیل و نیرو در شعاع قطع دوم استفاده می‌شود به‌طوری‌که در شعاع قطع اول، هم پتانسیل و هم نیرو پیوسته خواهند بود. شعاع قطع اول و دوم در پتانسیلهای غیرپیوندی به‌ترتیب برابر 8/1 و 0/2 نانومتر است.

نتایج و بحث

نتایج شبیه‌سازی اتمی مولکول DOPE در شکل ۲ نشان داده شده است. از ۴۰۰ نانوثانیه پایانی شبیه‌سازی اتمی، ساختار دولایه لیپیدی با چگالی ۳-(nm) 45/0 به دست آمد.

 

شکل ۲- ساختار به دست‌آمده از شبیه‌سازی اتمی. در تصویر اتمهای هیدروژن نشان داده نشده است.

توسعه پتانسیل درشت‌دانه از طریق ۴۰ مرحله تکرار بر اساس مسیر اتمی حاصله دنبال گردید که همگرایی بین توابع توزیع هدف و آنهایی که از شبیه‌سازی درشت‌دانه حاصل می‌شود، توسط مجموع میانگین مربع انحراف توابع توزیع در هر مرحله محاسبه می‌شود.

(4)  

در رابطه (۴)،  برابر تعداد کل برهمکنشهای پیوندی و غیر پیوندی و  تعداد مقادیر هر یک از جداول پتانسیل است. انحراف ( ) در ۲۰ مرحله ابتدایی معکوس بولتزمن تکرارشونده و ۲۰ مرحله نهایی مونت‌کارلو معکوس در شکل ۳ نشان داده شده است. شیب کاهش انحراف در ۲۰ تکرار پایانی بیشتر از ۲۰ تکرار ابتدایی بوده که نشان‌دهنده مزیت استفاده از روش مونت‌کارلو معکوس در برابر معکوس بولتزمن تکرارشونده در سیستمهای چندجزیی است. با این ‌حال، استفاده از روش معکوس بولتزمن تکرارشونده برای آغاز فرآیند تصحیح پتانسیل درشت‌دانه لازم است. در ۱۰ تکرار انتهایی به دلیل شروع شبیه‌سازی از ساختار تصادفی، افزایش نسبی در انحراف تکرارهای ۴۱ تا ۴۳ مشاهده می‌شود.

 

شکل ۳- میزان انحراف از توابع توزیع هدف در طول فرآیند تکرارشونده.

شکل ۴ (الف و ب) همگرایی در تابع توزیع و پتانسیل برهمکنش زاویه پیچشی CO-P-CO-C3 را در طول تکرار نشان می‌دهد. این پتانسیل پیچشی، چرخش حول پیوند P-CO در وسط زنجیره لیپیدی را کنترل می‌کند. همان‌طور که مشاهده می‌شود، تغییرات قابل توجهی در انرژی پتانسیل به‌منظور بازتولید تابع توزیع مرجع رخ داده است. استفاده از روش مونت‌کارلو معکوس منجر به اصلاح دقیق پتانسیلهای زاویه پیچشی و همگرای کامل تابع توزیع درشت‌دانه و مرجع می‌شود. در بقیه برهمکنشهای زوایای پیچشی نیز بخش اصلی تصحیح پتانسیل در طی فرآیند معکوس بولتزمن تکرارشونده رخ داده است. به رغم نیاز به اصلاح پتانسیل در تمامی برهمکنشهای زاویه پیچشی، اصلاح پتانسیلهای پیوندی و زاویه‌ای ناچیز بوده است. این موضوع نشان‌دهنده استحکام برهمکنشهای پیوندی و زاویه‌ای است تاجایی که استفاده از معکوس بولتزمن مستقیم برای تولید این پتانسیلها کافی به‌نظر می‌رسد.

 

 

شکل ۴- الف) همگرایی در تابع توزیع و ب) بهبود پتانسیل برهمکنش پیچشی CO-P-CO-C3 در طی فرآیند تکرارشونده.

 

شکل ۵ (الف و ب) نشان‌دهنده همگرایی در تابع توزیع شعاعی و بهبود در انرژی پتانسیل برای برهمکنش غیر پیوندی میان دانه‌های C4 و P را در طی فرآیند تکرارشونده است. در طی مراحل معکوس بولتزمن تکرارشونده ، علاوه بر بازتولید قله اصلی تابع توزیع شعاعی، دو قله اضافی در ناحیه 9/0 تا 1/1 نانومتر شکل گرفته است. قله‌های اضافی ناشی از این است که فرض اولیه در روش معکوس بولتزمن تکرارشونده، مستقل بودن برهمکنشها از یکدیگر بوده به‌طوری‌که پتانسیلهای درشت‌دانه بدون درنظرگرفتن همبستگی میانشان اصلاح می‌شوند. در واقع، به‌روز رسانی هر یک از مقادیر در جداول انرژی پتانسیل، نه تنها تابع توزیع متناظر با همان فاصله شعاعی را تغییر می‌دهد بلکه بر سایر مقادیر نیز تأثیر می‌گذارد. ازاین‌رو، دو چاه انرژی پتانسیل اضافی در ناحیه 9/0 تا 1/1 نانومتر شکل گرفته است. این نقص، در برخی دیگر از پتانسیلهای غیر پیوندی شامل C4-N،  C4-P، CdB-Nو CdB-P در پایان روند معکوس بولتزمن تکرارشونده دیده شد. ازآنجایی‌که این پتانسیلهای غیر پیوندی مربوط به تعامل سر لیپید با مرکز و دم مولکول لیپید است، توجه به همبستگی متقابل میان آنها در فرآیند اصلاح پتانسیل اجتناب‌ناپذیر است. در طی ۲۰ مرحله از روند مونت‌کارلو معکوس، چاههای پتانسیل اضافی به طور کامل محو شده و پتانسیل صافی تولید شده است. اگرچه استفاده از روش مونت‌کارلو معکوس برای برهمکنشهای پیوندی منجر به اصلاح ظریف در انرژی پتانسیل می‌شود، اما استفاده از آن در برهمکنشهای غیر پیوندی ضروری است.

 

 

شکل ۵- الف) همگرایی در تابع توزیع شعاعی و ب) اصلاح پتانسیل برهمکنش غیر پیوندی C4-P در طی فرآیند تکرارشونده.

 

به‌منظور بررسی پتانسیلهای حاصله در بازتولید ساختار دولایه تشکیل شده در شبیه‌سازی اتمی، نتایج شیبه‌سازی ۶۰ مولکول درشت‌دانه لیپید در چگالی متناسب با شبیه‌سازی اتمی در شکل ۶ آمده است. همان‌طورکه مشاهده می‌شود ساختار دولایه لیپیدی به‌درستی در ۱۰ نانوثانیه شبیه‌سازی درشت‌دانه حاصل گشته است. علت کوتاه‌بودن نتایج حاصل از شبیه‌سازی درشت‌دانه حذف حلال است. درعین‌حال، کاهش تعداد درجات آزادی سیستم در کنار اینکه پتانسیلهای درشت‌دانه دراصل، پتانسیل نیرو‌میانگین هستند، باعث کوتاه شدن زمان شبیه‌سازی فرآیندهای خودآرایی می‌شود. از این‌رو زمان شبیه‌سازیهای درشت‌دانه قابل تطابق با زمان حقیقی فرآیندهای زیستی نیست.

 

شکل ۶- شکل دولایه لیپیدی حاصل از ۶۰ مولکول DOPE.

نتایج شکل ۶ نشان‌دهنده پایداری چند لایه لیپیدی نیست بلکه، به دلیل اینکه با روشن کردن تصاویر دوره ای به‌منظور نمایش هر چه بهتر تشکیل دولایه لیپیدی، ساختار چند لایه نیز القاء شده است. در واقع پتانسیلهای حاصله قادر به مدل‌کردن دولایه لیپیدی هستند.

به‌منظور بررسی فرآیند همجوشی دولایه لیپیدی تعداد ۴۸۰ مولکول DOPE در محتوای آبی ۴۵ درصد وزنی- وزنی با شروع از یک ساختار چندلایه به مدت ۴۰ نانوثانیه شبیه‌سازی شد. نتایج ساختارهای شکل‌گرفته در طول زمان در شکل ۷ نشان داده شده است. در حین نزدیک شدن صفحات لیپیدی به یکدیگر برای همجوشی به یک دولایه پایدار، تراکم موضعی لیپیدها تغییر می‌کند که در محدوده زمانی 5/8 نانوثانیه قابل مشاهد است. در ادامه، با تکمیل همجوشی غشایی، استوانه‌هایی توخالی موسوم به ذرات لیپیدی تشکیل می‌شوند. در واقع آن دسته از لیپیدهایی که به‌دلیل تراکم لیپیدی، فرصت رفتن به یکی از دو وجه ساختار دولایه را پیدا نکرده‌اند، یک استوانه توخالی در داخل غشاء تشکیل می‌دهند. در واقع ذرات لیپیدی ساختارهای موقتی هستند که در طی فرآیند همجوشی دولایه‌های لیپیدی تشکیل می‌شوند(35).

 

شکل ۷- همجوشی چند لایه لیپیدی ۴۸۰ مولکول DOPE در طول زمان.

نتیجه گیری

در این تحقیق از رویکرد درشت‌دانه‌سازی سامان مند برای تولید یک مدل درشت‌دانه بدون حلال برای مولکول DOPE با ۱۴ دانه به‌ازای هر لیپید استفاده شده است. مدل درشت‌دانه بر اساس شبیه‌سازی اتمی یک ساختار دولایه لیپیدی در محتوای آبی ۴۵ درصد وزنی- وزنی انجام گرفت. اگر چه از روش معکوس بولتزمن تکرارشونده برای اصلاح پتانسیلهای درشت‌دانه استفاده گردید، اما پتانسیلهای درشت‌دانه نهایی با استفاده از روش مونت‌کارلو معکوس با دقت بالایی تولید گردیدند. پتانسیلهای حاصله نه تنها قادر به تشکیل و پایدار سازی دولایه لیپیدی بودند، بلکه همجوشی غشاهای لیپیدی و تشکیل ساختار گذار ذرات لیپیدی نیز، به‌درستی مدل شد. در مقالات وُث و دیگران (16 و 20) سعی بر استخراج پتانسیلهای درشت‌دانه برای لیپیدهای ‌DOPE و DOPC با استفاده از روش تطبیق نیروها شده است. همچنین در مقاله هیل و مَک‌گلینچی (12) نیز پتانسیلهای تعمیم یافته‌ای برای مولکولهای لیپید ‌DOPE، POPC و POPG با استفاده از روش تطبیق نیرو ارائه شده است. نتایج همگرایی در توابع توزیع شعاعی و پیوندی در این کار به مراتب بهتر از گزارشات قبلی است. از طرفی شکل‌گیری همجوشی غشایی با استفاده از پتانسیلهای فعلی صورت پذیرفته است درحالی که در گزارشات قبلی به خواص دولایه غشایی بیشتر تاکید شده و فرآیند همجوشی بررسی نشده است. از این‌رو پتانسیلهای حاصله برای بررسی رفتار فاز  مولکول DOPE و فرآیندهای مرتبط با این فاز لیپیدی مناسب می‌باشند.

تشکر

این تحقیق بخشی از پایان‌نامه سعید مرتضی‌زاده است. نویسندگان تمایل دارند که از حمایتهای دانشگاه تربیت مدرس و مرکز پژوهشهای بنیادی تشکر و قدردانی نمایند.

1- صفرزاده et al.,2015, "بررسی نحوه تأثیر جهش ها بر غیر فعال شدن آنزیم پیرازین آمیداز با شبیه سازی دینامیک مولکولی," مجله پژوهش‌های سلولی و مولکولی, vol. 28, no. 2, pp. 266-278.
2- صدیفیان, مارنانی, ر., and منش, ر.,2018, "بررسی نفوذ داروهای آسپیرین و ایبوپروفن در غشاء دولایه لیپیدی به کمک شبیه سازی دینامیک مولکولی," مجله پژوهش‌های سلولی و مولکولی, vol. 31, no. 2, pp. 317-328.
 
3- Abraham, M. J. et al.,2015, "GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers," SoftwareX, vol. 1, pp. 19-25.
4- Brini, E., Algaer, E. A., Ganguly, P., Li, C., Rodríguez-Ropero, F., and van der Vegt, N. F.,2013, "Systematic coarse-graining methods for soft matter simulations–a review," Soft Matter, vol. 9, no. 7, pp. 2108-2119.
5- Brown, W. M., Wang, P., Plimpton, S. J., and Tharrington, A. N.,2011, "Implementing molecular dynamics on hybrid high performance computers–short range forces," Computer Physics Communications, vol. 182, no. 4, pp. 898-911.
6- Brown, W. M., Kohlmeyer, A., Plimpton, S. J., and Tharrington, A. N.,2012, "Implementing molecular dynamics on hybrid high performance computers–Particle–particle particle-mesh," Computer Physics Communications, vol. 183, no. 3, pp. 449-459.
7- Bussi, G., Donadio, D., and Parrinello, M.,2007, "Canonical sampling through velocity rescaling," The Journal of chemical physics, vol. 126, no. 1, p. 014101.
8- Ermilova, I. and Lyubartsev, A. P.,2016, "Extension of the Slipids Force Field to Polyunsaturated Lipids," The Journal of Physical Chemistry B, vol. 120, no. 50, pp. 12826-12842.
9- Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G.,1995, "A smooth particle mesh Ewald method," The Journal of chemical physics, vol. 103, no. 19, pp. 8577-8593.
10- Fricker, G. et al.,2010, "Phospholipids and lipid-based formulations in oral drug delivery," Pharmaceutical research, vol. 27, no. 8, pp. 1469-1486.
11- Hess, B.,2008, "P-LINCS: A parallel linear constraint solver for molecular simulation," Journal of Chemical Theory Computation, vol. 4, no. 1, pp. 116-122.
12- Hills Jr, R. D. and McGlinchey, N.,2016, "Model parameters for simulation of physiological lipids," Journal of computational chemistry, vol. 37, no. 12, pp. 1112-1118.
13- Humphrey, W., Dalke, A., and Schulten, K.,1996, "VMD: visual molecular dynamics," Journal of molecular graphics, vol. 14, no. 1, pp. 33-38.
14- Ingólfsson, H. I. et al.,2014, "The power of coarse graining in biomolecular simulations," Wiley Interdisciplinary Reviews: Computational Molecular Science, vol. 4, no. 3, pp. 225-248.
15- Izvekov, S. and Voth, G. A.,Feb 24 2005, "A multiscale coarse-graining method for biomolecular systems," J Phys Chem B, vol. 109, no. 7, pp. 2469-73.
16- Izvekov, S. and Voth, G. A.,2009, "Solvent-free lipid bilayer model using multiscale coarse-graining," The Journal of Physical Chemistry B, vol. 113, no. 13, pp. 4443-4455.
17- Jämbeck, J. P. and Lyubartsev, A. P.,2012, "An extension and further validation of an all-atomistic force field for biological membranes," Journal of chemical theory computation, vol. 8, no. 8, pp. 2938-2948.
18- Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. J. T. J. o. c. p.,1983, "Comparison of simple potential functions for simulating liquid water," vol. 79, no. 2, pp. 926-935.
19- Li, J. et al.,2015, "A review on phospholipids and their main applications in drug delivery systems," Asian journal of pharmaceutical sciences, vol. 10, no. 2, pp. 81-98.
20- Lu, L. and Voth, G. A.,2009, "Systematic coarse-graining of a multicomponent lipid bilayer," The Journal of Physical Chemistry B, vol. 113, no. 5, pp. 1501-1510.
21- Lyubartsev, A. P. and Laaksonen, A.,Oct 1995, "Calculation of effective interaction potentials from radial distribution functions: A reverse Monte Carlo approach," Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics, vol. 52, no. 4, pp. 3730-3737.
22- Marrink, S.-J. and Mark, A. E.,2004, "Molecular view of hexagonal phase formation in phospholipid membranes," Biophysical journal, vol. 87, no. 6, pp. 3894-3900.
23- Marrink, S. J., De Vries, A. H., and Mark, A. E.,2004, "Coarse grained model for semiquantitative lipid simulations," The Journal of Physical Chemistry B, vol. 108, no. 2, pp. 750-760.
24- Martínez, L., Andrade, R., Birgin, E. G., and Martínez, J. M.,2009, "PACKMOL: a package for building initial configurations for molecular dynamics simulations," Journal of computational chemistry, vol. 30, no. 13, pp. 2157-2164.
25- Mirzoev, A. and Lyubartsev, A. P.,2013, "MagiC: software package for multiscale modeling," Journal of chemical theory computation, vol. 9, no. 3, pp. 1512-1520.
26- Orsi, M. and Essex, J. W.,2013, "Physical properties of mixed bilayers containing lamellar and nonlamellar lipids: insights from coarse-grain molecular dynamics simulations," Faraday discussions, vol. 161, pp. 249-272.
27- Parrinello, M. and Rahman, A.,1980, "Crystal structure and pair potentials: A molecular-dynamics study," Physical Review Letters, vol. 45, no. 14, p. 1196.
28- Plimpton, S.,1995, "Fast parallel algorithms for short-range molecular dynamics," Journal of computational physics, vol. 117, no. 1, pp. 1-19.
29- Reith, D., Pütz, M., and Müller‐Plathe, F.,2003, "Deriving effective mesoscale potentials from atomistic simulations," Journal of computational chemistry, vol. 24, no. 13, pp. 1624-1636.
30- Riniker, S., Allison, J. R., and van Gunsteren, W. F.,2012, "On developing coarse-grained models for biomolecular simulation: a review," Physical Chemistry Chemical Physics, vol. 14, no. 36, pp. 12423-12430.
31- Rosenberger, D., Hanke, M., and van der Vegt, N. F.,2016, "Comparison of iterative inverse coarse-graining methods," The European Physical Journal Special Topics, vol. 225, no. 8-9, pp. 1323-1345.
32- Rühle, V., Junghans, C., Lukyanov, A., Kremer, K., and Andrienko, D.,2009, "Versatile object-oriented toolkit for coarse-graining applications," Journal of Chemical Theory Computation, vol. 5, no. 12, pp. 3211-3223.
33- Shell, M. S.,2008, "The relative entropy is fundamental to multiscale and inverse thermodynamic problems," The Journal of chemical physics, vol. 129, no. 14, p. 144108.
34-Shelley, J. C., Shelley, M. Y., Reeder, R. C., Bandyopadhyay, S., and Klein, M. L.,2001, "A coarse grain model for phospholipid simulations," The Journal of Physical Chemistry B, vol. 105, no. 19, pp. 4464-4470.
35- Siegel, D.,1984, "Inverted micellar structures in bilayer membranes. Formation rates and half-lives," Biophysical journal, vol. 45, no. 2, pp. 399-420.
36- Sodt, A. J. and Head-Gordon, T.,May 28 2010, "An implicit solvent coarse-grained lipid model with correct stress profile," J Chem Phys, vol. 132, no. 20, p. 205103.
37- Van Der Spoel, D. and van Maaren, P. J.,2006, "The origin of layer structure artifacts in simulations of liquid water," Journal of chemical theory computation, vol. 2, no. 1, pp. 1-11.
38- Wang, Z. J. and Deserno, M.,Sep 2 2010, "A systematically coarse-grained solvent-free model for quantitative phospholipid bilayer simulations," Phys Chem B, vol. 114, no. 34, pp. 11207-20.