Wpływ protonowania histydyny i stanów rotamerycznych na wirtualne przesiewanie M. tuberculosis RmlC

W celu oceny wpływu protonowania histydyny i stanów rotamerycznych na predykcyjną wydajność receptorów, przeprowadziliśmy wirtualne przesiewanie (VS) dla enzymu Mtb RmlC w oparciu o wyniki poprzedniego badania high-throughput screening (HTS). Poniżej, najpierw przeanalizujemy typowe oddziaływania współkrystalicznego liganda TRH, aby zbadać zależność położenia liganda od protonowania histydyny. W dalszej części przedstawimy kontekst analizy wydajności wzbogacania i mocy predykcyjnej różnych modeli receptorów, omawiając interakcje z receptorem, aby pokazać wpływ różnych stanów protonacji histydyny na VS. Wreszcie, porównujemy przewidywane wartości pKa obliczone przez kilka popularnych pakietów obliczeniowych pKa do stanów protonacji receptora o najlepszej mocy predykcyjnej.

Docking of TRH

Docking współkrystalicznego liganda TRH z powrotem do 36 modeli receptorów został przeprowadzony w celu pokazania zależności pozowania, lub orientacji liganda względem receptora, od protonacji histydyny i stanów rotamerycznych. Intuicyjne chemicznie wzory wiązań wodorowych dla współrzędnych krystalicznych His62 i His119, przedstawione na Rys. 2b, sugerują potencjalne znaczenie wiązań wodorowych w dokowaniu TRH. Dokowanie tego liganda pozwoliło na wstępne zbadanie zależności pozowania od możliwych sieci wiązań wodorowych z receptorem.

Zmienny stan protonowania histydyny ma wyraźny wpływ na przewidywanie pozowania dla wyznaczonego liganda współkrystalicznego. RMSD pozycji dokowania samodokowanego TRH do współrzędnych krystalicznych dla różnych stanów protonacyjnych i rotamerycznych His62 i His119 wahał się od 2,91 do 5,44 Å. Stan protonacyjny obu histydyn z najlepszym średnim RMSD to HIE, który zgadza się z najbardziej prawdopodobnymi stanami protonacyjnymi współrzędnych krystalicznych TRH. Ponadto, we wszystkich przypadkach algorytm dokowania poprawnie przewiduje położenie pirofosforanu liganda, ale duże odchylenia od współrzędnych krystalicznych wynikają głównie z odwrócenia cząsteczek tymidyny i ramnozy wokół pirofosforanu, co skutkuje różnymi wzorami wiązania wodorowego pomiędzy TRH i dwiema histydynami. Wskazuje to na znaczenie sieci wiązań wodorowych z His62 i His119 w przewidywaniu pozy współkrystalicznego liganda TRH. Dlatego też, po zbadaniu zależności pozy od wiązań wodorowych dostarczanych przez dwie histydyny, rozszerzyliśmy nasze badania, aby spojrzeć systematycznie na ranking związków w VS i jak wpływają na to stany protonacyjne i rotomeryczne histydyny.

Virtual screening

Dokowanie molekularne zostało przeprowadzone w celu zbadania zależności rankingu związków od stanów protonacyjnych i rotomerycznych histydyny. Zestaw ligandów zawierał dziesięć aktywnych i 2000 nieaktywnych wybranych losowo z HTS. Zauważyliśmy, że wyniki Tanimoto wskazują, że większość naszych wabików ma niskie podobieństwo do związków aktywnych. Taki zestaw wabików stanowi mniejsze wyzwanie dla algorytmu dokowania, a predykcyjna wydajność VS może ulec pogorszeniu, gdy używane są wabiki o większym podobieństwie do substancji czynnych. Jednakże, celem tego badania nie było sprawdzenie predykcyjnej wydajności algorytmu dokowania per se, ale jak stany protonacji histydyny wpływają na względną wydajność w VS.

Dokowane aktywne ligandy i analog produktu zostały zbadane wstępnie w celu scharakteryzowania ważnych interakcji w miejscu wiązania RmlC. We wszystkich modelach receptora, hydrofobowe interakcje pi-pi przyczyniają się znacząco do wyniku dokowania aktywnych związków w miejscu aktywnym RmlC. Pierwotny związek z HTS, SID7975595, zajmuje wysoką pozycję w większości modeli receptorowych, od 8 do 51 miejsca w 26 z 36 receptorów. Chociaż istnieje tylko ograniczone podobieństwo strukturalne pomiędzy SID7975595 i współkrystalicznym ligandem TRH, pierścień trójcykliczny SID7975595 łatwo zastępuje cząsteczkę tymidynową TRH, podczas gdy pierścień benzimidazolonowy zastępuje cząsteczkę ramnozy, zapewniając strukturalną podstawę inhibicji. Jak pokazano na rycinie 3, w hydrofobowej interakcji pomiędzy substancjami aktywnymi a receptorem często biorą udział Tyr132 i Tyr138 z łańcucha A oraz Phe26 z łańcucha B (należy zauważyć, że część łańcucha B wchodzi do miejsca aktywnego łańcucha A). Poprzez oddziaływanie z istotnymi resztami miejsca wiązania i uniemożliwienie cząsteczkom wody dostępu do Phe26 i Tyr132, substancje czynne zapewniają liczne kontakty hydrofobowe w celu osiągnięcia wysokiego powinowactwa wiązania. Jak omówiono w Sivendran i wsp. , podstawienie grupy etylowej przyłączonej do azotu na pierścieniu trójcyklicznym SID7975595 grupą allilową (np. związek aktywny 77074) jeszcze bardziej zwiększa powinowactwo wiązania poprzez utworzenie jeszcze ściślejszego uszczelnienia hydrofobowego. Dla porównania, zastąpienie tej grupy mniejszą grupą metylową lub atomem wodoru skutkuje niższym powinowactwem wiązania. Oprócz opisanych powyżej kontaktów hydrofobowych, niektóre z substancji aktywnych tworzą również wiązania wodorowe z Ser51, Arg59 i Arg170. Rysunek opisujący interakcje zadokowanych substancji czynnych można znaleźć w Zasobach Online 3.

Co ciekawe, substancje aktywne generalnie nie osiągają polarnych interakcji z His62 i His119. Jak pokazano na Rys. 3, tlen karbonylowy i dwa nitrogeny benzimidazolonu w SID7975595 są zwrócone w kierunku przeciwnym do His62 i His119. Kierunek aromatycznych hydrogenów substancji czynnych często nie pozwala na udział w sieci wiązań wodorowych z dwiema histydynami. Niemniej jednak, różne stany protonacyjne i rotameryczne tych histydyn wpływają na wyniki VS poprzez ich interakcje z wabikami.

Ocena różnic w rankingu

Nierzadko zdarza się, że tylko najlepszy 1% przesiewanych związków może być testowany eksperymentalnie w badaniu VS, ze względu na ograniczone zasoby. Dlatego metryka współczynnika wzbogacenia (EF)1%, która odzwierciedla wydajność wzbogacania bazy danych w górnym 1% (20 zadokowanych związków) biblioteki, staje się szczególnie istotna w ocenie mocy predykcyjnej VS. EF1% waha się od 0 do 80 dla 36 modeli receptorów (Tabela 1), wskazując, że wyniki VS są wrażliwe na protonację i stany rotameryczne His62 i His119 RmlC. Niemniej jednak, 28 z 36 receptorów klasyfikuje więcej niż osiem substancji czynnych w górnych 10% w VS, co odzwierciedla EF10% (Tabela 1), co sugeruje, że większość receptorów jest w stanie rozróżnić substancje czynne i wabiki, gdy bierze się pod uwagę większą część (10%) bazy danych. Wyniki EF sugerują również, że modele receptorów z HIP62 lub HIP119 mają tendencję do słabej wydajności wzbogacania, prawdopodobnie z powodu rozległych sieci wiązań wodorowych z przynętami, jak omówiono później.

Obszar pod krzywą charakterystyki operacyjnej odbiornika (AUC) dla każdego modelu receptora został oceniony w celu przedstawienia wydajności wzbogacania modeli przy różnych stanach protonacji i rotameryzacji His62 i His119. Jak pokazano na Rys. 4a i w Tabeli 1, wartości AUC wszystkich modeli receptorów wahają się od 0,868 do 0,996, wskazując na ogólnie dobrą wydajność predykcyjną (AUC równe 0,5 odpowiada brakowi rozróżnienia między substancjami aktywnymi i wabikami). Ogólnie rzecz biorąc, wynik AUC jest komplementarny do oceny EF dla wydajności predykcyjnej receptora. Podsumowując Tabelę 1, Rys. 4c pokazuje, jak zakres wydajności receptora zależy od dwóch stanów protonacji histydyny i rotameryzacji. Biorąc pod uwagę zakres 25-75% AUC (Rys. 4c, zaznaczony grubszymi liniami), modele His62 wykazują większe zróżnicowanie w zależności od stanu His119. Z drugiej strony, modele His119 mają bardziej spójną wydajność niezależnie od stanów protonacji His62, z wyjątkiem stanu HIP. Wskazuje to, że różne stany protonacji His62 mają mniejszy wpływ niż stany His119 na wydajność receptora w naszym badaniu przesiewowym.

Silniejszą zależność wzbogacenia od stanów protonacji His119 obserwuje się w modelach HIE62 i HIP62. W przypadku stanu HIE62, modele z odwróconym HIP119 (model 6) i odwróconym HIE119 (model 2) dają najwyższą wydajność receptora. Modele 3 i 5 z HID119 i HIP119, odpowiednio, prowadzą do najgorszego wzbogacenia. Badając, dlaczego stan HIE62 ma największą zmienność w AUC, można stwierdzić, że His62 ma albo stosy pi-pi, albo nie ma żadnych oddziaływań z ligandami i tworzy tylko kilka wiązań wodorowych z wysoko postawionymi przynętami. Zatem wydajność receptora zależy od oddziaływań His119 z decoyami. Widać to również badając szeroki zakres wydajności AUC modeli HIP62. Sieci wiązań wodorowych z przynętami zostaną omówione w dalszej części rozdziału.

W celu oceny statystycznej istotności różnicy wartości AUC pomiędzy parą modeli receptorów, wykonaliśmy dwustronny test p na poziomie 95% dla hipotezy zerowej, że para ma statystycznie porównywalne wartości AUC, wobec hipotezy alternatywnej, że ich różnica w wartościach AUC i mocy predykcyjnej jest statystycznie istotna. Wynik przedstawiono w Tabeli 1 Zasobów Online, przy czym podkreślono wartości p mniejsze niż 0,05. Receptory mają średnio ponad 16 wartości p mniejszych niż 0,05, co świadczy o czułości VS na protonację histydyny i stany rotameryczne. Jak można się spodziewać, receptory o najbardziej znaczących różnicach odpowiadają modelom o najwyższych (model 6) lub najniższych wartościach AUC (modele 3, 29 i 5). Model 6 statystycznie lepiej szereguje substancje aktywne w stosunku do wabików niż 26 innych receptorów w zespole. Modele 3, 29 i 5 są wyraźnie gorsze w rankingu substancji czynnych niż odpowiednio 29, 25 i 31 innych receptorów.

Analizę ilościową oddziaływań wiązania wodorowego przeprowadzono dla górnego 1% (20 zadokowanych związków) każdego wyniku VS w celu uwzględnienia licznych oddziaływań wiązania wodorowego z resztami miejsca wiążącego, często obserwowanych w przypadku przynęt. Wyniki wskazują na odwrotną korelację pomiędzy udziałem wiązania wodorowego a wydajnością receptora. Rysunek 4b przedstawia średni udział procentowy wiązań wodorowych w każdym modelu receptora dla 1% zadokowanych związków. Procent wiązania wodorowego jest definiowany jako udział terminu wiązania wodorowego Glide XP w całkowitym wyniku dokowania. Porównanie Rys. 4a, b ujawnia odwrotną zależność pomiędzy procentem wiązania wodorowego a AUC z R2 równym 0,42 (y = -56,18x + 67,95, korelacja jest wykreślona w Online Resource 4). Odwrotna zależność jest powszechnie obserwowana w przypadku modeli z HIP119, odwróconym HIP119 lub HID62, gdzie wysoki procent wiązań wodorowych powodował słabe wzbogacenie. Na przykład, model receptora 29 z HIP62 i HIP119, gdzie obie histydyny prezentują donory wiązań wodorowych skierowane w stronę miejsca aktywnego, ma jedno z najgorszych AUC z powodu wysokiego odsetka wiązań wodorowych w najlepszych trafieniach.

Najwyraźniej, potencjał wiązania wodorowego His119 często określa wydajność receptora. Na przykład, model z HID62 i HIP119 był odchyleniem wśród modeli HID62 na Rys. 4c, z zauważalnie niskim wzbogaceniem w porównaniu z ogólną dobrą wydajnością pozostałych pięciu modeli HID62. Modele HID62 mają wysoką medianę AUC wynoszącą 0,989, pomimo częstego wiązania wodorowego z decoyami z HID62. Jest to spowodowane tym, że stany His119 osiągają niewiele oddziaływań wiązania wodorowego z decoyami. Jedynie w przypadku stanu HIP119 model HID62 tworzy wiązania wodorowe z wieloma decoyami, co skutkuje stosunkowo niskim AUC. Obserwacja ta zgadza się z omówioną powyżej silniejszą zależnością wydajności receptora od stanów protonacji His119. Zasób Online 4 opisuje rozkład AUC i procent wiązania wodorowego wraz z kierunkiem donora lub akceptora wiązania wodorowego z dwóch histydyn skierowanych w stronę receptora.

Powyższe analizy podkreślają utrudniający wpływ wiązania wodorowego do wabików na moc predykcyjną VS, ze względu na różne współrzędne dwóch histydyn o różnych stanach protonacyjnych i rotamerycznych. Rozrzut obserwowanej korelacji z R2 równym 0.42 jest prawdopodobnie przypisany kilku przyczynom, w tym chemicznej naturze zbioru danych wabików, jak również niewielkim różnicom w geometrii każdego receptora podczas minimalizacji w początkowym przygotowaniu białka. Poprzez wyraźne pokazanie wrażliwości wyników wirtualnego screeningu na różne stany protonacyjne i rotameryczne histydyn w miejscu aktywnym, podkreślamy, że należy zachować ostrożność podczas przygotowywania współrzędnych atomowych receptora dla VS, szczególnie biorąc pod uwagę ogólne właściwości badanych ligandów. Dotyczy to zarówno uwzględnienia wiązania wodorowego z ligandem współkrystalicznym i jego wpływu na preparatykę białka, jak i kompleksowej analizy sieci proksymalnych wiązań wodorowych. Do tego momentu porównaliśmy wyniki z różnych pakietów w odniesieniu do naszych wyników VS i omówiliśmy je dalej.

Dokowanie wabików

Różne czynniki prowadzą do różnic w rankingu receptorów, szczególnie w odniesieniu do wabików. Generalnie, wabiki, które uplasowały się wyżej niż substancje aktywne, miały dużą masę cząsteczkową i większy potencjał do tworzenia wiązań wodorowych z receptorem. W tej sekcji, dalej analizujemy częste wzorce interakcji zaobserwowane pomiędzy przynętami i receptorem, ze szczególnym uwzględnieniem modeli receptorów o słabym wzbogaceniu.

Przynęty mają tendencję do posiadania większej masy cząsteczkowej i więcej struktur pierścieniowych niż substancje czynne (Tabela 2). Powoduje to, że dekoksydy zajmują wyższe pozycje w rankingu, ze względu na oddziaływania hydrofobowe przy braku wiązań wodorowych z receptorem. Rysunek 5a przedstawia oddziaływania hydrofobowe uzyskane dzięki dużemu nieaktywnemu związkowi 16952387 w modelu receptora 19. Związek ten często plasuje się w pierwszej piątce w wielu przebiegach VS ze względu na znaczące oddziaływania typu pi-pi stacking z Phe26, Tyr132 i Tyr138. Tendencja ta jest często obserwowana w wirtualnych badaniach przesiewowych, gdzie większe cząsteczki zajmują lepsze pozycje w rankingu w wyniku rozległych interakcji z receptorem.

Wydajność wzbogacania jest szczególnie niska dla receptorów dostarczających obfitych sieci wiązań wodorowych do wabików. Interakcje poprzez His62 i His119 nie były szeroko obserwowane dla substancji czynnych, dlatego związki o większym wkładzie entalpicznym błędnie plasują się korzystniej. Przykład pokazany na Rys. 5b przedstawia interakcję nieaktywnego związku 17388064 w modelu receptora 3 (AUC 0,868), uplasowanego na pierwszym miejscu. W tym receptorze, który jest najgorszy w rankingu związków na podstawie AUC, związek 17388064 tworzy dwa wiązania wodorowe z dwoma histydynami, jedno między jego wodorem hydroksylowym i δ-nitrogenem HIE62, a drugie między jego tlenem hydroksylowym i wodorem na δ-nitrogenie HID119. Ten związek ma pięć donorów wiązań wodorowych i dziewięć akceptorów, duża liczba w porównaniu z odpowiednimi średnimi tych z wabików i substancji czynnych (Tabela 2). Dlatego też, z wysokim udziałem wiązania wodorowego w całkowitym wyniku 34,7 ± 6,62 %, ten związek wabik jest często obserwowany w celu utworzenia co najmniej jednego wiązania wodorowego z którąkolwiek z dwóch histydyn, osiągając w ten sposób wysokie rangi w wielu przebiegach VS.

Dwa inne modele receptora, model 29 z HIP62 i HIP119 oraz model 5 z HIE62 i HIP119, wykazują podobne wzory interakcji z wabikami jak model 3. Te trzy modele mają najniższe wartości AUC, ze średnią 0,870 wśród nich. Jak omówiono powyżej, ich wartości AUC różnią się znacząco od innych receptorów, odzwierciedlając subtelną zależność pomiędzy wiązaniami wodorowymi uzyskanymi poprzez His62 i His119 a słabym wzbogaceniem. Dodatkowa rycina opisująca sieci wiązań wodorowych pomiędzy wabikami i receptorami znajduje się w Online Resource 5.

pKa prediction for His62 and His119

Nasze wyniki wyraźnie pokazują wrażliwość wirtualnego screeningu na protonację histydyny i stany rotameryczne. W wielu obliczeniowych badaniach biofizycznych, stany protonacyjne reszt podatnych na miareczkowanie są określane przy użyciu różnych programów przewidywania pKa. Aby ocenić wydajność tych programów w identyfikacji modelu receptora o najlepszej mocy predykcyjnej w dokowaniu, porównaliśmy wyniki predykcji pKa dla His62 i His119 z PROPKA, Maestro, H++ i MCCE, jak pokazano w Tabeli 3, z obliczonymi wartościami pKa.

Po pierwsze, PROPKA 3.1 przewiduje, że zarówno His62 jak i His119 są neutralne niezależnie od obecności TRH podczas preparatyki. Program ten nie potrafi jednak przypisać stanów rotamerycznych histydyn. Dlatego stan HID, flipped HID, HIE lub flipped HIE musi być określony ręcznie. Podobnie jak PROPKA, program H++, wykorzystujący jedno-strukturalną elektrostatykę continuum, również stwierdza, że obie histydyny są neutralne, chociaż przewidywane wartości pKa różnią się od tych z PROPKA. Program MCCE, który opiera się na wielostrukturalnej elektrostatyce continuum, przewiduje, że His62 jest neutralny, podczas gdy His119 jest protonowany.

Następnie użyliśmy kreatora przygotowania białka w Maestro, aby obliczyć pKa His62 i His119 z i bez TRH. Należy zauważyć, że Maestro jest w stanie zmieniać stany rotameryczne, podczas gdy PROPKA tego nie potrafi. Ostatnia aktualizacja umożliwia Maestro wykorzystanie PROPKA do przewidywania pKa zamiast Epik. W przypadku Epik, Maestro przewiduje zarówno His62 jak i His119 w podwójnie protonowanych stanach, niezależnie od obecności TRH. Co ciekawe, model receptora odpowiadający temu wielohistydynowemu stanowi ma najgorszą moc predykcyjną z AUC równym 0,869. Kiedy używana jest PROPKA, HIP62 i HIE119 są przewidywane dla kompleksu białko-TRH, a HIE62 i HIE119 dla białka apo. Te dwa przewidywania PROPKA w Maestro odpowiadają modelom o umiarkowanej wydajności wzbogacania, z AUC wynoszącym 0,971 dla modelu 25 (HIP62 i HIE119) i 0,942 dla modelu 1 (HIE62 i HIE119), odpowiednio.

Zważywszy, że powyższe przewidywania wykonane przez różne oprogramowanie różnią się znacznie od siebie, należy zachować ostrożność przy stosowaniu tych wyników jako wytycznych do przygotowania białka do wirtualnego przesiewania. Bez dokładnej wiedzy na temat rzeczywistego stanu protonacji receptora, jak również badanych ligandów, trudno jest rozwiązać ten problem. Dlatego sugerujemy, że analiza na małą skalę, jak ta przeprowadzona w tym badaniu, i porównanie z danymi eksperymentalnymi, jeśli są dostępne, może zapewnić dokładniejszy opis stanów protonacyjnych i rotamerycznych mianowanych reszt w receptorach białkowych dla przyszłych badań przesiewowych na większą skalę. Alternatywnie, model, który obejmuje jawne włączenie alternatywnych stanów protonacyjnych i rotamerycznych łańcuchów bocznych podczas dokowania, potencjalnie z informacją przechowywaną w siatce, jak to ma miejsce w przypadku rotowalnych hydroksyli i tioli w programie Glide, może być warty kontynuowania. Badanie wyników w odniesieniu do stanów protonacji i rankingów opartych na interakcjach z histydynami powinno być dokładnie przeanalizowane przed przystąpieniem do badań eksperymentalnych.

Dodatkowo, elastyczność receptora prawdopodobnie wpłynie na stany protonacji jonizowalnych reszt. Chociaż nie było to tutaj bezpośrednio badane, poza minimalizacją każdego receptora po przypisaniu stanów protonacyjnych, elastyczność białka jest wyraźnie ważna dla projektowania i rozwoju leków. Rozważanie przestrzeni konformacyjnej i protonacyjnej w połączeniu staje się szybko niewykonalne przy użyciu metod fizycznych, takich jak te opisane tutaj, ale ulepszone metody próbkowania okazują się obiecujące w radzeniu sobie z takimi trudnościami. Obejmuje to symulacje dynamiki molekularnej o stałym pH, dla których pH jest zewnętrzną zmienną termodynamiczną, używaną do ślepego przewidywania wartości pKa reszt ulegających miareczkowaniu. Efektywne zastosowanie wyników tych symulacji do projektowania molekularnego jest przedmiotem ciągłego zainteresowania. Zespoły równowagowe pochodzące z takich symulacji mogą być stosowane w połączeniu z dokowaniem jako aplikacja złożonego schematu relaksacyjnego, w którym wirtualny screening jest prowadzony z zespołem różnie protonowanych struktur, w celu poprawy wyników wzbogacania. Uwzględnienie elastyczności receptora w przygotowaniu celu doprowadzi do szerszego próbkowania przestrzeni konformacyjnej i protonacyjnej, zwiększając w ten sposób wydajność VS.