Projektowanie leków

Modelowanie w ICM
Mod. biomolekularne
Projektowanie leków
Numeryka
Wizualizacja
Inna działalność

Spis treści

Wstęp

Ciągły wzrost mocy obliczeniowej komputerów, postępy bioinformatyki oraz rozwój strukturalnych baz danych białek i małych molekuł sprawiają, że komputerowe wspomaganie projektowania leków (Computer Aided Drug Design) staje się coraz wydajniejszym i coraz częściej stosowanym narzędziem.

Współczesne metody, choć nie pozwalają jeszcze na pewne i wiarygodne przewidzenie ostatecznej struktury działających leków, umożliwiają wybór obiecujących kandydatów spośród milionów molekuł dostępnych w bazach danych oraz ich optymalizację pod kątem oddziaływania z receptorem. Pozwala to obniżyć koszty i czas trwania procesu projektowania nowego leku.

Schemat wykorzystania poszczególnych metod komputerowego projektowania leków

Rysunek przedstawia schemat wykorzystania podstawowych metod komputerowych we wstępnych etapach projektowania leków. To, czy projektowanie zakończy się sukcesem w dużej mierze zależy od jakości danych wejściowych. Etapy zaznaczone kolorem różowym przyczyniają się do wzrostu ryzyka niepowodzenia. Konieczność skonstruowania modelu receptora (np. poprzez modelowanie homologiczne) w przypadku braku jego struktury krystalograficznej, bądź nieznajomość miejsca wiążącego wymagają wykorzystania dodatkowych założeń i przybliżeń, przez co zwiększają prawdopodobieństwo uzyskania nieprawidłowych wyników. Generując ligandy de novo częściej otrzymuje się struktury nie będące lekami, niż wtedy, gdy ogranicza się do modyfikacji znanego wcześniej związku wiodącego (lead compound).

Zależnie od wiedzy, którą dysponuje się przystępując do projektowania leków wyróżnić można dwa podstawowe nurty:

  • projektowanie w oparciu o strukturę znanych ligandów (ligand based design), które obejmuje takie etapy jak budowa farmakoforu, analiza QSAR, budowa pseudo-receptora
  • projektowanie w oparciu o znaną strukturę receptora (receptor based design), które obejmuje analizę miejsca wiążącego, analizę zmienności konformacyjnej receptora, konstrukcję więzów farmakoforycznych

Podział ten dotyczy raczej początkowych etapów projektowania leków, zmierzających do stworzenia modelu (farmakoforu lub modelu miejsca wiążącego) warunkującego właściwości przestrzenne i fizykochemiczne projektowanych molekuł. Kolejne kroki takie jak: dokowanie, przeszukiwanie baz danych, generowanie nowych ligandów mogą w zasadzie przebiegać w podobny sposób niezależnie od etapów początkowych.

W optymalnej sytuacji znana jest zarówno grupa ligandów o zbadanej aktywności jak i kilka struktur krystalograficznych receptora. W takim przypadku projektowanie oparte na znajomości struktury receptora zdecydowanie warto jest uzupełnić o analizę farmakoforu i model QSAR.


Budowanie farmakoforu

Farmakofor jest modelem opisującym relacje przestrzenne między elementami wspólnymi dla ligandów oddziałujących z danym receptorem [1,2]. Tworzenie farmakoforu jest kluczowym etapem projektowania leków w przypadku, gdy znana jest struktura co najmniej kilku aktywnych ligandów, a brakuje danych dotyczących budowy receptora.

Przykładowy farmakofor: wspólne elementy strukturalne i ich wzajemne ułożenie dla czterech aktywnych ligandów

Farmakofor może być wykorzystany w wielu aspektach projektowania leków:

  • do przeszukiwania baz danych
  • do generowania ligandów de novo
  • do optymalizacji molekuł wiodących
  • do znalezienia względnych pozycji molekuł badanych metodami 3D-QSAR

Obecnie istnieje wiele programów komputerowych służących do generowania farmakoforów na podstawie zbioru ligandów. Mimo, iż wykorzystują one różnorodne algorytmy wszystkie dzielą zadanie na podobne etapy.

Etapy wstępne

Przygotowanie danych

Do budowy farmakoforu niezbędna jest znajomość struktury chemicznej zbioru aktywnych ligandów danego receptora. Najlepiej jeśli ligandy te stanowią maksymalnie aktywną i różnorodną grupę. Niezbędnym warunkiem jest jednak aby wszystkie oddziaływały z tym samym miejscem wiążącym. Niektóre programy (DISCO, CLEW) umożliwiają również wykorzystanie informacji o strukturach nieaktywnych ligandów, co pozwala wyodrębnić cechy których należy unikać. Maksymalna liczba wykorzystywanych ligandów jest różna dla poszczególnych programów i zależy też od rozmiarów ligandów (w szczególności od liczby wiązań mogących podlegać rotacji). Może ona wynosić od kilku do kilkuset cząsteczek.

Przygotowując ligandy należy zwrócić uwagę na prawidłowe przypisanie typów i ładunków (jeśli używany program tego wymaga) poszczególnym atomom. Od tego zależeć będzie jak zostaną one rozpoznane w następnych etapach analizy. Bardzo istotne są również stany protonacyjne grup miareczkowanych. To, czy dana grupa zostanie rozpoznana jako donor bądź akceptor wiązania wodorowego może zależeć od obecności protonu.

Przeszukiwanie przestrzeni konformacyjnej

Z reguły ligandy są giętkimi molekułami i konformacja w jakiej łączą się z receptorem nie jest znana. Budowanie farmakoforu wymaga zatem zbadania przestrzeni konformacyjnej analizowanych cząsteczek. Część programów na wstępie generuje różne konformacje i zapisuje je do wykorzystania w kolejnych etapach. Inne natomiast uwzględniają giętkość molekuł w trakcie tworzenia farmakoforu.

Analiza całej przestrzeni konformacyjnej jest niemożliwa ze względu na jej rozmiar. Ogranicza się ją zatem do konformacji o najniższych energiach potencjalnych bądź też dokonuje się klasteryzacji i rozważa potem tylko konformacje reprezentatywne dla utworzonych grup.

Wybór elementów i zakodowanie struktury molekuł

Farmakofor zbudowany jest jedynie z elementów odgrywających istotną rolę w oddziaływaniach międzycząsteczkowych. Mogą to być atomy szczególnych typów (np. atom siarki), określone grupy chemiczne (np. grupa karboksylowa, aminowa, pierścień fenylowy), bądź też grupy o szczególnych funkcjach (np. donor lub akceptor wiązania wodorowego, grupa o charakterze zasadowym lub kwasowym, grupa hydrofobowa).

Z każdej molekuły wyodrębniane są tylko te elementy, które mogą służyć do budowy farmakoforu. W dalszych etapach są one reprezentowane w sposób symboliczny, najczęściej jako punkt skojarzony ze środkiem ciężkości atomów tworzących dany element. Każda rozważana konformacja jest następnie zapisywana w uproszczonej postaci. Z reguły stosowane są takie formy zapisu jak: zbiór punktów w przestrzeni trójwymiarowej, graf (węzły reprezentują poszczególne elementy, a połączenia ich wzajemne odległości) lub macierz odległości między wszystkimi parami elementów.

Generowanie farmakoforów

Analiza danych przygotowanych w etapach wstępnych sprowadza się najczęściej do znalezienia maksymalnego farmakoforu wspólnego dla wszystkich badanych molekuł. Jest to oparte na założeniu, że wszystkie ligandy oddziałując z tym samym miejscem wiążącym muszą dzielić jeden, wspólny farmakofor. Założenie to nie zawsze się sprawdza. Istnieją np. wysoce aktywne ligandy, nie zawierające któregoś z elementów farmakoforu.

Innym podejściem jest poszukiwanie najmniejszego zbioru cech farmakoforycznych (elementów i ustalonych odległości między nimi), który ma co najmniej zadane m cech wspólnych z każdym badanym ligandem. Nie wymagane jest zatem, aby każda molekuła zawierała dokładnie wszystkie elementy farmakoforu.

Do analizy danych wykorzystywane są różne metody, zależne od sposobu zakodowania struktury molekuł. Wynikiem ich pracy jest z reguły kilka propozycji farmakoforów, z których po odpowiedniej ocenie wybierany jest ostateczny.

Proste przeszukiwanie

polega na badaniu możliwych kombinacji pasujących do siebie cech (rozumianych jako para elementów i odległości między nimi). We wstępnym etapie generowane są wszystkie możliwe farmakofory zawierające jedynie po dwa elementy (w ustalonej odległości od siebie) wspólne dla wszystkich ligandów. Następnie rozważane są farmakofory trójelementowe poprzez dodanie trzeciej cechy do znalezionych wcześniej farmakoforów dwuelementowych. Proces ten jest powtarzany do momentu znalezienia farmakoforu o maksymalnej liczbie cech wspólnych dla wszystkich ligandów.

Algorytmy genetyczne

We wstępnym etapie generowane są "chromosomy" zawierające informacje o aktualnym położeniu elementów danej molekuły w danej konformacji. Jedna z molekuł, z reguły ta o najmniejszej liczbie elementów wybierana jest jako wzorcowa. Kolejne "pokolenia chromosomów" generowane są poprzez losowe zmiany (translacje, rotacje) położeń pozostałych molekuł oraz ewentualne zmiany ich konformacji (jeśli zmienność konformacyjna nie została uwzględniona we wstępnych etapach). Przeżywalność poszczególnych "chromosomów" zależy od poprawności dopasowania zakodowanych w nich elementów do struktury molekuły wzorcowej. Efektem "ewolucji" jest grupa chromosomów reprezentująca najlepsze dopasowanie do molekuły wzorcowej. Na podstawie analizy wspólnych elementów generowany jest farmakofor.

Znajdowanie największego podgrafu

jest wykorzystywane przez programy zapisujące struktury molekuł w postaci grafów.

Wynikiem działania opisanych algorytmów są wstępne propozycje farmakoforów. Z nich wybierany jest ostateczny model na podstawie różnorodnych kryteriów. Znaczenie ma liczba elementów w proponowanych farmakoforach. Generalnie im więcej tym lepiej. Nie wszystkie jednak elementy są równie ważne - z reguły im rzadziej spotykany jest dany element tym jest on ważniejszy. Np. centra naładowane są ważniejsze niż grupy hydrofobowe.

W przypadku metod dopuszczających aby elementy farmakoforu nie były wspólne dla wszystkich ligandów z badanej grupy, znaczenie ma liczba molekuł dla których dany farmakofor jest częścią wspólną. Np. farmakofor zbudowany z trzech elementów wspólnych dla wszystkich ligandów może okazać się bardziej istotny niż farmakofor czteroelementowy reprezentatywny dla wszystkich z wyjątkiem jednego ligandów.

Wybrane programy


QSAR

QSAR (Quantitative Structure - Activity Relationships) [8] jest to metoda polegająca na analizie zależności między aktywnością biologiczną grupy ligandów danego receptora, a ich cechami fizykochemicznymi. Dzięki niej możliwe jest zbadanie które właściwości tych ligandów mają największy wpływ na aktywność oraz sformułowanie matematycznego modelu tej zależności - równania QSAR. Znajomość tego równania pozwala:

  • wyciągnąć wnioski co do mechanizmów interakcji ligandów z receptorem
  • ocenić aktywność niezbadanych jeszcze związków
  • zaproponować najkorzystniejsze z punktu widzenia aktywności modyfikacje
  • uporządkować i pogrupować ligandy pod względem ich aktywności oraz cech fizykochemicznych

Aby możliwe było przeprowadzenie analizy statystycznej i sformułowanie równania niezbędne jest wyrażenie parametrów fizykochemicznych badanych związków w postaci liczbowej. Służą do tego tzw. deskryptory.

Deskryptory

Aktywność biologiczna związków chemicznych zależy od szeregu ich cech takich jak: rodzaj podstawników, rozmiar cząsteczki, ładunek, polaryzowalność czy rozpuszczalność w wodzie i tłuszczach. Deskryptory są funkcjami pozwalającymi przypisać tym cechom wartości liczbowe. Ze względu na charakter właściwości, które opisują można podzielić je na różne grupy:

  • właściwości strukturalne, np: masa molowa, objętość molekularna, pole powierzchni, momenty bezwładności, liczba wiązań pojedynczych (rotujących), parametry steryczne Tafta E_S
  • właściwości elektronowe, np: ładunek całkowity, moment dipolowy, polaryzowalność, energia HOMO, energia LUMO, stałe Hammeta, współczynniki indukcyjne i rezonansowe
  • właściwości termodynamiczne, np: współczynnik podziału między wodę i oktanol \log P, hydrofobowość podstawników (\pi_X), refraktywność molowa (MR)
  • właściwości przestrzenne pól - patrz 3D QSAR

Oprócz podanych powyżej podstawowych deskryptorów stosowanych jest jeszcze wiele innych np. opartych na analizie konformacyjnej molekuł bądź zastosowaniu teorii grafów w odniesieniu do połączeń międzyatomowych. Spośród klasycznych deskryptorów niektóre mogą wymagać wyjaśnienia:

Współczynnik podziału \log P

jest zdefiniowany jako logarytm stosunku stężeń równowagowych danej substancji w wodzie i oktanolu:

P=\frac{[lek]_{oktanol}}{[lek]_{woda}}

Jest on miarą hydrofobowości molekuły i koreluje z jej zdolnością do przenikania przez błony biologiczne. Ma duże znaczenie np. w przypadku leków działających na ośrodkowy układ nerwowy, dla których istotnym parametrem jest zdolność do przechodzenia przez barierę krew-mózg.

Chociaż jest on zdefiniowany w oparciu o dane doświadczalne programy do analizy QSAR z reguły zawierają algorytmy pozwalające oszacować jego wartości na podstawie budowy chemicznej danego związku. Parametr \log P uzyskany w wyniku takich obliczeń nosi nazwę CLOGP.

Współczynniki \pi

opisują hydrofobowość poszczególnych grup funkcyjnych. Są zdefiniowane następująco:

\pi_X=\frac{\log P_{RX}}{\log P_{RH}}

gdzie P_{RH} - wartość P dla wzorcowej molekuły z atomem wodoru w miejscu grupy funkcyjnej, P_{RX} - wartość P dla molekuły z grupą funkcyjną X. W analizie QSAR z reguły zakłada się ich addytywność w obrębie danej molekuły i stosuje \sum_X \pi_X jako jeden deskryptor.

Refraktywność molowa

dana wzorem:

MR=\frac{n^2-1}{n^2+1}\frac{MW}{d}

gdzie n - współczynnik załamania, MW - masa molowa, d - gęstość molowa, odzwierciedla kombinację właściwości elektronowych (współczynnik załamania) oraz strukturalnych (objętości cząsteczki).

Parametr steryczny Tafta

jest zdefiniowany tylko dla niektórych podstawników i odzwierciedla właściwości otoczenia molekularnego R (przede wszystkim jego rozmiar), do którego dany podstawnik jest przyłączony. Oparty jest na analizie wpływu otoczenia R na stałą reakcji k:

CADD-Taft.png

i przedstawia się wzorem: E_S=\frac{\log k_{R CO_2CH_3}}{\log k_{CH_3 CO_2 CH_3}}

gdzie k_{CH_3 CO_2 CH_3} - stała reakcji dla molekuły wzorcowej, k_{R CO_2CH_3} - stała reakcji dla molekuły zawierającej R.

Stała Hammeta

\sigma jest miarą wpływu podstawnika R na gęstość elektronową w pierścieniu aromatycznym. Pozwala rozróżnić te podstawniki, które "oddają" elektrony do pierścienia od tych, które "ściągają" gęstość elektronową na siebie. Oparta jest na badaniu następujących reakcji dysocjacji

CADD-Hammet.png

Przedstawia się wzorem:

\sigma_X=\frac{\log k_{R=X}}{\log k_{R=H}}

Wartość stałej Hammeta zależy również od pozycji podstawnika względem grupy dysocjującej. Stąd też dla wybranego podstawnika rozważać można trzy stałe - dla położeń orto, meta i para.

Analiza deskryptorów

Głównym zadaniem analizy QSAR jest wygenerowanie równania opisującego zależność aktywności biologicznej ligandów od wartości poszczególnych deskryptorów. Najczęściej poszukuje się równania liniowego postaci:

\log A = a_1 D_1 + a_2 D_2 +a_3 D_3 + ...

gdzie: A oznacza aktywność biologiczną z reguły przyjmowaną jako odwrotność stężenia C ligandu powodującego określony efekt (IC_{50}, LD_{50} itp), a_i wagę z jaką uwzględniane są wartości D_i (lub czasem kwadraty ich wartości) i-tego deskryptora. W klasycznej analizie QSAR do wyznaczenia wartości współczynników a_i wykorzystuje się metodę wielokrotnej regresji.

Należy pamiętać jednak o tym, że w metodzie wielokrotnej regresji liczba szukanych współczynników (a więc i liczba wykorzystywanych deskryptorów) musi być znacznie mniejsza niż liczba punktów (czyli molekuł o znanej aktywności) do których dopasowujemy dane. W praktyce jeden współczynnik powinien przypadać na co najmniej pięć molekuł. Dodatkowo wartości wielu deskryptorów mogą być ze sobą skorelowane, zatem ich jednoczesne stosowanie nie ma sensu.

Problemy te można częściowo ominąć dzięki metodzie analizy czynników głównych (Principal Component Analysis). Polega ona na znalezieniu wektorów własnych macierzy kowariancji deskryptorów. Następnie jako "nowe" deskryptory, będące liniowymi kombinacjami oryginalnych deskryptorów, wybiera się te, które odpowiadają kilku największym wartościom własnym macierzy kowariancji. Tak utworzone deskryptory nie są wzajemnie skorelowane oraz najlepiej opisują zmiany aktywności.

Innym podejściem jest wykorzystanie algorytmów klasteryzacji. Dzięki nim można pogrupować ligandy względem podobnych wartości w przestrzeni deskryptorów. Pozwala to ujawnić właściwości sprzyjające korzystnemu oddziaływaniu z receptorem.

3D-QSAR

Metody oparte na analizie opisanych deskryptorów nie pozwalają na uwzględnienie pełnej informacji jaką niesie ze sobą znajomość struktury przestrzennej cząsteczki. Informację tę wykorzystują natomiast metody bazujące na analizie właściwości pól molekularnych - CoMFA (Comparative Molecular Field Analysis).

Ich generowanie przebiega następująco: cząsteczkę o wybranej konformacji otacza się siatką sześcienną. Węzłom siatki przypisuje się wartości energii oddziaływania umieszczanej w nich cząstki próbnej (np. ładunku punktowego lub grupy CH{_3}) z badaną molekułą. Uzyskuje się w ten sposób opis przestrzennych właściwości pól (w tym przypadku odpowiednio elektrostatycznego lub sterycznego, odpowiadającego oddziaływaniom van der Waalsa). Stosując różnorodne cząstki próbne można generować też inne pola np: wiązań wodorowych, oddziaływań hydrofobowych, lipofilowe.

Wartości danego pola w każdym z węzłów traktuje się następnie jako niezależne deskryptory. Oczywiście deskryptory takie są ze sobą silnie skorelowane oraz jest ich znacznie więcej niż badanych molekuł. Do ich analizy wykorzystuje się metodę częściowych najmniejszych kwadratów (Partial Least Squares -PLS) lub algorytmy genetyczne. Wynikiem są takie kombinacje liniowe deskryptorów, które maksymalizują stopień objaśnienia zmienności aktywności.

CADD-Comfa.png

Wynik analizy CoMFA - przykładowa molekuła na tle pól elektrostatycznego i sterycznego. Kolory pól oznaczają: czerwony niebieski - obszary, w których preferowane są podstawniki dodatnio/ujemnie naładowane, żółtyielony - obszary, w których preferowane są podstawniki o dużym/małym rozmiarze.

Metoda CoMFA wymaga, aby wszystkie analizowane molekuły były układane względem siatki w takich pozycjach i orientacjach jakie przyjmują w miejscu wiążącym receptora. Osiąga się to poprzez dopasowanie każdej z nich do wspólnego farmakoforu.

Wybrane programy

Znajdowanie miejsca wiążącego

Znajomość topologii miejsca wiążącego receptora ma istotne znaczenie w procesie projektowania leków. Nie zawsze jednak dysponuje się strukturą krystalograficzną receptora z przyłączonym ligandem, na podstawie której można jednoznacznie zdefiniować miejsce wiążące. Najczęściej potrzeba poszukiwania miejsca wiążącego zachodzi w następujących przypadkach:

  • znana jest jedynie struktura krystalograficzna natywnej formy receptora
  • znana jest jedynie sekwencja aminokwasowa receptora
  • poszukuje się alternatywnego miejsca wiążącego

W takich sytuacjach można skorzystać z jednej z opisanych poniżej metod [10]. Należy jednak pamiętać, że nawet poprawne wytypowanie miejsca wiążącego nie gwarantuje sukcesu - niekiedy w wyniku związania liganda receptor podlega dużym zmianom konformacyjnym, które trudno przewidzieć bez wykonywania dokładnych symulacji.

Metody oparte na homologii

W sytuacji, gdy poszukuje się miejsca wiążącego dla naturalnych ligandów danego receptora istnieje duża szansa, że jest ono podobne do innych, znanych już miejsc wiążących w pokrewnych białkach. Pierwszym krokiem jest porównanie sekwencji badanego receptora z sekwencjami tych białek. Układ aminokwasów w miejscach wiążących pozostaje często zachowany w procesie ewolucji, stąd też znalezienie obszarów o małej zmienności sekwencyjnej może sugerować, że to właśnie one tworzą miejsce wiążące.

W przypadku, gdy sekwencja badanego receptora okaże się unikalna, można jeszcze liczyć na podobieństwo strukturalne do białek z dobrze opisanym miejscem wiążącym. Wykorzystać można bazę danych białek podobnych strukturalnie. Na przykład baza FSSP [11,12] (Families of structurally similar proteins) jest dostępna bezpłatnie dla użytkowników akademickich. Do wyszukiwania białek o podobnych strukturach służyć może program Dali [12].

Metody oparte na modelu solwatacyjnym

Ta grupa metod służy do odnajdywania zagłębień w powierzchni molekularnej białka. Generalnie działają one według następującego algorytmu: cała powierzchnia białka zostaje "oblepiona" sferami o rozmiarze podobnym do rozmiaru cząsteczki wody. Te sfery, w których otoczeniu o określonym promieniu jest mniej niż zadane N atomów białka, zostają usunięte. Następnie białko wraz z pozostałymi sferami zostaje oblepione kolejną warstwą sfer i powtarzana jest operacja usuwania. Cykl ten wykonywany jest do momentu, gdy liczba sfer znajdujących się przy powierzchni białka przestanie się zmieniać.

CADD-Solvation.png

Kolejne etapy metody solwatacyjnej (zmodyfikowane na podstawie [10]).

Sfery, które nie zostały usunięte tworzą grupy wypełniające zagłębienia w powierzchni molekularnej. Grupy te służą jako punkty odniesienia do wyboru miejsc wiążących.

Metody oparte na siatkach

Metody te są dość podobne do opisanych powyżej. Struktura białka otaczana jest punktami leżącymi w węzłach sieci sześciennej. Spośród tych punktów eliminowane są te, w których otoczeniu (sześciennym obszarze o zadanym boku) nie ma ani jednego atomu białka. Dzięki temu pozostawiane są tylko punkty leżące blisko powierzchni białka. W kolejnym kroku usuwane są te punkty, w których sąsiedztwie (o innym niż poprzednio rozmiarze) znajduje się mniej niż określona liczba atomów białka. Następnie punkty łączone są w grupy. Te z nich, które zawierają zbyt mało punktów są eliminowane, a resztę pozostawia się bez zmian, lub łączy w większe grupy jeśli odległość między nimi jest mniejsza od zadanego progu.

Metody oparte na siatkach lepiej niż solwatacyjne (choć w dużej mierze jest to kwestia doboru parametrów) nadają się do znajdowania głębokich i wąskich szczelin. Ich wadą jest natomiast zależność od położenia i orientacji białka względem siatki.

Algorytmy służące do znajdowania miejsc wiążących dostarczają z reguły kilku rozwiązań. To, które z nich wybrać zależy od intuicji badacza. Pomocna może być również analiza rezultatów otrzymanych różnymi metodami. Niekiedy wykonuje się też dokowania naturalnych ligandów do każdego z proponowanych miejsc wiążących i na podstawie ich wyników wybiera się najlepsze.

Wybrane programy


Przeszukiwanie baz danych małych cząsteczek

Przeszukiwanie molekularnych baz danych [13] umożliwia wybranie spośród milionów zgromadzonych w nich molekuł tych, które spełniają zadane kryteria. Pozwala to na ograniczenie liczby cząsteczek analizowanych w kolejnych etapach bardziej kosztowymi metodami, np. za pomocą dokowania.

Początkowo przeszukiwanie umożliwiało znalezienie molekuł zawierających jedynie zadane elementy strukturalne - deskryptory (2D search). W miarę wzrostu mocy obliczeniowych w kryteriach wyszukiwania zaczęto uwzględniać relacje przestrzenne pomiędzy poszczególnymi deskryptorami (3D search), możliwość zmian konformacyjnych molekuł (3D flexible search) oraz więzy wynikające ze struktury miejsca wiążącego (receptor site search).

Przeszukiwanie dwuwymiarowe

Ten rodzaj przeszukiwania pozwala znaleźć molekuły zawierające określone elementy strukturalne: grupy chemiczne, rodzaje wiązań lub atomów, donory lub akceptory wiązań wodorowych. Wzajemne położenia tych elementów, rozmiary molekuł ani ich konformacje nie są brane pod uwagę.

Przykładowe zapytanie 2D i jego wynik

Na etapie generowania bazy danych cechy każdej molekuły kodowane są w postaci wektora binarnego - ciągu jedynek i zer opisujących odpowiednio obecność lub brak danej cechy. Wektory takie mogą składać się z kilkuset pozycji.

Tworzone przez użytkownika zapytanie również tłumaczone jest na ciąg zer i jedynek. Przeszukiwanie bazy danych polega na porównywaniu wektorów binarnych, dzięki czemu jest szybko i wydajnie realizowane przez komputer.

Gdy zapytanie ma rozbudowaną postać często pożądanym wynikiem jest znalezienie nie tylko tych molekuł, które zawierają dokładnie wszystkie jego elementy, ale również molekuł wykazujących jedynie pewien stopień podobieństwa (similarity search). Podobieństwo takie można obliczać na kilka sposobów. Najczęściej korzysta się ze współczynnika Tanimoto (Tanimoto coefficient), który przyjmuje wartości od 0 gdy porównywane struktury nie mają wspólnych deskryptorów do 1, gdy wszystkie deskryptory są jednakowe. Współczynnik ten przedstawia się wzorem:

S_T=\frac{N_{QT}}{N_Q+N_T-N_{QT}}

gdzie: N_Q - liczba deskryptorów zapytania, N_T - liczba deskryptorów badanej molekuły, N_{QT} - liczba wspólnych deskryptorów. W wielu programach domyślna wartość współczynnika Tanimoto jest ustawiona na 0.85.

Przeszukiwanie trójwymiarowe

Ten rodzaj przeszukiwania pozwala uwzględnić również relacje przestrzenne między poszczególnymi deskryptorami. Zapytanie ma postać farmakoforu.

Przykładowe zapytanie 3D i jego wynik

Przeszukiwanie trójwymiarowe wymaga znajomości struktur przestrzennych molekuł. Zwykle jednak bazy danych zawierają struktury dwuwymiarowe lub jedynie opis połączeń międzyatomowych np. zakodowany w systemie SMILES. W tej sytuacji można skorzystać z programów automatycznie generujących struktury trójwymiarowe. Do najbardziej popularnych należą: CORINA, CONCORD, Converter, ChemX, Molgeo.

Ich tworzenie odbywa się na podstawie algorytmów opartych na wiedzy o najczęstszych konformacjach określonych grup atomów - wykorzystywanie do tego celu mechaniki molekularnej byłoby zbyt czasochłonne. Cząsteczka dzielona jest na fragmenty, których struktura odczytywana jest z odpowiednich tablic. Następnie fragmenty są łączone, przy czym analizuje się ich różne możliwe ustawienia tak, aby uniknąć zbyt bliskich kontaktów między atomami.

Metody tego typu pozwalają dla małych molekuł wygenerować konformacje zbliżone do doświadczalnych. Wyniki są gorsze dla cząsteczek zawierających duże nienasycone pierścienie bądź dla kompleksów metali.

Należy jednak brać pod uwagę fakt, że w rzeczywistości cząsteczki nieustannie zmieniają konformację, a zwłaszcza dzieje się to podczas kontaktu z receptorem. Przeszukiwanie baz danych w oparciu o sztywne, trójwymiarowe modele cząsteczek (zarówno w konformacjach ustalonych doświadczalnie jak i wygenerowanych komputerowo) prowadzi do pominięcia wielu dobrych kandydatów.

Ominąć ten problem można wykorzystując algorytmy uwzględniające zmienność konformacyjną badanych molekuł (3D flexible search). Generowane są różne konformacje dzięki dopuszczeniu rotacji wokół wybranych wiązań pojedynczych i dla każdej z nich sprawdzana jest zgodność z zapytaniem. Choć daje dużo lepsze wyniki niż przeszukiwanie baz sztywnych molekuł jest to metoda bardzo czasochłonna. Wykorzystuje ją np. program UNITY.

Przeszukiwanie w oparciu o strukturę miejsca wiążącego

W przypadku, gdy znana jest struktura receptora, a nieznany jest farmakofor ligandów możliwe jest skonstruowanie zapytania w oparciu o cechy miejsca wiążącego.

Zapytanie takie składa się z definicji obszarów, w których powinny znaleźć się elementy liganda o określonych cechach, np. miejsca akceptorowe lub donorowe dla wiązań wodorowych, grupy o charakterze hydrofobowym, atomy określonego typu itp. Istotne jest również zdefiniowanie obszarów, w których ligand nie powinien się znaleźć.

CADD-ReceptorSearch.png

Przykładowe zapytanie. Po lewej stronie: struktura molekularna miejsca wiążącego służąca jako punkt odniesienia do skonstruowania zapytania (na zielono - ligand ze struktury krystalograficznej). Po prawej stronie: tak zapytanie widziane jest przez program przeszukujący. Fioletowe kule - obszary wzbronione dla atomów liganda, czerwona kula - obszar, w którym powinna się znaleźć grupa hydrofobowa, obszary jasno fioletowe i zielone - miejsca, w których powinny się znaleźć odpowiednio donory bądź akceptory wiązań wodorowych

W przypadku tego typu zapytań nie należy oczekiwać aby ligand zawierał wszystkie zdefiniowane cechy jednocześnie. Spowodowałoby to odrzucenie większości interesujących molekuł. Warto skorzystać z możliwości wskazania minimalnej oraz maksymalnej liczby spełnionych jednocześnie kryteriów (partial match). Opcja taka istnieje np. w programie UNITY.

Przeszukiwanie w oparciu o strukturę receptora wymaga uwzględnienia zmienności konformacyjnej ligandów jest zatem metodą bardzo czasochłonną w porównaniu z innymi rodzajami przeszukiwania. Należy stosować je zatem raczej do podzbioru bazy danych, wyodrębnionego dzięki wcześniejszym prostszym wyszukiwaniom. Wynik nie zawiera ponadto informacji o różnych możliwych ułożeniach poszczególnych ligandów w miejscu wiążącym ani oceny jakości wiązania - informacje te dostarcza dokowanie.

Wybrane programy i bazy danych

Zbiory baz danych:

Konwertery 2D \to 3D

Zaawansowane programy do zaawansowanego przeszukiwania 3D

Dokowanie ligandów do receptora

Celem dokowania jest znalezienie sposobu ułożenia liganda w miejscu wiążącym receptora (binding mode) oraz ocena siły wiązania[15]. Z reguły dla wybranego receptora analizie poddaje się wiele (do kilkunastu tysięcy) ligandów, stąd też algorytmy dokowania są oparte na kompromisie pomiędzy szybkością, a dokładnością oceny.

Przygotowanie układu

Do przeprowadzenia dokowania niezbędna jest znajomość pełnoatomowej struktury receptora oraz posiadanie biblioteki potencjalnych ligandów. Potrzebne jest również wskazanie miejsca wiążącego - dokowanie do całego receptora "na raz" jest pozbawione sensu.

Samo miejsce wiążące warto wnikliwie przeanalizować. Należy zwrócić uwagę na miejsca donorowe i akceptorowe dla wiązań wodorowych oraz aminokwasy miareczkowalne - poprawne działanie programów zależy od ustawienia właściwych stanów protonacyjnych tych miejsc. Istotne jest również prawidłowe ustawienie atomów wodoru dodanych do struktury krystalograficznej - kierunek, pod którym przyłączony jest proton do atomu-donora warunkuje więzy geometryczne dla potencjalnego wiązania wodorowego. W centrach aktywnych enzymów często występują dwudodatnie jony. Jeśli podejmiemy decyzję o ich pozostawieniu w procesie dokowania, należy upewnić się, że program przypisze im właściwy ładunek (często wymaga to wczytania dodatkowych plików z topologiami).

Odrębnym zagadnieniem jest obecność cząsteczek wody w miejscu wiążącym - często są one niezbędne do "mostkowania" wiązań wodorowych. Gdy dana cząsteczka jest np. obecna w kilku strukturach krystalograficznych receptora dobrze jest pozostawić ją jako element miejsca wiążącego. Niektóre programy np. FlexX posiadają funkcję pozwalającą automatycznie umieszczać cząsteczki wody w miejscu aktywnym w trakcie dokowania i oceniać ich przydatność.

Gdy posiadamy wiedzę o kluczowych residuach w miejscu wiążącym warto jest zdefiniować więzy farmakoforyczne. Pozwalają one wskazać miejsca, w które musi zostać wpasowany określony element liganda (np. akceptor wiązania wodorowego lub centrum hydrofobowe). Taką możliwość daje np. program FlexX-Pharm. Należy pamiętać, że gdy definiujemy kilka więzów w jednym miejscu wiążącym nie należy raczej żądać ich jednoczesnego spełnienia - lepiej jest dopuścić częściowe spełnienie (partial match) np. od dwóch do trzech z pięciu zadanych. Korzystniej jest wyszukać więcej różnorodnych ligandów i poddać je analizie, niż liczyć na znalezienie kilku "idealnych" molekuł.

Przygotowując ligandy niezbędne jest upewnienie się, że posiadają one prawidłowe typy atomów oraz, że atomom tym przypisano właściwe ładunki.

Algorytmy dokowania

Proces wiązania liganda z receptorem prowadzi do zmian konformacyjnych obu molekuł. Uwzględnienie pełnej zmienności konformacyjnej jest praktycznie niemożliwe ze względu na ogromny koszt numeryczny. W najprostszym przybliżeniu efekt ten jest zaniedbywany, standardem jest uwzględnianie zmian konformacyjnych samych ligandów, obecnie intensywnie rozwijane są natomiast metody uwzględniające także giętkość receptora.

Dokowanie sztywnych ligandów - przykładowe algorytmy [14]

Metoda grafów

Polega na zbudowaniu grafu, którego węzłami są kombinacje par pasujących do siebie cech ligandu i receptora (np. wnęki hydrofobowej receptora z grupą hydrofobową ligandu). W przypadku, gdy odległość (geometryczna) między miejscami zawierającymi dane cechy jest taka sama w ligandzie, jak i w receptorze odpowiadające im węzły grafu są łączone. Najlepsze dopasowanie liganda i receptora odpowiada znalezieniu największego w pełni połączonego podgrafu (clique) wyjściowego grafu.

CADD-Graph.png

Metoda grafów - duże litery oznaczają cechy miejsca wiążącego (np. B - obszar hydrofobowy, A1 i A2 miejsca donorowe dla wiązań wodorowych), małe - cechy ligandu (b - grupa hydrofobowa, a1, a2, a3 akceptory wiązania wodorowego). Np. A1a2 i Bb stanowią pary komplementarnych cech. Po prawej stronie graf zbudowany dla tego przypadku. Np. odległość B - A2 jest taka sama jak b - a2, stąd połączenie węzłów Bb i A2a2.

Metoda wykorzystująca pole siłowe

Polega na wypełnieniu obszaru miejsca aktywnego sferami stycznymi do powierzchni molekularnej, a następnie próbie dopasowania położeń grupy atomów liganda do środków tych sfer metodą najmniejszych kwadratów. Po każdorazowym dopasowaniu obliczana jest energia oddziaływania ligand - receptor. Za najlepsze dopasowanie uważa się takie, które osiągnęło najniższą energię.

Metoda geometric hashing

Dla każdego liganda przeprowadzana jest następująca operacja: analizowane są wszystkie kombinacje niewspółliniowych trójek punktów (współrzędnych atomowych, położeń grup funkcyjnych lub miejsc występowania określonych cech). Każda taka kombinacja wyznacza trójwymiarowy układ odniesienia, w którym pozostałe punkty mają określone współrzędne. Są one zapamiętywane jako obrazy liganda dla każdej trójki punktów tworzących układy odniesienia. Miejsce wiążące receptora jest analizowane w podobny sposób, z tym, że punktami są pozycje, w których powinny znaleźć się komplementarne elementy ligandów. Następnie oceniane jest podobieństwo wygenerowanych obrazów ligandów i miejsc wiążących, tak by znaleźć najlepsze nałożenie układów odniesienia.

Dokowanie giętkich ligandów

  • dokowanie zespołów konformacji - polega na wygenerowaniu zbioru różnych konformacji dla każdego liganda, a następnie dokowaniu każdej z nich jako osobnej sztywnej molekuły
  • fragmentacja ligandów - w pierwszym etapie ligand jest rozkładany względem wiązań pojedynczych na kilka sztywnych fragmentów. W jednym z wariantów metody dokowany jest główny, największy fragment, a pozostałe są dobudowywane do niego tak, aby odtworzyć układ wiązań. W innym wariancie dokowane są wszystkie fragmenty jednocześnie i następnie sprawdzana jest możliwość odtworzenia molekuły
  • algorytmy stochastyczne - oparte na metodzie Monte Carlo. Kolejne losowe położenia oraz konformacje liganda są akceptowane, bądź odrzucane zgodnie z kryterium Metropolisa w oparciu o obliczoną energię oddziaływania z receptorem.
  • algorytmy genetyczne - położenie oraz konformacja liganda zapisywane są w postaci ciągu liczb pełniącego rolę chromosomu. Po wygenerowaniu początkowej puli stanów układu kolejne ich "pokolenie" tworzone jest dzięki wymianie części chromosomów między poszczególnymi stanami, bądź dzięki "mutacjom punktowym" poszczególnych chromosomów. Następnie dla każdego stanu obliczana jest energia oddziaływania ligand-receptor. Im jest ona korzystniejsza tym większy jest udział danego stanu w generowaniu kolejnego pokolenia.

Uwzględnianie giętkości receptora [16]

  • ruchome fragmenty receptora - części receptora mające istotne znaczenie dla oddziaływania z ligandem (np. określone łańcuchy boczne) traktuje się jako ruchome na podobnej zasadzie jak elementy ligandów w metodzie fragmentacji. Podobnie jak ligand biorą one udział w procedurze dokowania do pozostałej, sztywnej części receptora. Zwykle nie rozważa się wszystkich możliwych konformacji wybranych łańcuchów bocznych receptora, lecz ogranicza się do je takich, które znajdują się w bazach rotamerów (bazy te zawierają zbiór najbardziej prawdopodobnych konformacji łańcuchów bocznych aminokwasów).
  • zespoły konformacji receptora - do dokowania wykorzystuje się kilka konformacji receptora (np. pochodzących z różnych struktur krystalograficznych). Na ich podstawie generowane są struktury receptora będące różnymi kombinacjami stanów wyjściowych. Następnie do każdej z nich, traktowanej jako sztywną, dokowane są ligandy.

Funkcje oceniające energię wiązania

Dokowanie dostarcza nie tylko zbioru sensownych konformacji układu ligand-receptor, ale pozwala również uszeregować je pod względem powinowactwa (m.in. to odróżnia je od zwykłego przeszukiwania baz danych).

Ze względu na stosowane przybliżenia oraz fakt, iż w trakcie dokowania nie są badane dynamiczne właściwości układów, ocena energii swobodnej wiązania liganda do receptora na podstawie "pierwszych zasad" jest niemożliwa. Również ocena samej energii potencjalnej wiązania w oparciu o pola siłowe jest niedokładna ze względu na ograniczone próbkowanie przestrzeni konformacyjnej. Co więcej duży wkład do energii swobodnej wiązania mają efekty solwatacyjne, których ocena ilościowa jest wciąż dużym problemem teoretycznym.

Z powyższych przyczyn nie udało się jak dotąd skonstruować jednej funkcji pozwalającej wiarygodnie obliczyć energię wiązania. Do jej oceny wykorzystuje się natomiast cały zbiór tzw. funkcji oceniających (scoring functions)[18]. Część z nich oparta jest na elementach pól siłowych, a część wykorzystuje potencjały empiryczne.

Funkcje oparte na polach siłowych

Energia swobodna (w zasadzie entalpia swobodna) wiązania G_{bind} może być przedstawiona w następujący sposób:

G_{bind}=G_{nb}+G_{conf} -\Delta G_{solv,elec}- \Delta G_{solv,np}

gdzie: G_{nb} to energia oddziaływań niewiążących ligand - receptor (elektrostatycznych, wiązań wodorowych i Van der Waalsa), G_{conf} - energia danej konformacji liganda, \Delta G_{solv,elec} - zmiana energii związana z desolwatacją ładunków liganda liczona np. przy wykorzystaniu modelu Borna, a \Delta G_{solv,np} - zmiana energii związana z desolwatacją grup nienaładowanych, z reguły liczona w oparciu o zmianę powierzchni dostępnej dla wody.

Funkcje empiryczne

stanowią zdecydowaną większość wśród funkcji oceniających. Z reguły poszczególnym wkładom do energii wiązania przypisuje się określoną wartość skalowaną przez odpowiednie funkcje, zależne od poprawności spełnienia kryteriów wystąpienia danych oddziaływań.

Na przykład funkcja wykorzystywana w programie Ludi ma następującą postać:

G_{bind}=\Delta G_{HB} \sum_{H-bond}{f(\Delta R,\Delta \alpha)}+ \Delta G_{i... ...ph} \sum_{hydroph}{\vert A_{hydroph}\vert}+ \Delta G_{rot} N_{rot}+\Delta G_0

gdzie: \Delta G_{HB} - wkład od pojedynczego wiązania wodorowego, f(\Delta R,\Delta \alpha) - funkcja, której wartość zależy od jakości spełnienia geometrycznych kryteriów wiązania, \Delta G_{ion} oraz \Delta G_{hydroph} wkłady opisujące oddziaływanie odpowiednio grup naładowanych i hydrofobowych, A_{hydroph} funkcja proporcjonalna do pola powierzchni oddziałujących grup hydrofobowych, \Delta G_{rot} wkład "entropowy", wynikający z ograniczenia rotacji fragmentów liganda wokół wiązań pojedynczych, będącego skutkiem związania z receptorem, \Delta G_0 - stała addytywna.

Inna grupa funkcji empirycznych jest oparta na znajomości prawdopodobieństwa, że dane dwie oddziałujące grupy znajdą się w zadanej odległości. Prawdopodobieństwo to jest wyznaczane na podstawie analizy dostępnych struktur (np. krystalograficznych) kompleksów ligand-receptor. Wykorzystując zależność:

A_{ij}(r)=-k_BT \ln p_{ij}(r)

gdzie k_B - stała Boltzmanna, T - temperatura, p_{ij}(r) - prawdopodobieństwo znalezienia grup ij w odległości r, można wyznaczyć wartość potencjału średniej siły (PMF - potential of mean force) A_{ij}(r) oddziaływania grup ij. Na takiej zasadzie oparta jest np. funkcja PMF SCORE wykorzystywana przez wiele programów dokujących.

Consensus scoring

Jak łatwo się domyślić po dużej liczbie istniejących funkcji oceniających żadna z nich nie zapewnia prawidłowej oceny energii swobodnej wiązania. W takiej sytuacji "rzeczywiste" powinowactwo liganda do receptora szacuje się biorąc pod uwagę wyniki kilku lub nawet kilkudziesięciu funkcji oceniających na raz.

W najprostszym przypadku wykorzystuje się rodzaj średniej ważonej po poszczególnych wynikach:

F=\sum_{i} \lambda_i f_i

gdzie F - estymator energii swobodnej, f_i oraz \lambda_i wartości poszczególnych funkcji oceniających oraz ich wagi.

Wyniki dokowania

Dokowanie dużej bazy ligandów dostarcza ogromnych ilości danych. Dla każdego liganda generowany jest zbiór kilkudziesięciu możliwych położeń w miejscu aktywnym oraz tabela zawierająca wartości funkcji oceniających dla każdego z nich.

CADD-Score.png

Przykładowy wynik dokowania: tabela zawierająca wartości różnych funkcji oceniających dla każdego położenia liganda. Im bardziej ujemna wartość "Total Score" tym lepsza konformacja.

Często niestety zdarza się, że prawidłowe umiejscowienie liganda zostaje znalezione, ale ginie wśród innych, "fałszywie pozytywnych" wyników o lepszej punktacji. Jeśli to możliwe warto w grupie dokowanych ligandów zawrzeć takie, których położenie jest znane z danych krystalograficznych. Pozwala to na kalibrację programu dokującego.

W przypadku dużej liczby dokowanych ligandów trzeba przyjąć kryteria, według których automatycznie wybrane zostaną struktury do dalszej analizy. Można np. zdecydować się na wybranie kilku najwyżej ocenionych wyników dla każdego z wąskiej grupy najlepszych ligandów, bądź też wybrać po jednej konformacji dla większej liczby molekuł. Zadanie to ułatwiają programy do analizy wyników dokowań np. Cscore.

Wybrane programy

Projektowanie nowych ligandów i modyfikacja istniejących

Omawiane do tej pory metody, takie jak przeszukiwanie baz danych bądź dokowanie, opierają się na badaniu istniejących już ligandów w celu znalezienia spośród nich najlepszych kandydatów na "molekuły wiodące" (lead compounds). Te z kolei po serii udoskonaleń mają szansę stać się prototypem leku. Warunkiem niezbędnym do zastosowania tych metod jest posiadanie bazy danych struktur analizowanych cząsteczek, co nie zawsze jest łatwe do spełnienia.

Alternatywnym rozwiązaniem jest generowanie nowych ligandów (de novo drug design) [19] w oparciu o znaną strukturę receptora lub farmakofor. Znajduje ono również zastosowanie w procesie modyfikowania istniejących cząsteczek w celu zwiększenia ich powinowactwa do receptora.

Należy zdawać sobie sprawę, że istnieje szacunkowo 10^{60} - 10^{100} możliwych do zsyntetyzowania molekuł mogących być potencjalnymi lekami. Do tego każda z nich może występować w wielu różnych konformacjach. Systematyczne przeszukanie tej przestrzeni jest niemożliwe. Projektowanie nowych ligandów w dużej mierze opiera się zatem na algorytmach stochastycznych. W ogólności polegają one na losowych zmianach struktury i położenia tworzonego liganda, a następnie ocenie jego dopasowania do zdefiniowanego wcześniej miejsca wiążącego (lub farmakoforu). Odrzucanie zmian niekorzystnych prowadzi do powstawania coraz lepiej dopasowanego liganda.

Analiza miejsca wiążącego

Analiza miejsca wiążącego ma na celu zdefiniowanie więzów przestrzennych regulujących powstawanie liganda. Kluczowe jest znalezienie obszarów, w których mogą wystąpić korzystne oddziaływania receptora z ligandem. Z reguły brane są pod uwagę oddziaływania elektrostatyczne, hydrofobowe oraz miejsca akceptorowe i donorowe dla wiązań wodorowych. Wiązania wodorowe dzięki swoim właściwościom kierunkowym w szczególnie dużym stopniu warunkują specyficzność danego miejsca wiążącego. Wynikiem analizy miejsca wiążącego jest mapa przestrzenna obszarów, w których korzystne będzie umieszczenie określonych elementów liganda.

Część programów tworzy taką mapę bazując na rozpoznawaniu zdefiniowanych wcześniej fragmentów molekularnych. Na przykład trafiając na grupę CO w którymś z aminokwasów tworzących miejsce wiążące wykorzystują one zapisane parametry do określenia pozycji, w której powinien pojawić się atom-donor dla wiązania wodorowego. Następnie dopuszczając pewne odchylenia od idealnej geometrii wiązania definiują położenie odpowiedniego obszaru względem rozważanej grupy funkcyjnej.

CADD-Ludi-1.png

Parametry \alpha, \omega i R opisują geometrię wiązania wodorowego dla grupy CO. Na ich podstawie generowany jest zbiór wektorów, wyznaczających miejsca, w których powinna znaleźć się grupa NH, aby stworzyć dobre wiązanie wodorowe.

W innych programach zastosowana jest metoda oparta na siatce przestrzennej. Węzły siatki wyznaczają zbiór punktów wypełniających miejsce wiążące. W punktach tych umieszczane są cząstki próbne oddziałujące w określony sposób z receptorem (np. za pośrednictwem oddziaływań elektrostatycznych). Każdemu punktowi przypisywana jest następnie wartość energii oddziaływania znajdującej się w nim cząstki próbnej dzięki czemu uzyskiwana jest przestrzenna mapa oddziaływań.

CADD-Ludi-2.png

Przykładowy opis obszarów oddziaływań wygenerowany przez program LUDI. Czerwono-białe wektory - miejsca korzystne dla akceptorów wiązań wodorowych, granatowo-niebieskie wektory - miejsca korzystne dla donorów wiązań wodorowych, białe kule - obszary korzystne dla grup hydrofobowych.

W przypadku, gdy struktura miejsca wiążącego jest nieznana, jako postawa do konstrukcji mapy oddziaływań może posłużyć farmakofor wygenerowany w oparciu o zbiór znanych ligandów.

Konstruowanie ligandów

Nowe cząsteczki konstruowane są ze zbioru podstawowych cegiełek. Cegiełkami takimi mogą być pojedyncze atomy, częściej jednak wykorzystywane są większe fragmenty molekularne (np. grupy funkcyjne, pierścienie aromatyczne, kilkuatomowe "łączniki" itp.). Większość programów posiada stały zbiór takich fragmentów, część natomiast umożliwia wygenerowanie go lub uzupełnienie w oparciu o analizę molekularnych baz danych - występujące w nich cząsteczki są "cięte" na kawałki, które po odpowiednim opisaniu (wyliczeniu wartości różnych deskryptorów, zaznaczeniu miejsc, w których mogą łączyć się z innymi fragmentami) trafiają do bazy fragmentów.

Poszczególne fragmenty są wybierane do budowy nowej molekuły w oparciu o ich komplementarność do odpowiednich obszarów miejsca wiążącego. Ich właściwe umiejscowienie odbywa się dzięki wykorzystaniu różnorodnych funkcji oceniających, bardzo podobnych do tych, których używa się w programach służących do dokowania.

Ligand jest budowany i optymalizowany stopniowo. Z reguły wykorzystywane są metody stochastyczne, takie Monte Carlo lub algorytmy genetyczne, ewentualnie połączone z okresową minimalizacją energii potencjalnej. Poszczególne kroki polegać mogą na przesunięciu lub obróceniu całej cząsteczki bądź jej części, dołączeniu, usunięciu lub zamianie wybranych fragmentów. Najczęściej stosowane są następujące strategie:

  • łączenie grup funkcyjnych: obszary miejsca wiążącego istotne dla oddziaływania z ligandem wypełniane są komplementarnymi fragmentami, które w kolejnych krokach łączone są w jedną molekułę
  • rozbudowa zarodka: w jednym z kluczowych obszarów miejsca wiążącego umieszczany jest pojedynczy fragment, do którego dołączane są kolejne tak, aby zapewnić maksymalną liczbę korzystnych oddziaływań
  • chaotyczne mutacje: po umieszczeniu elementu początkowego (fragmentu lub całej cząsteczki) w miejscu aktywnym dopuszczane są wszelkie możliwe kroki: translacje, rotacje, mutacje, rozbudowa lub usunięcie wybranego fragmentu. Kroki te wykonywane w ramach algorytmu genetycznego lub metody Monte Carlo zadaną liczbę razy
  • łączenie węzłów siatki: miejsce wiążące wypełniane jest regularną siatką, w której odległości między węzłami odpowiadają średniej długości wiązania kowalencyjnego. Punkty znajdujące się w obszarach najważniejszych oddziaływań są łączone najkrótszymi możliwymi drogami przez pozostałe węzły siatki. Te węzły, które wchodzą w skład sieci połączeń są zachowywane jako atomy, a wygenerowana tak, z początku niefizyczna "molekuła", jest następnie optymalizowana np. metodą Monte Carlo. Optymalizacja obejmuje zmiany typów atomów, ich translacje, zamianę atomu na grupę funkcyjną itp.

Wiele programów pozwala na wybór jednej z dwóch globalnych strategii budowania liganda: nakierowanej na otrzymanie maksymalnie zróżnicowanych lub maksymalnie dopasowanych cząsteczek. W ramach pierwszej z nich wykonywane są przede wszystkim kroki polegające na mutacjach strukturalnych budowanej cząsteczki, w drugiej nacisk kładziony jest na optymalizację konformacji i położenia w miejscu wiążącym.

Ocena ligandów zaprojektowanych de novo

Powinowactwo liganda do miejsca wiążącego jest na bieżąco szacowane w trakcie jego budowy za pomocą funkcji oceniających podobnych lub identycznych do tych, które stosowane są w dokowaniu. Niestety zdecydowana większość wygenerowanych ligandów mimo bardzo dobrego dopasowania do miejsca wiążącego nie ma szansy stać się działającymi lekami. O możliwości ich praktycznego wykorzystania decyduje bowiem również wiele innych czynników, a przede wszystkim: łatwość syntezy chemicznej oraz grupa właściwości biologicznych określana mianem ADME T (absorption, distribution, metabolism, excretion (ADME) and toxicity (T))[20,21].

Łatwość syntezy chemicznej

Respektowanie walencyjności poszczególnych atomów i wygenerowanie cząsteczki o poprawnej liczbie wiązań chemicznych nie gwarantuje, że będzie ona stabilna lub możliwa do zsyntetyzowania. Niewłaściwe molekuły można wyeliminować sprawdzając, czy nie zawierają one elementów znajdujących się w bazie niepożądanych fragmentów strukturalnych.

Metoda zwiększająca prawdopodobieństwo uzyskania łatwej do zsyntetyzowania struktury polega na łączeniu elementów budulcowych w oparciu o zestaw znanych reakcji chemicznych. Niektóre programy (TOPAS, SYNOPSIS) dysponują bazą fragmentów molekularnych, zawierającą informacje o preferowanych połączeniach, uzyskaną dzięki dekompozycji istniejących cząsteczek polegającej na odwróceniu szeregu reakcji, które służyły do ich syntezy.

Inne programy (CAESA, SEEDS) pozwalają na ocenę już kompletnych ligandów. W oparciu o analizę połączeń międzyatomowych znajdowane są ewentualne drogi syntezy jak również identyfikowane są fragmenty mogące ją utrudniać.

ADME T

Ocena parametrów farmakodynamicznych (opisujących los leku w organizmie) oraz toksyczności jest sprawą trudną, a istniejące programy pozwalają na uzyskanie jedynie zgrubnych wyników.

Wypadkowy efekt działania leków podawanych doustnie zależy od czynników takich jak: absorpcja w przewodzie pokarmowym, stopień dezaktywacji w wyniku przejścia wchłoniętego z jelit leku przez wątrobę (tzw. efekt pierwszego przejścia), zdolności leku do osiągnięcia docelowego miejsca działania (zależnej np. w przypadku leków działających na ośrodkowy układ nerwowy od zdolności przechodzenia przez barierę krew-mózg), a wreszcie od sposobu i wydajności jego eliminacji z organizmu. Część z tych procesów może dodatkowo zależeć od specyficznych interakcji danego leku z różnymi białkami np. białkami transportowymi lub enzymami dezaktywującymi.

Obecnie próbę oceny właściwości farmakodynamicznych leku in silico podejmuje się na dwa sposoby:

  • modelując procesy zachodzące w organizmie
  • budując modele QSAR w oparciu o analizę statystyczną istniejących danych

Podejście oparte na modelowaniu, wykorzystuje do oceny biodostępności symulacje procesów wchłaniania i dystrybucji leku do krwi i tkanek (np. program GastroPlus), w których jako danych wejściowych używa się cech fizykochemicznych leku (np. lipofilności, rozpuszczalności w wodzie, pK grup miareczkowalnych, postaci farmaceutycznej). Istnieją również próby oceny metabolizmu i wydajności dezaktywacji (np. przez cytochrom P450) bazujące na analizie powinowactwa leku do kluczowych białek, z wykorzystaniem ich struktury krystalograficznej.

Metody oparte na modelach QSAR mają długą tradycję. Bazują na takich samych deskryptorach jak modele opisujące aktywność leku, z tym, że parametry, których zmienność pragnie się wyjaśnić są związane z jego biodostępnością (np. stężenie leku w płynie mózgowo-rdzeniowym).

Stosowanym powszechnie kryterium, pozwalającym z dużym prawdopodobieństwem odrzucić nieefektywne molekuły jest zbiór reguł Lipińskiego (Lipinski rules). Oparty jest on na analizie statystycznej właściwości istniejących leków i głosi, że cząsteczka chemiczna ma niewielkie szanse aby być lekiem, jeśli spełnia co najmniej dwa spośród poniższych kryteriów:

  • ma masę > 500Da
  • ma CLOGP (obliczoną wartość logP) > 5
  • zawiera więcej niż 5 miejsc donorowych dla wiązań wodorowych
  • zawiera więcej niż 10 miejsc akceptorowych dla wiązań wodorowych

Modele QSAR wykorzystywane są w podobny sposób do oceny toksyczności leków. Inne podejście oparte jest na systemach eksperckich. Bazują one na analizie struktury potencjalnego leku w poszukiwaniu fragmentów o ustalonych i opisanych wcześniej (zgromadzonych w bazach danych) właściwościach.

Ponieważ rozwijające się metody projektowania leków de novo dostarczają coraz większej liczby potencjalnie sensownych struktur możliwość teoretycznej oceny ich właściwości biologicznych staje się coraz bardziej istotna. Stąd też w ostatnich latach kładzie się duży nacisk na rozwój tego rodzaju oprogramowania.

Niestety, ze względu na złożoność problemu oraz niewielką ilość dostępnych danych nadających sie do budowania szczegółowych modeli QSAR zdolność predykcyjna obecnych narzędzi w zakresie właściwości ADMET wciąż wymaga udoskonalenia.

Wybrane programy

Budowanie liganda:

Ocena ADMET:




Autor: Piotr Setny (piosto @ icm.edu.pl)

Spis literatury

  1. D. Schneidman-Duhovny, R. Nussinov, H. Wolfson, Current Medicinal Chemistry in press.
  2. G. M. A. Milne, Pharmacophore and Drug Discovery. In Encyclopedia of Computational Chemistry; John Wiley & Sons: Chichester, UK:, (1998)
  3. S. Handschuh, M. Wagener, J. Gasteiger, Superposition of Three-Dimensional Chemical Structures Allowing for Conformational Flexibility by a Hybrid Method, J. Chem. Inf. Comput. Sci., 38:220-232 (1998)
  4. J. Holliday, P. Willet, Using a genetic algorithm to identify common structural features in sets of ligands ,Journal of Molecular Graphics and Modelling, 15:221-232 (1997)
  5. P. W. Finn, L. E. Kavraki, J. C. Latombe, R. Motwani, C. Shelton, S. Venkata-subramanian, A. Yao, RAPID: Randomized Pharmacophore Identification for Drug Design, J. Computational Geometry: Theory and Applications, 10:263-272 (1998)
  6. D. Dolata, A. Parrill, W. Walters, CLEW: the Generation of Pharmacophore Hypotheses through Machine Learning, SAR and QSAR in Environmental Research, 9:53-81 (1998)
  7. X. Chen, A. Rusinko III, A. Tropsha, S. Young, Automated pharmacophore identification for. large chemical data sets, J. Chem. Inf. Comput. Sci., 39:887-896 (1999)
  8. http:www.netsci.orgcienceompchem - zbiór artykułów internetowych
  9. M. Hendlich, F. Rippman, G. Barnickel, LIGSITE: automatic and efficient detection of potential small molecule-binding sites in proteins, J. Mol. Graphics, 15:359­363 (1997)
  10. R. C. Willis, Surveying the binding site, Modern Drug Discovery, 28-34, September 2002
  11. ftp:ftp.embl-heidelberg.de/pub/databases
  12. L. Holm, C. Sander, Dali/FSSP classification of three-dimensional protein folds; Nucleic Acids Res., 25(1):231-234 (1997)
  13. H. Xu, Retrospect and prospect of virtual screening in drug discovery, Curr. Op. Med Chem., 2(12):1305-20 (2002)
  14. I. Muegge, M. Rarey, Small Molecule Docking and Scoring, Rev. Comput. Chem., 17:1-60 (2001)
  15. D. B. Kitchen, H. Decornez, J. R. Furr, J. Bajorath, Docking and scoring in virtual screening for drug discovery: methods and applications, Nat Rev Drug Discov, 3(11):935-949 (2004)
  16. I. L. Alberts, N.P. Todorov, P. M. Dean, Receptor flexibility in de novo ligand design and docking, J Med Chem, 48(21):6585-6596 (2005)
  17. J.Y. Trosset, H. A. Scheraga, Prodock: Software package for protein modeling and docking, J. Chem. Phys., 20(4):412-427 (1999)
  18. R. Wang, Y. Lu, S. Wang, Comparative evaluation of 11 scoring functions for molecular docking; J Med Chem. 46(12):2287-303 (2003)
  19. G. Schneider, U. Fechner, Computer-based de novo design of drug-like molecules, Nat Rev Drug Discov 4:649-63 (2005)
  20. I. O. Tudor, Virtual Screening in Lead Discovery: A Viewpoint, Molecules, 7:51­62 (2002)
  21. H. van de Waterbeemd, E. Gifford, ADMET in silico modelling: towards prediction paradise?, Nat Rev Drug Discov., 2(3):192-204 (2003)
  22. H. M. Vinkers et al., SYNOPSIS: SYNthesize and OPtimize System in Silico, J Med Chem., 46(13):2765-73 (2003)
  23. G. Schneider et al, De novo design of molecular architectures by evolutionary assembly of drug-derived building blocks, J Comput Aided Mol Des., 14(5):487-94 (2000)
  24. T. Honma et al, Structure-based generation of a new class of potent Cdk4 inhibitors: new de novo design strategy and library design, J Med Chem., 44(26):4615-27 (2001)