Розробка технології передбачення спорідненості органічних сполук до молекулярно-біологічних мішеней для рецепторно-орієнтованого скринінгу in silico

Размещено на http://allbest.ru

НАЦІОНАЛЬНА АКАДЕМІЯ НАУК УКРАЇНИ

ІНСТИТУТ МОЛЕКУЛЯРНОЇ БІОЛОГІЇ ТА ГЕНЕТИКИ

УДК 579.69 : 577.12

Автореферат

дисертації на здобуття наукового ступеня

кандидата біологічних наук

Розробка технології передбачення спорідненості органічних сполук до молекулярно-біологічних мішеней для рецепторно-орієнтованого скринінгу in silico

03.00.20 - біотехнологія

ЯКОВЕНКО Олександр Ярославович

Київ - 2008

Дисертацією є рукопис. Роботу виконано у відділі комбінаторної хімії Інституту молекулярної біології та генетики НАН України (м. Київ).

Науковий керівник:кандидат хімічних наук, старший науковий співробітник Бджола Володимир Григорович, Інститут молекулярної біології та генетики НАН України, старший науковий співробітник відділу комбінаторної хімії.

Науковий консультант:доктор біологічних наук, професор, член-кореспондент НАН України Говорун Дмитро Миколайович, Інститут молекулярної біології та генетики НАН України, заступник директора інституту з наукової роботи, завідувач відділу молекулярної та квантової біофізики.

Офіційні опоненти: доктор біологічних наук, професор Тукало Михайло Арсентійович, Інститут молекулярної біології та генетики НАН України, заступник директора інституту з наукової роботи, завідувач відділу ензимології білкового синтезу;

кандидат біологічних наук Нипорко Олексій Юрійович, Інститут клітинної біології та генної інженерії НАН України, науковий співробітник відділу клітинної біології та біотехнології.

Захист дисертації відбудеться 28 жовтня 2008 р. о 10⁰⁰ годині на засіданні спеціалізованої вченої ради Д 26.237.01 Інституту молекулярної біології та генетики НАН України за адресою: 03143, м. Київ, вул. Академіка Заболотного, 150. З дисертацією можна ознайомитись у бібліотеці Інституту молекулярної біології та генетики НАН України за адресою: 03143, м. Київ, вул. Академіка Заболотного, 150.

Автореферат розіслано _23_ вересня 2008 року.

Вчений секретар спеціалізованої вченої ради,

кандидат біологічних наук О. В. Підпала

ЗАГАЛЬНА ХАРАКТЕРИСТИКА РОБОТИ

Актуальність теми. Теоретичні методи обчислювальної біології моделюють перебіг фізико-хімічних процесів у природі і дозволяють передбачувати результати реальних біологічних експериментів. Так, зокрема, сучасні комп'ютерні розрахунки спорідненості низькомолекулярних органічних молекул до молекулярно-біологічних мішеней спроможні частково замінити біологічний скринінг хімічних сполук віртуальним скринінгом (ВС) in silico. Застосування ВС для, наприклад, ідентифікації сполук із бажаними біологічними властивостями суттєво пришвидшує та здешевлює розробку нових ліків (Muthas at al., 2008).

Останнім десятиліттям завдяки прогресу обчислювальної техніки та стрімкому накопиченню експериментальних даних про структуру біомолекул найефективнішим та найпоширенішим комп'ютерним методом ідентифікації молекул із бажаними біологічними властивостями є високопродуктивний (ВП) рецепторно-орієнтований (РО) ВС.

Це один із методів обчислювальної біології, який забезпечує потокове передбачення in silico спорідненості низькомолекулярних сполук до молекулярно-біологічних мішеней моделюванням структури їхніх просторових комплексів та подальшої оцінки енергії міжмолекулярних взаємодій, які відповідають за їхню стабілізацію.

Сучасні квантово-хімічні методи ab initio здатні розраховувати молекулярну спорідненість з будь-якою бажаною точністю, але обчислювальна складність не дозволяє застосовувати їх для ВС хімічних баз великої розмірності.

Популярні нині методи так званих оціночних функцій ВС наближено моделюють міжмолекулярні взаємодії, проте їхня обчислювальна складність є задовільною. Оскільки точність оціночних функцій ВС залежить більшою мірою від базових наближень, на яких вони ґрунтуються, ніж від їхньої обчислювальної складності, то це, в принципі, дозволяє розробити достатньо точні методи при збереженні швидкодії, прийнятної для ВС великих масивів.

Таким чином, розробка нових ефективних моделей молекулярної спорідненості органічних сполук до молекулярно-біологічних мішеней для потреб потокового РОВС in silico є актуальною біотехнологічним завданням.

Зв'язок роботи з науковими програмами, планами, темами. Дисертаційна робота виконувалася за конкурсними тематиками «Розробка інгібіторів протеїнкіназ CDK2 та СК2 як потенційних ліків проти онкологічних захворювань» (номер держреєстрації №0106U000543, 2004-2006 рр.) та «Оптимізація інгібіторів протеїнкінази СК2 та дослідження їхньої біологічної активності на культурах ракових клітин» (номер держреєстрації 0107U004939, 2007-2009 рр.), а також за фундаментальною тематикою «Вивчення протеїнкіназ як молекулярних мішеней для розробки терапевтичних засобів методами комбінаторної хімії та комп'ютерного моделювання» (номер держреєстрації №0107U003345, 2008-2012 рр.).

Мета і завдання дослідження. Мета дисертаційної роботи - створення на базі нових фізико-хімічних моделей оригінальних методів потокового РО передбачення афінності органічних сполук до молекулярно-біологічних мішеней з поліпшеними прогностичними можливостями та розробка обчислювально простих і ефективних алгоритмів їхньої комп'ютерної реалізації.

Для досягнення цієї мети було поставлено та вирішено такі завдання:

1. Розробити модель потокової комп'ютерної оцінки спорідненості малих органічних молекул до молекулярно-біологічних мішеней у водному середовищі.

2. Створити ефективну модель для комп'ютерної оцінки енергії міжмолекулярної взаємодії низькомолекулярних органічних сполук з молекулярно-біологічними мішенями.

3. Запропонувати метод комп'ютерної оцінки енергії, необхідної для перенесення низькомолекулярної сполуки із водного середовища у неполярну, гідрофобну кишеню центру її зв'язування з мішенню.

4. Розробити обчислювально нескладний і водночас ефективний алгоритм передбачення афінності органічної молекули до молекулярно-біологічної мішені, який адекватно враховує ентропійні ефекти комплексоутворення.

5. Перевірити прогностичні можливості розроблених алгоритмів на репрезентативній виборці сполук, для яких відомі експериментальні константи зв'язування з відповідними молекулярно-біологічними мішенями.

6. Спроектувати та випробувати лабораторний програмний комплекс для РО ВПВС in silico.

Об'єкт дослідження: технологія передбачення афінності низькомолекулярних органічних сполук до молекулярно-біологічних мішеней, включаючи алгоритми їхньої комп'ютерної реалізації.

Предмет дослідження: точність і обчислювальна складність алгоритмів передбачення афінності низькомолекулярних органічних сполук до молекулярно-біологічних мішеней, адекватність моделей та методів моделювання міжмолекулярних взаємодій.

Методи дослідження: неемпіричні методи квантової хімії; числові методи пошуку екстремумів функцій багатьох змінних: симплекс-метод, Ньютона та l-BFGS; методи оцінки суперпозицій множин: включень-виключень і функцій розподілу молекулярних станів; геометричні алгоритми оцінки суперпозиції точок та багатомірних поверхонь у багатомірному просторі: імітація невиродженості, тріангуляція Делоне, -сфери, геометричні примітиви обчислення площ і об'ємів перетину сфер; метод Кіргофа обчислення струмів у електричних колах; метод моделювання потенціальної енергії емпіричними силовими полями; метод оцінки енергії десольватації із застосуванням неявного розчинника GBSA; мови програмування.

Наукова новизна одержаних результатів:

1. Вдосконалено метод релаксації електронегативностей за Кіргофом: сформульовано та параметризовано нову модель емпіричного розрахунку точкових атомних зарядів органічних сполук шляхом релаксації електронегативностей за формалізмом Кіргофа.

2. Створено нове силове поле для оцінки енергії міжмолекулярних взаємодій низькомолекулярних органічних сполук із молекулярно-біологічним мішенями на основі розробленого методу розрахунку точкових атомних зарядів.

3. Розроблено нову модифікацію GBSA-методу обчислення енергії, необхідної для перенесення низькомолекулярної органічної сполуки з водного середовища у неполярну кишеню центру зв'язування біополімеру.

4. Запропоновано оригінальну модель передбачення афінності низькомолекулярних органічних сполук до молекулярно-біологічних мішеней, яка адекватно враховує зміну ентропії при утворенні комплексу.

5. Для потреб ВПВС in silico створено обчислювально прості алгоритми комп'ютерного передбачення афінності низькомолекулярних органічних сполук до молекулярно-біологічних мішеней та перевірено їх ефективність на репрезентативній вибірці лігандів, для яких відомі експериментально встановлені константи зв'язування.

6. Створено лабораторний програмний комплекс для потокової оцінки спорідненості низькомолекулярних органічних сполук до молекулярно-біологічних мішеней. Із застосуванням розробленої технології ідентифіковано нові високоефективні інгібітори протеїнкінази СК2.

Практичне значення одержаних результатів. У дисертації розроблено нові моделі та обчислювально нескладні алгоритми з поліпшеною точністю передбачення афінності низькомолекулярних органічних сполук до центрів їхнього зв'язування з біополімерами. Ці розробки можуть бути використані:

- для РО ВПВС у фармації;

- як основа для подальшого вдосконалення модельних засад оцінки енергії міжмолекулярних взаємодій та передбачення молекулярної афінності;

- для створення новітніх програмно-обчислювальних комплексів РОВС in silico із залученням додаткових алгоритмів молекулярного докінгу.

Особистий внесок здобувача. Автором розроблено модель та проведено кодування алгоритмів розрахунку функцій розподілу молекулярних станів, переформульовано теорію, проведено параметризацію нової моделі та кодування алгоритмів розрахунку точкових атомних зарядів. Дисертантом власноруч модифіковано GBSA-метод оцінки десольватаційних енергій органічних сполук та створено алгоритм його комп'ютерної реалізації.

Дисертантом особисто проведено всі обчислювальні експерименти, пов'язані з молекулярним моделюванням in silico, а також створено лабораторний програмний комплекс для потокової оцінки спорідненості низькомолекулярних органічних сполук до молекулярно-біологічних мішеней з адекватним врахуванням ентропійних ефектів для потреб ВПВС із залученням додаткового алгоритму молекулярного докінгу DOCK.

Постановку наукових завдань та інтерпретацію отриманих результатів здійснено спільно з науковим керівником к.х.н., с.н.с. В. Г. Бджолою, науковим консультантом д.б.н., проф., чл.-кор. НАН України Д. М. Говоруном та зав. відділу комбінаторної хімії д.х.н., проф. С. М. Ярмолюком. Теорію обчислень часткових атомних зарядів створено спільно з групою обчислювальної хімії кафедри органічної хімії хімічного факультету Московського державного університету імені М.В.Ломоносова (к.х.н. О. О. Оліференко та ін.).

Апробація результатів дисертації. Матеріали дисертаційної роботи доповідалися на 4-ій Міжнародній конференції «Інгібітори протеїнкіназ 2005» (Варшава, Польща, 2005), Конференції молодих вчених в Інституті молекулярної біології та генетики НАН України (Київ, Україна, 2007) та І-ій Всеросійській конференції «Многомасштабное моделирование процессов и структур в нанотехнологиях» (Москва, Росія, 2008).

Результати дисертаційної роботи опубліковано у 4 статтях у наукових фахових журналах, а також представлені на 3 наукових конференціях у вигляді усних доповідей, стендових презентацій та тез.

Структура та обсяг дисертації. Дисертація складається із вступу, огляду літератури, матеріалів та методів дослідження, експериментальної частини, аналізу та узагальнення результатів роботи, висновків, переліку використаних літературних джерел, який нараховує 115 найменувань. Дисертація містить 5 рисунків та 13 таблиць. Загальний обсяг дисертації складає 154 сторінок машинописного тексту.

ОСНОВНИЙ ЗМІСТ РОБОТИ

атомний органічний протеїнкіназа

Матеріали і методи. Для обчислення констант інгібування ферменту СК2 органічними сполуками використано їхні кристалографічні тривимірні моделі 1zog, 1om1, 1f0q, 1m2r, 1m2p, 1m2q і 1j91 (www.pdb.org).

Для оптимізації геометрії і розрахунку електричних дипольних та квадрупольних моментів застосовували неемпіричні методи квантової хімії на рівні теорії HF/6-311++G(d,p) за допомогою вільно поширюваного програмного пакету "GAMESS".

Програма DOCK 4.0 використана для побудови моделей молекулярних комплексів органічних сполук з молекулярно-біологічною мішенню.

Власноруч створені програмні продукти використано для:

- мінімізації значень функцій багатьох змінних методом деформованого багатогранника (симплексу), методом найкоротшого спуску, квазіньютонівським методом l-BFGS та класичним методом Ньютона;

- розрахунку атомних зарядів та електричних моментів органічних сполук;

- оптимізації геометрії молекулярних систем методом атомних силових полів;

- розрахунку площ поверхонь і об'ємів молекул, доступних і недоступних для розчинника (води) відповідно;

- обчислень енергії десольватації низькомолекулярних органічних сполук;

- визначення функцій розподілу станів лігандів за обертально-поступальними ступенями вільності;

- розрахунку розподілу станів лігандів за внутрішніми ступенями вільності;

- визначення сталих зв'язування органічних сполук з молекулярно-біологічними мішенями на основі відомих просторових моделей їхніх молекулярних комплексів.

Програми “Origin”, “Excel” і “Mathematica” застосовано на стадії розробки математичного формалізму та інтерпретації отриманих in silico результатів.

РЕЗУЛЬТАТИ ДОСЛІДЖЕННЯ ТА ЇХНЄ ОБГОВОРЕННЯ

Моделювання міжмолекулярних електростатичних взаємодій. Типовою моделлю електростатичних взаємодій для ВПВС є кулонівська взаємодія точкових зарядів, що моделюють атомні ядра. Через значні розміри біомолекулярних систем високоточні методи квантової хімії для моделювання міжмолекулярних взаємодій у ВПВС практично не застосовують.

Натомість використовують метод силових полів (СП) із кулонівською моделлю електростатичних взаємодій. Такий підхід є обчислювально простим і потребує лише відомостей про розподіл електричного заряду на атомах сполук, що аналізуються.

Однак якісні моделі електричного розподілу відомі лише для незначної кількості сполук, таких як амінокислоти, нуклеїнові кислоти та деякі інші біохімічно важливі молекули, тоді як висопродуктивний ВС покликаний аналізувати органічні сполуки, відомості про які, зазвичай, зводяться лише до їхніх структурних формул. Для знаходження розподілу електричного заряду на атомах довільної органічної молекули, виходячи з її структурної формули, автором розроблено та параметризовано модель релаксації електронегативностей за Кіргофом.

Цей метод є емпіричним. На противагу неемпіричним квантово-хімічним методам, які оперують складними обчисленнями функцій електронної густини, емпіричні методи базуються на емпірично підібраних числових параметрах. Вони обчислювально прості, але потребують належної параметризації. Метод релаксації електронегативностей представляє розподіл електричної густини як дискретний набір точкових зарядів, локалізованих на атомах молекули.

Метод релаксації електронегативностей вперше запропоновано колективом кафедри органічної хімії хімічного факультету Московського державного університету імені М.В.Ломоносова та детально описано у літературі (Oliferenko at al., 2006). Однак його теоретична частина містить низку недоліків. Так, зокрема, для коректного відтворення електронейтральності застосовують додаткові функції перенормування розрахованих зарядів, що є причиною систематичної похибки методу та зниження його швидкодії.

Нами вдосконалено формалізм методу релаксації електронегативностей, звівши його до розв'язку рівнянь електростатики і переформулювавши процедуру релаксації електронегативностей у термінах розв'язку класичних рівнянь Лапласа. Хоча формалізм Кіргофа залишився ключовим математичним апаратом моделі, запропоновані модифікації дали змогу усунути систематичну похибку перенормування та спростити обчислювальну складність, усунувши поправочні функції.

Формалізм розрахунку точкових зарядів методом релаксації електронегативностей за Кіргофом. Якщо органічна сполука не має заряду, то застосовують такий формалізм:

² ( r ) = 0. ( 1 )

Якщо ж сполука має конститутивний заряд (прикладом може бути кислота або основа в іонізованому стані), то формалізм набуває складнішого вигляду:

² ( r ) = -, ( 2 )

де ² - оператор Лапласа, - електронегативність у формі електричного потенціалу, r - радіус-вектор атома, - густина конститутивного заряду.

Ефективна жорсткість хімічних зв'язків у запропонованій нами квадратичній моделі релаксації електронегативностей є аналогом провідності у теорії електричних кіл. Ефективні хімічні жорсткості враховуються за допомогою матриці вагових коефіцієнтів при конструюванні дискретної форми оператора Лапласа:

L = B^T W B, ( 3 )

де L - дискретна форма оператора Лапласа, задана матрицею розміру nn, B - матриця інцидентності молекулярного графа розміру nm, W - матриця вагових коефіцієнтів каналів розміром mm, n - кількість атомів, а m - кількість ковалентних зв'язків аналізованої сполуки. Для розв'язку рівнянь (1) та (2) застосовується рівняння (3), яке перетворює точкові потенціали електричного поля (електронегативності) у дискретні струми (точкові заряди на атомах). Таким чином, для вектору точкових зарядів отримано розрахункову формулу:

q = ( L^-1 - I ) , ( 4 )

де q - вектор точкових атомних зарядів, І - одинична діагональна матриця розміру nn, а - вектор атомних електронегативностей.

Параметризація моделі. Для параметризації розробленої моделі застосовано неемпіричні квантово-хімічні методи на рівні теорії HF/6-311++G(d,p). Оскільки квантово-хімічні методи оперують функціями електронної густини, а не точковими атомними зарядами, то прямо оцінити похибку відтворення розподілу електронної густини розробленим нами методом неможливо. З іншого боку, цікавішою для практичного використання є точність відтворення електростатичних міжмолекулярних взаємодій. Для оцінки точності відтворення електростатичних міжмолекулярних взаємодій нами застосовано дипольні та квадрупольні електричні моменти. Електричні моменти визначають енергію міжмолекулярної електричної взаємодії і їх можна розрахувати як з квантово-хімічних функцій розподілу електронної густини, так і на основі точкових зарядів. Для розрахунку електричних моментів на основі квантово-хімічних функцій розподілу електронної густини використано вбудовану команду програми "GAMESS" (IEMOM=2), для розрахунку електричних моментів на основі точкових атомних зарядів використано такі рівняння:

( 5 )

де d і Q - декартові компоненти дипольного та квадрупольного електричних моментів відповідно; q - точковий заряд атома; x, y та z - декартові координати атомів у системі центру мас сполуки; індекс i пробігає всі атоми сполуки.

Для параметризації моделі розподілу атомних зарядів використано оціночну функцію:

S = N^-1 ( d^a_i - d^e_i )² + Ѕ ( Q^a_ii - Q^e_ii )²,

де індексом е позначено параметри, розраховані неемпіричними методами квантової хімії, а індексом а - компоненти моментів, розраховані оригінальним емпіричним методом, індекс і - індекс сумування, а N - кількість сполук у навчальній вибірці.

Якість параметризації моделі. Середня похибка відтворення запропонованим нами методом дипольних моментів у навчальній вибірці для моделі дипольно-квадрупольної параметризації склала 0,28 D, що майже збігається з типовою експериментальною похибкою.

Величина лінійної кореляції емпіричних і неемпіричних дипольних моментів доволі висока (r = 0,948), в той час як до симплекс-оптимізації цей параметр складає r = 0,630. Є декілька структур - індол, пірол, нітробензол і піридазин, дипольні моменти яких відтворюються з похибкою, що перевищує 1 D. Вірогідною причиною цього є значний перерозподіл -електронної густини, який важко моделюється без залучення додаткових алгоритмів. Незважаючи на декілька нетипових випадків, електричні моменти всіх проаналізованих нами хімічних класів молекул розраховувалися із задовільною точністю.

Так, дипольні моменти алканів, спиртів, альдегідів, кетонів та сірковмісних сполук відтворено з середньою похибкою <0,2 D. Дипольні моменти галогеноалканів, естерів, амінів, карбонових кислот та амідів розраховано з середньою похибкою <0,3 D, а аналогічні показники ароматичних вуглеводнів, ароматичних гетероциклів та фосфоровмісних сполук оцінено з середньою похибкою <0,4 D.

Запропонований нами метод відтворює з прийнятною точністю і квадрупольні моменти. Загалом середня похибка по навчальній вибірці відтворення компонентів квадрупольного тензора Q_xx, Q_yy і Q_zz склала 0,171, 0,165 і 0,186 В відповідно. Через менший внесок квадрупольного моменту у величину електричного потенціалу похибка його експериментального визначення перевищує теоретичні похибки. При розгляді середніх похибок передбачення за класами сполук виявилося, що їхні величини для диполів та квадруполів не збігаються.

Такі хімічні класи сполук, як галогеноалкани, аміни та карбоксикислоти характеризуються відмінно відтвореним дипольним моментом, але точність відтворення квадрупольного моменту для них невисока. З іншого боку, для амідів, нітро- та ціаносполук, а також для ароматичних гетероциклів квадрупольні моменти розраховуються з більшою точністю, ніж дипольні. Таким чином, можна зробити висновок, що загальна тенденція параметризації виявляється у тому, що більші помилки передбачення дипольних моментів хімічних сполук супроводжуються підвищеною точністю відтворення їхніх квадрупольних моментів і навпаки.

Порівняння з іншими методами розрахунку точкових зарядів. Ми порівняли заряди на атомах, розраховані вдосконаленим нами методом релаксації електронегативностей за Кіргофом для широкої виборки органічних молекул, з аналогічними літературними даними, отриманими низкою неемпіричних методів, а саме - аналізом заселеностей за Мулікеном (MPA), природнім аналізом заселеностей (NPA), AIM зарядами Батдера, методом узагальнених атомних полярних тензорів (GAPT), методом CHELPG у параметризації проти ESP-зарядів. Крім того, для порівняння брали також точкові заряди, розраховані програмою PETRA за емпіричним методом Гастейгера PEOE, який вважається нині найточнішим емпіричним методом.

Таблиця 1

Порівняння методів розрахунку розподілу електронної густини

Метод	r	a	b
PEOE	0,978	-0,0034	0,6313
MPA	0,962	0,0025	1,2127
NPA	0,943	-0,0214	1,5582
AIM	0,880	0,1141	2,4469
GAPT	0,857	0,0999	1,6613
CHELPG	0,927	-0,0113	1,1594

Порівняно з даними, отриманими методом РЕОЕ, нашому підходу притаманна більша поляризація атомних зарядів, як у неемпіричних методів. Це ілюструє статистичний аналіз лінійної кореляції отриманих нами зарядів з аналогічними даними, отриманими іншими методами. Кореляційні регресії наведено у таблиці 1 у відповідності з рівнянням y = a + bx, де y - заряди, отримані альтернативними методами, x - авторські заряди; а a та b -коефіцієнти, розраховані методом найменших квадратів.

Коефіцієнти b більші за одиницю для зарядів, розрахованих методами MPA, NPA, AIM, GAPT та CHELPG: це узгоджується з відомим фактом завищення абсолютних значень точкових зарядів відомими неемпіричними методами. За значеннями коефіцієнтів b авторські заряди відрізняються лише на 10% та 20% від даних неемпіричних методів CHELPG та MPA відповідно. В той же час малі значення коефіцієнта a вказують на низьку похибку у визначенні авторських зарядів, що відповідає високій кореляції 0,927 та 0,962 з даними методів CHELPG та MPA відповідно.

Це доводить, що принаймні для наведеної вибірки із 13 молекул запропонований алгоритм розрахунку атомних зарядів задовільно відтворює точкові заряди, отримані на основі аналізу функцій електронної густини, обчислених неемпіричними методами CHELPG та MPA.

На противагу задовільній узгодженості з методами CHELPG та MPA, атомні заряди, розраховані методами AIM та GAPT, відтворюються авторським методом гірше (коефіцієнти кореляції 0,880 та 0,857 відповідно). Однак AIM та GAPT - це методи, результати яких вкрай погано відтворюються емпіричними схемами. Кореляція з NPA-розрахунками займає проміжне значення між попередніми двома групами методів, що ґрунтуються на аналізі функцій електронної густини, розрахованих неемпіричними методами. Порівняння авторських зарядів з результатами методу РЕОЕ показало, що вони доволі сильно корелюють між собою, але РЕОЕ розраховує менш полярні точкові заряди, ніж авторський: про це свідчить значення b=0,6313.

Обчислювальна складність розробленої моделі визначається найскладнішим її етапом - оберненням матриць розміру NхN і оцінюється щонайбільше, як N³, де N - кількість атомів сполуки. Усі відомі автору методи розрахунку точкових атомних зарядів мають значно вищу обчислювальну складність, ніж запропонований.

Таким чином, нами переформульовано теоретичну частину методу розрахунку атомних зарядів шляхом релаксації електронегативностей за Кіргофом та розроблено нову, обчислювально нескладну модель розрахунку атомних зарядів органічних молекул.

Модель належним чином параметризовано методом деформованого багатогранника та перевірено на похибку відтворення електричних моментів порівняно з неемпіричними методами.

Розроблений метод є універсальним і вирізняється з-поміж відомих емпіричних методів найбільшою швидкодією та найвищою точністю відтворення розподілу зарядів на атомах органічних молекул. Він є найприйнятнішим для застосування у ВПВС, оскільки потребує лише відомостей про структурну формулу органічної сполуки, а його швидкодія у задачах потокового ВПВС in silico сумісна з сучасним рівнем розвитку обчислювальної техніки.

Силове поле YFF. Для моделювання взаємодій у молекулярних системах за участі молекулярно-біологічних мішеней найчастіше застосовують метод молекулярно-механічних СП. Загальноприйнятою в методі СП є використання моделі точкових центрів взаємодії, що моделюють атоми. Цим методом молекулярні взаємодії моделюються трьома складовими: ковалентними (в межах однієї молекули), Ван-дер-Ваальсовими та електростатичними взаємодіями між атомами сусідніх молекул. Потенціали ковалентних та Ван-дер-Ваальсових взаємодій у методі СП визначаються структурою сполук та параметризацією СП. Моделлю електростатичних взаємодій у методі СП є кулонівська модель точкових зарядів.

Попри велику кількість молекуляро-механічних СП, жодне з них досі не параметризовано на точне відтворення електричних моментів. Оскільки електричні моменти характеризують потенціал електричного поля у просторі, то їхнє адекватне відтворення дозволяє оцінювати електростатичні взаємодії з високою точністю.

Таким чином, застосування методу розрахунку точкових атомних зарядів за схемою релаксації електронегативностей за Кіргофом у методі СП дозволяє поліпшити точність розрахунку як міжмолекулярних взаємодій у цілому, так і електростатичних взаємодій зокрема.

Базову параметризацію ковалентних та Ван-дер-Ваальсових взаємодій ми запозичили з СП, спеціально параметризованого на відтворення взаємодій у системах "біополімер - низькомолекулярний органічний ліганд" MMFF94 (Halgren, 1995). Комбінацією базової параметризації MMFF94 та авторської схеми розрахунку точкових зарядів створено нове СП YFF1. Застосування авторської схеми розрахунку точкових зарядів дозволяє приблизно на порядок покращити точність моделювання електростатичних взаємодій, заявлених у СП MMFF94 (середня похибка відтворення електричних моментів 0,8 - 1,2 D).

При створенні СП YFF1 використано атомну типізацію методу розрахунку атомних зарядів. Це СП спроектовано для автоматичної компіляції його параметрів, виходячи лише з молекулярного графу органічної сполуки.

З метою створення універсального СП при перенесенні параметрів MMFF94 також вилучено ті типи атомів, параметризація яких у MMFF94 спиралася на топологічну структуру особливих підграфів вихідного молекулярного графу.

Таким чином, комбінацією базової параметризації MMFF94 з авторською схемою розрахунку атомних зарядів нами створено нове універсальне СП YFF1, яке найкраще з-поміж інших СП оцінює енергію міжмолекулярних взаємодій у системах "молекулярно-біологічна мішень - низькомолекулярна органічна сполука". Створено програмний комплекс, який забезпечує автоматичну компіляцію довільної органічної хімічної структури у СП YFF1.

Визначення in silico енергії десольватації органічних сполук. Оцінка спорідненості ліганду до рецептора не може бути проведена лише на основі оцінки вільної енергії його взаємодії з рецептором. Для цього необхідно обов'язково враховувати різницю взаємодій дисоційованих та асоційованих молекулярних ансамблів з оточуючим середовищем - водою.

У задачах ВПВС, однак, можна спростити розгляд проблеми різниці взаємодії ансамблів з водою до питання про різницю між взаємодіями ліганду із водою та з неполярним середовищем кишені центру зв'язування біополімеру.

Взаємодія молекули з водним середовищем - це та енергія, яка виділиться (поглинеться) при перенесенні молекули з водного середовища у вакуум - так звана енергія десольватації. У експериментах in silico ВПВС розміри біологічного полімеру на кілька порядків перевищують розмір лігандів. Це дозволяє знехтувати різницею енергій десольватації окремих комплексів одного й того ж самого біополімеру з різними синтетичними низькомолекулярними лігандами.

Крім того, очевидно, що енергія десольватації таких комплексів буде приблизно рівною десольватаційній енергії вільного біополімеру. Таким чином, єдина складова, для якої енергія десольватації суттєво змінюється при комплексоутворенні, є десольватаційна енергія органічного ліганду. Саме різниця значень десольватаційних енергій низькомолекулярного ліганду визначає спорідненість, зумовлену взаємодією з водним середовищем.

Поширеною моделлю для швидкої оцінки десольватаційної енергії є розрахунок енергії, необхідної для перенесення електричного поля, створеного лігандом, з гомогенного середовища з діелектричною проникністю розчинника ( для води ~ 80 ) у середовище з діелектричною проникністю вакууму ( чи середовище з діелектричною проникністю близькою до одиниці ).

Очевидно, що така модель є чутливою до точності опису розподілу електричного поля навколо ліганду. Для корекції високої десольватаційної енергії, притаманної незарядженим алканам, модель доповнюють корекцією на площу поверхні атомів, доступну для молекул розчинника. Така модель розрахунку десольватаційних енергій отримала назву Generalized Born Surface Area (GBSA). Ми її обрали за основу для розрахунків десольватаційної енергії ліганду у водному розчині.

Нами модифіковано формалізм обчислення десольватаційних енергій порівняно з моделлю GBSA. Це пов'язано із спробами поліпшення швидкодії алгоритму через введення обчислення ефективного радіуса Борна з параметрів об'ємів атомів, які обчислювально тісно пов'язані з обчисленням площі поверхні сполуки, доступної для розчинника.

Використання точних площ та об'ємів, розрахованих за допомогою сучасних алгоритмів обчислювальної геометрії (тріангуляції Делоне, імітації невиродженості, примітивів обчислення площ і поверхонь перетинів множин сфер), замість їхніх наближених значень дозволяє спростити формалізм обчислення радіуса Борна і уникнути використання коефіцієнтів геометричної регресії. Слід зауважити, що обчислення доступної для розчинника поверхні є складовою частиною розрахунку недоступного для розчинника об'єму і, таким чином, обидві процедури можна об'єднати.

Врешті-решт запропонована формула для розрахунку ефективного атомного радіуса багатоатомної сполуки набуває вигляду:

,

де V_j точний об'єм j-го атома, що моделюється об'єднанням сфер різного діаметра, r_i^vdw Ван-дер-Ваальсів радіус і-го атома, R_solv - радіус сфери, що моделює молекулу розчинника, r_ij - відстань між атомами і та j. Сума береться по всіх n атомах досліджуваної молекули. Знаючи ефективний радіус, легко розрахувати вільну енергію десольватації досліджуваної органічної молекули.

У запропонованому нами методі визначають дві незалежні складові десольватаційної енергії - енергію, необхідну для перенесення сполуки із середовища з діелектричною проникністю розчинника у середовище з діелектричною проникністю вакууму, та енергію, необхідну для впорядкування найближчого шару розчинника. Перша з них визначається через роботу електричного поля та ефективний радіус молекули (радіус Борна), а друга - через лінійну комбінацію площ атомів, доступних для розчинника, та статистичних коефіцієнтів регресії:

де _i - радіус Борна i-го атома, , r_ij - відстань між атомами i та j, _i - коефіцієнт регресії десольватаційної енергії по площі поверхні, доступній для розчинника i-го атома , A_i - площа i-го атома, доступна для розчинника, q_i і q_j - точкові заряди на атомах i та j відповідно, - діелектрична проникність розчинника. Суми взято по всіх n атомах молекули.

Модифіковану нами модель GBSA перевірено на тестовій вибірці, що містила >400 сполук різноманітної природи, для яких відомі експериментальні значення десольватаційної енергії. Для всіх цих молекул розраховано об'єми та площі поверхні, доступної для розчинника, точкові атомні заряди (за модифікованим методом Кіргофа) і теоретично передбачено десольватаційну енергію. Отримане значення коефіцієнта лінійної кореляції між передбаченими та експериментально визначеними значеннями енергії десольватації складає 0,84, що є на сьогодні найвищим серед відомих емпіричних методів сімейства GBSA.

Розрахунок молекулярної спорідненості на основі оцінки енергії міжмолекулярних взаємодій. Апарат статистичної фізики дозволяє визначати імовірності знаходження молекулярних ансамблів у певному стані за умови наявності відомостей про їхні взаємодії. Нами було застосовано формалізм функцій розподілу станів молекулярного ансамблю за ступенями вільності для обчислення молекулярних спорідненостей з урахуванням ентропійних ефектів. Для моделювання спорідненості у ВПВС обрано обчислювально найпростіше наближення - моделювання розподілів на основі розрахунку гармонійних коливань молекулярної системи навколо енергетичного мінімуму. Кількість копій молекулярної системи Z, що знаходиться у певній точці фазового простору х, визначається як: , де U - енергія молекулярної системи, а - обернений температурний фактор =1/kT. Процедура знаходження наближеного розв'язку цього рівняння - розкладання функції U в ряд до членів другого порядку включно поблизу її мінімуму. Існування аналітичного розв'язку дозволяє швидко розраховувати молекулярні спорідненості, коли відомі координати мінімуму функції U та міжмолекулярні взаємодії, тобто значення других похідних по усіх координатах фазового простору.

Автором запропоновано унікальну методику розрахунку функцій розподілу як для обертально-поступальних, так і для внутрішніх ступенів вільності молекулярних комплексів. Для моделювання молекулярних взаємодій використано формалізм СП YFF та метод GBSA.

Критерієм правильності прогнозу міжмолекулярних спорідненостей обрано передбачення експериментально встановлених сталих інгібування (Кі) для інгібіторів ферментів через незалежність цього параметру від умов його визначення, зокрема від концентрації субстрату.

Розрахунок статистичних сум для обертально-трансляційних ступенів вільності ліганду. Розподіл системи за обертально-поступальними ступенями вільності можна вивчати за допомогою функцій розподілу обертально-поступальних станів ліганду, жорстко зв'язавши початок координат з молекулярно-біологічною мішенню, враховуючи той факт, що обертання глобули біополімеру навколо своєї осі є ускладненим через значну товщину шару структурованого розчинника біля його поверхні. Відтак, змінними будуть лише координати ліганду, а для знаходження координат енергетичного мінімуму системи у такому наближенні достатньо оптимізації відносної позиції та орієнтації молекули ліганду в полі біомакромолекули. Оптимізація проводилася за допомогою СП YFF.

Апарат статистичної фізики, який надалі називається довільною моделлю розподілу молекулярних станів за обертально-поступальними ступенями вільності, визначає молекулярну спорідненість К_RT рівнянням:

,

де V - обернена молярна концентрація, Н - детермінант матриці других похідних функції U по обертально-поступальних координатах, індекси f та b позначають вільний та зв'язаний стан ліганду відповідно, U₀ - енергія незв'язаних молекул. Результати обчислень спорідненості за цим рівнянням для вибірки експериментально встановлених моделей тривимірних кристалічних структур комплексів ферменту СК2 з різними інгібіторами наведено в таблиці 2.

З таблиці видно, що квадратичне наближення рядом Тейлора є недостатнім для задовільного моделювання функції потенціальної енергії в активному центрі зв'язування біополімеру.

Отримані значення функції розподілу, попри високу точність відтворення енергії міжмолекулярних взаємодій, приблизно в квадрат разів перевищують відомі значення Кі.

Імовірною причиною похибки є складна структура центру зв'язування біополімеру та несиметричність функціоналів Ван-дер-Ваальсової взаємодії відносно мінімуму.

Таблиця 2

Моделі функцій розподілу за обертально-поступальними ступенями вільності

PDB код	Z, довільна модель розподілу, 10-6	Z, квазі-довільна модель розподілу, 10-6	Ki, експеримент, 10-6
1f0q	8,05	1,36	1,85
1j91	21,02	4,14	0,40
1m2p	10,65	1,89	0,78
1m2q	11,20	1,99	0,80
1m2r	7,81	1,32	0.42
1om1	4,93	0,78	0,17
1zog	9,35	1,62	0,07

Для вирішення цієї проблеми запропоновано таке наближення - знехтувати відштовхувальною гілкою експоненти Ван-дер-Вааль-сового потенціалу. Оскільки гілки модельної параболи в аналітичному розв'язку не відрізняються одна від одної.

Спочатку розраховують другі похідні по ступенях вільності відповідно до обраної моделі міжмолекулярних взаємодій - СП YFF. Потім береться корінь четвертого ступеня з детермінанту.

Запропоновану модель названо моделлю квазідовільного розподілу обертально-поступальних сту-пенів вільності молекулярних комплексів. Формально це записується як рівняння:

,

позначення якого збігаються з попередніми. Як видно з таблиці, застосування такої поправки привело до більш зважених значень Кі порівняно з експериментальними.

Оскільки в такій моделі не враховано розподіл за внутрішніми ступеням вільно-сті, то передбачені значення повинні бути дещо заниже-ними, що й отримано.

Розрахунок функцій розподілу молекулярних станів лігандів за внутрішніми ступенями вільності. Згідно формалізму YFF сполуки мають лише внутрішні ступені вільності, пов'язані з вільним обертанням їхніх фрагментів один відносно одного навколо одинарних валентних зв'язків при кімнатній температурі.

Для моделювання таких розподілів використано той же апарат статистичної фізики, що й для обертально-поступальних.

Нами розроблено комп'ютерний алгоритм, який конструює та розраховує матрицю других похідних по внутрішніх ступенях вільності для всіх пар двогранних кутів молекулярної системи. У зв'язку зі значною кількістю внутрішніх ступенів вільності білкової молекули розгляд її розподілу за внутрішніми ступенями вільності у роботі не проводили.

Для обчислення функцій розподілу за внутрішніми ступенями вільності К_IN застосовували довільний формалізм розподілу, описаний раніше. Для об'єднання функцій розподілу обертально-поступальних та внутрішніх координат використано відоме з літератури наближення:

.

Таблиця 3

Точність теоретичної оцінки Кі

PDB код	Теорія	Експе-римент
1f0q	-5,068	-5,733
1j91	-6,956	-6,398
1m2p	-5,456	-6,108
1m2q	-5,323	-6,097
1m2r	-6,451	-6,377
1om1	-7,220	-6,770
1zog	-8,189	-7,155

Розрахунки К_і інгібіторів протеїнкінази СК2. З метою кількісної оцінки точності передбачення молекулярної спорідненості (таблиця 3) нами розраховано Кі вибірки інгібіторів протеїнкінази СК2, для яких відомі експериментально встановлені тривимірні структури їхніх комплексів з ферментом. Дані теоретичного передбачення спорідненості та експериментальні афінності подані у шкалі десяткових логарифмів. Теоретично передбачені афінності сполук відрізняються між собою суттєвіше, ніж експериментальні значення Кі, однак впорядкування комплексів за відносними афінностями є правильним.

Для ілюстрації кореляції між передбаченими та експериментальними значеннями афінностей побудовано діаграму взаємозалеж-ності між десятковими логарифмами теоретично передба-чених афінностей сполук у вибірці та десятковими логарифмами експе-риментальних Кі. Для цієї діаграми стандартний коефіцієнт кореляції r=0,96.

На основі цих даних можна стверджувати, що метод теоретик-ного передбачення спорідненості низькомолекулярних органічних спо-лук до молекулярно-біологічних мішеней задовільно корелює з експе-риментальними даними, принаймні для використаної вибірки інгібіторів протеїнкінази СК2.

Ідентифікація нових інгібіторів протеїнкінази СК2. На основі розроблених методів та моделей створено технологію передбачення спорідненості низькомолекулярних органічних сполук до молекулярно-біологічних мішеней, яку реалізовано у вигляді лабораторного обчислювального комплексу ВПВС.

Із його застосуванням проведено скринінг in silico хімічної бази даних відділу комбінаторної хімії ІМБіГ НАН України, що нараховує близько 90 000 сполук. Із них віртуально відібрано 0,1% перспективних сполук. У цій вибірці біологічним скринінгом ідентифіковано низку нових інгібіторів протеїнкінази СК2, зокрема невідомий досі клас селективних інгібіторів - тетрагалогенізоіндолів. Кі найафіннішого із них складає 100нM. Швидкодія системи ВПВС у цьому експерименті склала 14 діб на 8-ми процесорному обчислювальному кластері відділу комбінаторної хімії ІМБіГ НАН України.

ВИСНОВКИ

На основі нових фізико-хімічних моделей створено оригінальні технології потокового РО передбачення спорідненості органічних сполук до молекулярно-біологічних мішеней in silico з поліпшеними прогностичними можливостями та обчислювально прості та ефективні алгоритми їхньої комп'ютерної реалізації.

1. Розроблено новий метод обчислення точкових атомних зарядів органічних сполук, що ґрунтується на методі релаксації електронегативностей за Кіргофом, та проведено його адекватну параметризацію. Доведено, що розроблений метод гарантує точність відтворення експериментальних величин дипольного та квадрупольного моментів органічних молекул із похибкою, що в середньому не перевищує 0,3 D та 0,2 B відповідно.

2. На основі розробленого методу обчислення атомних зарядів органічних молекул та силового поля MMFF94 створено нове силове поле YFF, яке забезпечує вищу точність оцінки електростатичної складової енергії міжмолекулярних взаємодій порівняно з відомими силовими полями, що використовуються у ВПВС.

3. Модифіковано метод GBSA розрахунку енергії десольватації органічних сполук. Одержані цим методом значення задовільно корелюють (r=0,84) з експериментальними величинами.

4. На основі функцій розподілу обертально-поступальних та внутрішніх ступенів вільності ліганду створено нову модель передбачення спорідненості органічних молекул до молекулярно-біологічних мішеней, яка враховує основні фізико-хімічні особливості просторової будови обох партнерів взаємодії. Розроблено обчислювально простий та ефективний алгоритм її комп'ютерної реалізації.

5. На основі розроблених моделей та алгоритмів створено лабораторний програмний комплекс для оцінки молекулярної спорідненості з урахуванням ентропійних ефектів у задачах РО ВПВС із залученням додаткового алгоритму молекулярного докінгу. З його використанням передбачено афінність органічних молекул-лігандів до активного центру протеїнкінази СК2 та ідентифіковано низку нових ефективних інгібіторів цього ферменту.

ПЕРЕЛІК НАУКОВИХ ПРАЦЬ, ОПУБЛІКОВАНИХ ЗА ТЕМОЮ ДИСЕРТАЦІЇ

1. Application of distribution function of rotation and translation degrees of freedom for CK2 inhibitors Ki estimation / O.Ya. Yakovenko, G.A. Golub, V.G. Bdzhola [et al.] // Ukrainica Bioorganica Acta. - 2006. - Vol.4, №2. - P.47-56.

Особистий внесок здобувача - розробка та програмна реалізація методу оцінки спорідненості лігандів до молекулярно-біологічних мішеней.

2. The new method of distribution integrals evaluations for high throughput virtual screening / O.Ya. Yakovenko, A.A. Oliferenko, G.A. Golub [et al.] // Ukrainica Bioorganica Acta. - 2007. - Vol.5, №1. - P.52-63.

Особистий внесок здобувача - розробка та програмна реалізація методів оцінки спорідненості з урахуванням ентропії та енергії десольватації.

3. Evaluation of 4,5,6,7-tetrahalogeno-1H-isoindole-1,3(2H)-diones as inhibitors of human protein kinase CK2 / A.G. Golub, O.Ya. Yakovenko, A.O. Prykhod'ko [et al.] // Biochim. Biophys. Acta. - 2008. - Vol.1784, №1. - Р.143-149.

Особистий внесок здобувача - проведення експериментів з віртуального скринінгу in silico.

4. Kirchhoff atomic charges fitted to multipole moments: implementation for a virtual screening system / O.Ya. Yakovenko, A.A. Oliferenko, V.G. Bdzhola [et al.] // J.Comput.Chem. - 2008. - Vol.29, №8. - P.1332-1343.

Особистий внесок здобувача - розробка та програмна реалізація методу релаксації електронегативностей за Кіргофом.

5. Implementation of orbital graph model for partial atom charges definition in virtual screening system : Conference thesis [“Inhibitors of Protein Kinases 2005”], (Warsaw, 25-29 June, 2005) / Polish Society for Cell Biology and the University of Warsaw, 2005 - P.126.

Особистий внесок здобувача - застосування методу релаксації електронегативностей за Кіргофом у лабораторній системі потокового віртуального скринінгу.

6. Development and improvement of receptor-based virtual screening approach : Abstract book of conference for young scientists, PhD students and students on molecular biology and genetics, ( Kyiv, 20-22 September 2007) / Inst. of Mol. Boil. and Gen.- K., 2007 - 206p.

Особистий внесок здобувача - застосування функцій розподілу молекулярних станів у лабораторній системі потокового віртуального скринінг.

7. Масштабируемое силовое поле для моделирования наноразмерных биологических структур Тезисы конференции [“Многомасштабное моделирование процессов и структур в нанотехнологиях”], (Москва, 14-12 марта 2008г.) / Министерство образования и науки Российской федерации, РНЦ «Курчатовский институт», Московский инженерно-физический институт, 2008 - С.226-227.

Особистий внесок здобувача - розробка та реалізація нового методу силового поля на основі методу релаксації електронегативностей Кіргофа.

АНОТАЦІЯ

Яковенко О.Я. Розробка технології передбачення спорідненості органічних сполук до молекулярно-біологічних мішеней для рецепторно-орієнтованого скринінгу in silico. - Рукопис.

Дисертація на здобуття наукового ступеня кандидата біологічних наук за спеціальністю 03.00.20 - біотехнологія. - Інститут молекулярної біології та генетики НАН України, Київ, 2008.

Мета дисертаційної роботи - створити на базі нових фізико-хімічних моделей технології потокового рецепторно-орієнтованого передбачення афінності органічних сполук до молекулярно-біологічних мішеней, включаючи алгоритми їхньої комп'ютерної реалізації. Розроблено та параметризовано новий емпіричний метод розрахунку розподілу електронної густини. На його основі створено нові модифікації методу GBSA для оцінки енергії десольватації та методу силового поля для оцінки енергії міжмолекулярних взаємодій. Ці методи застосовано для передбачення енергії взаємодії у модельній системі скринінгу in silico, що складається з ліганду, молекулярно-біологічної мішені та води. Для оцінки спорідненості з урахуванням ентропії розроблено та перевірено метод, що ґрунтується на функціях розподілу молекулярних станів. Із застосуванням цих технологій ідентифіковано ряд нових ефективних інгібіторів протеїнкінази СК2.

Ключові слова: спорідненість, енергія взаємодії, функції розподілу станів, ліганд, молекулярно-біологічна мішень, віртуальний скринінг.

АННОТАЦИЯ

Яковенко А.Я. Разработка технологий предсказания сродства органических соединений к молекулярно-биологическим мишеням для рецепторно-ориентированного скрининга in silico. - Рукопись.

Диссертация на соискание ученой ступени кандидата биологических наук по специальности 03.00.20 - биотехнология. - Институт молекулярной биологии и генетики НАН Украины, Киев, 2008.

Цель диссертационной роботы - создание на базе новых физико-химических моделей технологии потокового рецепторно-ориентированного предсказания афинности органических соединений к молекулярно-биологическим мишеням, включая алгоритмы их компьютерной реализации. Разработан и параметризирован новый эмпирический метод расчета распределения электронной плотности. На его основе предложены новые модификации метода GBSA для оценки энергии десольватации и метода силового поля для оценки энергии межмолекулярных взаимодействий. Эти методы применены для предсказания энергии взаимодействий в модельной системе высокопродуктивного скрининга in silico, которая состоит из лиганда, молекулярно-биологической мишени и воды. Для оценки сродства с учетом энтропии разработан и проверен метод, основанный на функциях распределения по молекулярным состояниям. С применением этих технологий идентифицирован ряд новых эффективных ингибиторов протеинкиназы СК2.