Obliczeniowy schemat postępowania w celu priorytetyzacji genów gospodarza powiązanych z metabolitami mikrobiologicznymi w zespole jelita drażliwego z przewagą zaparć

Podsumowanie

Niniejszy protokół integruje przewidywanie celów dla metabolitów mikrobiologicznych, transkryptomikę błony śluzowej odbytnicy, analizę oddziaływań białko-białko i wzbogacenia szlaków, dokowanie molekularne, symulacje dynamiki molekularnej oraz szacowanie wolnej energii wiązania metodą mechaniki molekularnej/powierzchniowej powierzchni Poisson–Boltzmann (MM-PBSA) w celu stworzenia uszeregowionej, służącej do generowania hipotez krótkiej listy kandydackich genów gospodarza związanych z metabolitami oraz strukturalnie priorytetyzowanych kompleksów białko-ligand do dalszych badań eksperymentalnych.

Streszczenie

Nie istnieje ustandaryzowany potok obliczeniowy do systematycznego priorytetyzowania genów gospodarza powiązanych z metabolitami mikrobiologicznymi oraz kompleksów białko-ligand z publicznie dostępnych baz danych chemicznych, genomicznych i strukturalnych. Niniejszy artykuł opisuje ośmioetapowy schemat postępowania, który przyjmuje zdefiniowany przez użytkownika zestaw metabolitów pochodzących z mikrobioty jelitowej i generuje rankingową listę kandydatów na powiązane z metabolitami geny gospodarza, wzbogaconych szlaków biologicznych oraz strukturalnie priorytetyzowanych kompleksów białko-ligand do dalszych badań eksperymentalnych. Potok integruje (i) chemoinformatyczne profilowanie metabolitów; (ii) przewidywanie potencjalnych celów z wielu baz danych przy użyciu narzędzia do przewidywania celów w oparciu o oddziaływania białko-chemikalia i ligandy oraz programu do dokowania molekularnego; (iii) analizę różnicowej ekspresji genów na podstawie publicznie dostępnych danych transkrypcyjnych; (iv) nakładanie celów na geny o różnicowej ekspresji; (v) konstrukcję sieci oddziaływań białko-białko i wzbogacanie szlaków; (vi) dokowanie molekularne z wykorzystaniem programu do dokowania molekularnego; (vii) symulację dynamiki molekularnej przez 200 ns przy użyciu silnika dynamiki molekularnej z polem siłowym białek stosowanym w symulacjach dynamiki molekularnej; oraz (viii) szacowanie wolnej energii wiązania metodą MM-PBSA. Jako przykład praktyczny, dziewięć metabolitów pochodzących z mikrobioty jelitowej lub modyfikowanych przez mikrobiotę, reprezentujących krótkołańcuchowe kwasy tłuszczowe, kwasy żółciowe, metabolity pochodne tryptofanu oraz urolitynę A, zostało przetworzonych z wykorzystaniem publicznego zestawu danych transkrypcyjnych błony śluzowej odbytnicy IBS-C GSE36701. Schemat postępowania wyłonił ranking 17 unikalnych przewidywanych genów powiązanych z metabolitami, które wykazywały różnicową ekspresję w tym zbiorze danych. Dokowanie, symulacja dynamiki molekularnej oraz analizy MM-PBSA pozwoliły na strukturalną priorytetyzację pięciu kompleksów metabolit-białko: kwas litocholowy-VDR, kwas litocholowy-NR1H4/FXR, kwas ursodeoksycholowy-NR1H4/FXR, tryptamina-HTR2A (symulowana w jawnej dwuwarstwie lipidowej 1-Palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC)) oraz urolityna A-CASP3. Protokół został zaprojektowany tak, aby można go było dostosować do innych zestawów metabolitów, zbiorów danych transkrypcyjnych chorób oraz klas celów; wszystkie wyniki są generującymi hipotezy przewidywaniami obliczeniowymi, które wymagają niezależnej replikacji transkrypcyjnej, walidacji na poziomie białka oraz funkcjonalnych testów odpowiedzi na ligandy przed wyciągnięciem wniosków przyczynowych lub terapeutycznych.

Wprowadzenie

Zespół jelita drażliwego z zaparciami (IBS-C) jest powszechnym funkcjonalnym zaburzeniem przewodu pokarmowego charakteryzującym się nawracającym bólem brzucha, zmianą rytmu wypróżnień, wzdęciami i zaparciami, a jego globalna rozpowszechnienie szacuje się na około 10–15% ogólnej populacji^1,2. Obecne terapie farmakologiczne, w tym sekregogi, prokinetyki i leki rozkurczowe, mogą poprawić poszczególne objawy u części pacjentów; jednak odpowiedź na leczenie pozostaje niejednorodna, a trwałe remisje osiąga się rzadko, co odzwierciedla złożoną, wieloczynnikową patobiologię tego schorzenia^1,3,4. W związku z tym niezbędne jest pełniejsze mechanistyczne zrozumienie tego, w jaki sposób sygnały mikrobioty jelitowej są przekazywane na poziomie błony śluzowej, aby sformułować sprawdzalne hipotezy dotyczące nowych celów terapeutycznych.

Mikrobiota jelitowa przyczynia się do homeostazy dolnego odcinka przewodu pokarmowego poprzez produkcję i biotransformację chemicznie zróżnicowanych metabolitów, w tym krótkołańcuchowych kwasów tłuszczowych (SCFA), wtórnych kwasów żółciowych, związków pochodnych tryptofanu oraz metabolitów pochodnych polifenoli, takich jak urolityny^5,6,7,8. Cząsteczki te komunikują się z komórkami gospodarza poprzez szeroki i nie w pełni scharakteryzowany repertuar celów molekularnych, który wykracza daleko poza kanoniczne błonowe receptory wykrywające metabolity, obejmując receptory jądrowe, enzymy cytozolowe, białka modyfikujące histony, prekursory hormonów peptydowych oraz wewnątrzkomórkowe białka sygnałowe⁹. U pacjentów z zespołem jelita drażliwego (IBS) udokumentowano zmiany w składzie społeczności mikroorganizmów jelitowych oraz w profilach metabolitów, co stanowi biologiczne uzasadnienie dla zbadania, czy geny gospodarza związane z reaktywnością na metabolity mikrobiologiczne ulegają zaburzeniom transkrypcyjnym w błonie śluzowej odbytnicy w przebiegu IBS-C¹⁰.

Panel dziewięciu metabolitów zdefiniowano a priori, aby zapewnić zwarty, zróżnicowany chemicznie i biologicznie interpretowalny zestaw małych cząsteczek pochodzących z mikrobioty jelitowej lub modyfikowanych przez nią. Wybór oparto na pięciu kryteriach: reprezentacji głównych klas metabolitów mikrobiologicznych zaangażowanych w sygnalizację gospodarz-mikrobiota; znanej lub prawdopodobnej ekspozycji błony śluzowej dalszych odcinków jelit; dostępności jednoznacznych identyfikatorów PubChem i struktur kanonicznych; wielkości cząsteczki i łatwości analizy strukturalnej dla przewidywania celu na podstawie ligandów oraz dokowania; oraz wcześniejszego prawdopodobieństwa sygnalizacji nabłonkowej, neuroimmunologicznej, enteroendokrynnej, receptorów jądrowych lub związanej z motoryką w przebiegu IBS-C. Wybrany panel obejmował maślan i propionian jako SCFA; kwas chenodeoksycholowy, kwas litocholowy i kwas ursodeoksycholowy jako kwasy żółciowe; tryptaminę, kwas indolo-3-propionowy i kwas indolo-3-mlekowy jako metabolity pochodne tryptofanu oraz urolitynę A jako polifenolowy metabolit pochodzący z mikrobioty jelitowej^5,6,7,8,9,10.

Większość wcześniejszych badań obliczeniowych i eksperymentalnych analizowała poszczególne pary metabolit–receptor lub metabolit–enzym w izolacji, co jest podejściem, które nie oddaje rozproszonego i zbieżnego charakteru sygnalizacji metabolitów mikrobiologicznych w obrębie ścieżek gospodarza^9,11. Integracja wielu etapów analitycznych zapewnia wzajemnie wzmacniającą się moc filtrowania, której żaden pojedynczy etap nie może zapewnić samodzielnie. Obliczeniowa predykcja celu w oparciu o wyselekcjonowane bazy danych pozwala uzyskać szeroki zestaw kandydackich białek gospodarza dla każdego metabolitu. Przecięcie tych danych z transkrypomicznymi danymi istotnymi dla danej choroby znacznie filtruje ten zestaw, pozostawiając jedynie kandydatów, których transkrypty uległy zmianie w kontekście choroby. Następnie analiza wzbogacenia ścieżek oraz analiza sieci oddziaływań białko-białko przypisuje zredukowaną listę kandydatów do znanych modułów biologicznych. Dokowanie molekularne zapewnia wstępną obliczeniową ocenę komplementarności kieszeni wiążącej dla każdego kompleksu kandydackiego, a dodatkowa 200 ns symulacja dynamiki molekularnej (MD) z dekompozycją wolnej energii wiązania MM-PBSA wprowadza wymiar termodynamiczny i rozdzielczość czasową do priorytetyzacji strukturalnej, co nie jest możliwe na podstawie samych wyników dokowania. Wykonywanie każdego kroku niezależnie, bez systematycznej integracji i sekwencyjnego filtrowania, doprowadziłoby do powstania list kandydatów zbyt obszernych, aby były one możliwe do zweryfikowania eksperymentalnie, a także uniemożliwiłoby wykrycie zbieżnej architektury ścieżek.

W ramach niniejszego protokołu przez termin „gen związany z metabolitem” (MAG) rozumiemy gen ludzki, którego produkt białkowy został wskazany jako potencjalny cel molekularny jednego lub więcej metabolitów pochodzących z mikrobioty jelitowej przez co najmniej jedną kuratowaną bazę danych przewidywań obliczeniowych, a którego transkrypty wykazują różnicową ekspresję w zbiorze danych transkrypomicznych związanych z chorobą, wykorzystanym do zademonstrowania schematu postępowania. Ta definicja operacyjna celowo obejmuje receptory błonowe oraz receptory jądrowe, enzymy cytozolowe, białka sygnałowe, prekursory hormonów peptydowych i inne białka wewnątrzkomórkowe. Oznaczenie jako MAG nie stanowi dowodu eksperymentalnego na to, że metabolit wiąże się, tworzy kompleks białko–ligand, aktywuje receptor, zmienia ilość białka lub wywołuje chorobę, lecz jest nominacją opracowaną obliczeniowo w celu sformułowania hipotezy, która wymaga walidacji eksperymentalnej.

Niniejszy protokół opisuje pełny, ośmioetapowy schemat obliczeniowy (Rysunek 1) z wystarczającym poziomem szczegółów operacyjnych, aby umożliwić niezależne powtórzenie badania, adaptację do innych paneli metabolitów lub zbiorów danych o chorobach oraz rozszerzenie na inne konteksty interakcji między gospodarzem a mikrobiotą. Schemat ten jest wyraźnie zdefiniowany jako ramy generowania hipotez i priorytetyzacji strukturalnej, które operują wyłącznie na publicznie dostępnych zasobach omicznych i strukturalnych i nie służą do wnioskowania o zmienionych stężeniach metabolitów, stanach aktywacji receptorów, zmianach w ekspresji białek, aktywności sygnalizacji w dół lub znaczeniu klinicznym wyłącznie na podstawie wyników obliczeniowych. W niniejszej pracy demonstrujemy protokół na przykładzie dziewięciu metabolitów pochodzących z mikrobioty jelitowej lub modyfikowanych przez mikrobiotę oraz publicznego zbioru danych transkrypomicznych błony śluzowej odbytnicy w IBS-C (GSE36701), w celu identyfikacji MAG-ów i priorytetyzacji kompleksów metabolit-białko do późniejszych badań eksperymentalnych.

Protokół

W analizie wykorzystano wyłącznie ogólnodostępne, zanonimizowane dane transkrypcyjne z zestawu GSE36701 oraz ogólnodostępne bazy danych chemicznych, białkowych i strukturalnych. Dostęp do baz danych uzyskano między styczniem a majem 2026 roku. Każda późniejsza data dostępu została odnotowana w osobnej Tabeli Materiałów.

1. Projekt badania, wymagania sprzętowe i programistyczne

1. Zdefiniuj schemat pracy przed rozpoczęciem analizy. Wykorzystaj osiem etapów: wybór metabolitów, przewidywanie celów, analiza różnicowej ekspresji, nakładanie się genów różnicowo eksponowanych (DEG) z celami, analiza interakcji białko-białko (PPI) / wzbogacenia szlaków, dokowanie molekularne, symulacja MD oraz szacowanie MM-PBSA.
2. Odnotuj, że dokowanie, MD i MM-PBSA są wyłącznie narzędziami do priorytetyzacji strukturalnej. Nie interpretuj tych wyników jako dowodów eksperymentalnych na wiązanie, aktywację receptora, zmianę zawartości białka, skuteczność terapeutyczną lub przyczynowość choroby.
3. Potwierdź sprzęt obliczeniowy przed uruchomieniem symulacji MD. Użyj 64-bitowego systemu operacyjnego Linux, procesora 6-rdzeniowego lub lepszego, platformy akceleracji GPU z ≥8 GB VRAM lub równoważnej platformy akceleracji GPU z co najmniej 8 GB VRAM, minimum 32 GB RAM oraz co najmniej 200 GB wolnej przestrzeni dyskowej na każdy system MD.
4. Zapisz podstawowe oprogramowanie: silnik dynamiki molekularnej, program do dokowania molekularnego ligandów z metabolitów do białek docelowych, narzędzie do konwersji formatów plików chemicznych¹², program do generowania trójwymiarowych ligandów, narzędzie do przygotowania ligandów, zestaw narzędzi do przygotowania danych wejściowych do dokowania, program do dokowania molekularnego ligandów z metabolitów do białek docelowych, ogólne środowisko programistyczne, środowisko obliczeń statystycznych z frameworkiem oprogramowania bioinformatycznego oraz pakiet do analizy różnicowej ekspresji genów.
5. Zapisz narzędzia do analizy strukturalnej: internetowe narzędzie do konstruowania systemów błonowych, usługę parametryzacji ligandów zgodną z CHARMM, narzędzie do obliczania energii wiązania w mechanice molekularnej/ciągłym rozpuszczalniku, bibliotekę konwersji topologii i parametrów molekularnych, program do trójwymiarowej wizualizacji molekularnej oraz narzędzie do wizualizacji molekularnej i tworzenia dwuwymiarowych diagramów oddziaływań 2021 (patrz Tabela Materiałów w celu uzyskania linków do pobierania i informacji o wersjach).
6. Zapisz w Tabeli Materiałów/pliku środowiskowym dokładne identyfikatory wersji użytego pola siłowego białek do symulacji dynamiki molekularnej, CGenFF, CHARMM-GUI, frameworka R/bioinformatycznego oraz narzędzia do obliczania energii wiązania w mechanice molekularnej/ciągłym rozpuszczalniku. Brakujące identyfikatory oznacz jako „not recoverable”; nie wnioskuj ich samodzielnie.

2. Wybór metabolitów i charakterystyka chemoinformatyczna

1. Zdefiniuj panel metabolitów przed przystąpieniem do przewidywania celów. Uwzględnij maślan (PubChem CID: 264), propionian (CID: 1032), kwas chenodeoksycholowy (CID: 10133), kwas litocholowy (CID: 9903), kwas ursodeoksycholowy (CID: 31401), tryptaminę (CID: 1150), kwas indole-3-propionowy (CID: 3744), kwas indole-3-mlekowy (CID: 92904) oraz urolitynę A (CID: 5488186).
2. Pobierz kanoniczne ciągi SMILES (Simplified Molecular Input Line Entry System) oraz numery PubChem CID dla każdego metabolitu. Przed przewidywaniem celów zweryfikuj synonimy i zduplikowane struktury. Zapisz końcowe identyfikatory w głównym arkuszu metabolitów.
3. Wprowadź kanoniczne ciągi SMILES do narzędzia internetowego do przewidywania właściwości fizykochemicznych i ADME¹³ (patrz Tabela materiałów). Odnotuj masę cząsteczkową, topologiczną polarną powierzchnię cząsteczki (TPSA), konsensus logP, liczbę donatorów wiązań wodorowych, liczbę akceptorów wiązań wodorowych, liczbę wiązań rotowalnych, przewidywany stopień wchłaniania w przewodzie pokarmowym, przewidywania dotyczące glikoproteiny P oraz alerty Lipinskiego, Vebera, Ghose'a, Egana, Muegge'a i PAINS.
4. Zachowaj metabolity z pomyślnie rozpoznaną strukturą, masą cząsteczkową ≤500 Da i brakiem alertów PAINS. Odnotuj każde niespełnione kryterium oraz decyzję o zachowaniu lub wykluczeniu metabolitu.
5. Przypisz stany jonizacji przed dalszym przewidywaniem celów i dokowaniem. Dla maślanu i propionianu zastosuj zdeprotonowane karboksylany, dla kwasów żółciowych neutralne formy kwasów karboksylowych, dla tryptaminy protonowany amon, a dla pozostałych metabolitów formy neutralne.

3. Przewidywanie potencjalnych ludzkich celów

1. Otwórz narzędzie do przewidywania celów interakcji chemiczno-białkowych¹⁴ (patrz Tabela materiałów). Wprowadź nazwę każdego metabolitu lub numer PubChem CID, wybierz Homo sapiens (ID taksonomiczne: 9606) i ustaw minimalną łączną punktację interakcji na poziomie ≥0,700.
2. Nadaj priorytet kanałom dowodów z eksperymentów i kuratorowanych baz danych w przewidywaniu celów interakcji chemiczno-białkowych. Pobierz pełną tabelę asocjacji białkowych dla każdego metabolitu.
3. Otwórz program do dokowania molekularnego¹⁵ (patrz Tabela materiałów). Prześlij każdy kanoniczny ciąg SMILES z wybranym gatunkiem Homo sapiens i zachowaj cele o prawdopodobieństwie ≥0,70.
4. Połącz wyniki przewidywania celów interakcji chemiczno-białkowych oraz programu do dokowania molekularnego jako zbiór sumy dla każdego metabolitu. Zachowaj każdy cel spełniający próg którejkolwiek z baz danych i usuń dokładnie powtarzające się wpisy symboli genów.
5. Znormalizuj wpisy białek do symboli genów zatwierdzonych przez Komitet Nazewnictwa Genów HUGO (HGNC), korzystając z mapowania identyfikatorów białek na ustandaryzowane symbole genów zatwierdzone przez HGNC lub zintegrowanej bazy informacji o genach ludzkich (patrz Tabela materiałów). Rozwiąż aliasy, przestarzałe symbole i adnotacje izoform, aby przypisać jeden symbol genu do jednego białka.
6. Sklasyfikuj każdy cel jako receptor błonowy, receptor jądrowy, enzym, wewnątrzkomórkowe białko sygnałowe, hormon peptydowy, białko powiązane z hormonami lub inne białko wewnątrzkomórkowe. Zapisz klasę w tabeli celów.

4. Zbiór danych transkryptomicznych i analiza różnicowej ekspresji genów

1. Uzyskaj dostęp do zestawu danych GSE36701 poprzez internetowe narzędzie NCBI do analizy różnicowej ekspresji genów^16,17 (patrz Tabela materiałów). Odnotuj, że zestaw danych zawiera dane o ekspresji z biopsji błony śluzowej odbytnicy grup z IBS-C, zespołem jelita drażliwego z przewagą biegunek (IBS-D), IBS poinfekcyjnym oraz grupą zdrowych ochotników¹⁸.
2. Przeszukaj bazy GEO oraz ArrayExpress w celu znalezienia niezależnej kohorty walidacyjnej. Użyj kombinacji haseł: IBS-C, constipation-predominant irritable bowel syndrome, rectal mucosa, colonic mucosa, biopsy, transcriptome, microarray oraz RNA-seq. Zapisz nazwy repozytoriów, terminy wyszukiwania, datę wyszukiwania oraz informację, czy zidentyfikowano porównywalny zestaw danych walidacyjnych.
3. Uruchom internetowe narzędzie do analizy różnicowej ekspresji genów z rekordu GSE36701 (patrz Tabela materiałów). Przypisz 18 próbek IBS-C do grupy IBS-C, 40 zdrowych ochotników do grupy kontrolnej, a próbki IBS-D oraz IBS poinfekcyjnego pozostaw bez przypisania.
4. Przeprowadź analizę różnicowej ekspresji, korzystając z frameworka pakietu do analizy różnicowej ekspresji genów z korekcją wskaźnika fałszywych odkryć (FDR) metodą Benjamini-Hochberga¹⁹. Pobierz pełną tabelę wyników zawierającą ID sondy, symbol genu, pełną nazwę genu, logFC, AveExpr, umiarkowaną statystykę t, surową wartość P oraz skorygowaną wartość P.
5. Zredukuj sondy do wpisów na poziomie genów. Usuń sondy pozbawione symboli genów; w przypadku powielonych symboli zachowaj sondę z najniższą wartością FDR; w przypadku remisów wybierz większą bezwzględną wartość logFC.

5. Analiza nakładania się celów i genów różnicowo wyrażanych oraz ocena statystyczna

1. Wykonaj część wspólną każdej przewidywanej listy celów specyficznych dla danego metabolitu z listą DEG na poziomie genów przy FDR < 0,05. Zapisz nakładające się geny, metabolit pochodzenia, logFC, skorygowaną wartość P oraz kierunek ekspresji.
2. Połącz listy nakładania specyficzne dla metabolitów w nieredundantną listę MAG. Oblicz całkowitą liczbę przewidywanych celów, liczbę nakładania specyficznych dla metabolitów oraz całkowitą liczbę unikalnych MAG.
3. Ocenia spójność kierunkową na poziomie sond dla genów z wieloma sondami. Oznacz każdy gen, dla którego sondy wykazują rozbieżność w kierunku ekspresji.
4. Sporządź tabelę kontyngencyjną do dokładnego testu Fishera, wykorzystując całkowitą liczbę wpisów po agregacji genów, całkowitą liczbę DEG, całkowitą liczbę unikalnych przewidywanych celów oraz zaobserwowane MAG. Oblicz jednostronną wartość P, iloraz szans oraz 95% przedział ufności przy użyciu implementacji dokładnego testu Fishera.
5. Jeśli tło DEG przekracza 50%, przedstaw nakładanie jako opisowe, a nie jako niezależnie zwalidowane wzbogacenie. Jednolitną regulację w dół traktuj jako opisowy wzorzec kierunkowy, chyba że przeprowadzono oddzielny test kierunkowości.

6. Analiza sieci oddziaływań białko-białko i wzbogacenie szlaków

1. Prześlij pełną, unikalną listę MAG do narzędzia do budowy sieci oddziaływań białko-białko i wzbogacania szlaków²⁰ (patrz Tabela materiałów). Wybierz gatunek Homo sapiens i ustaw minimalną wartość wyniku oddziaływania na 0.700.
2. Wyeksportuj połączoną sieć budowy oddziaływań białko-białko i wzbogacania szlaków oraz pełną tabelę oddziaływań. Jeśli analiza tekstowa (text-mining) generuje nienaturalnie gęstą topologię, odznacz tę opcję i zachowaj kanały eksperymentalne, koekspresyjne oraz bazodanowe.
3. Wygeneruj podsieci klas metabolitów dla MAG powiązanych z SCFA, kwasami żółciowymi oraz tryptaminą/serotoniną. Użyj tych samych ustawień organizmu i poziomu ufności, co w narzędziu do budowy sieci oddziaływań białko-białko i wzbogacania szlaków.
4. Przeprowadź analizę budowy sieci oddziaływań białko-białko i wzbogacania szlaków w oparciu o Kyoto Encyclopedia of Genes and Genomes (KEGG)²¹, Reactome²² oraz Procesy Biologiczne Gene Ontology (GO)^23,24. Zastosuj korektę Benjamini–Hochberg BH FDR <0.05 i wyeksportuj wszystkie tabele wzbogacenia.

7. Dokowanie molekularne

1. Pobierz eksperymentalnie określone struktury receptorów z Protein Data Bank RCSB PDB (Research Collaboratory for Structural Bioinformatics)²⁵ (patrz Tabela materiałów). Dla pięciu priorytetowych kompleksów białko-ligand wykorzystaj VDR/1DB1, NR1H4/FXR/3DCT, CASP3/2DKO oraz HTR2A/6A93.
2. Przygotuj każdy receptor, zachowując łańcuch A i usuwając cząsteczki wody, ligandy kokrystalizowane, kofaktory, jony oraz rekordy HETATM niebędące białkami. W przypadku struktury 6A93 przed przygotowaniem receptora usuń segment fuzji z lysozymem T4.
3. Dodaj wodory polarne, przypisz ładunki Gasteigera i zapisz każdy receptor w formacie PDBQT za pomocą zestawu narzędzi do struktur molekularnych. Przed konwersją do PDBQT sprawdź stany protonacji histydyny w miejscu wiązania i udokumentuj wybrane stany.
4. Wygeneruj trójwymiarową strukturę każdego liganda w narzędziu do konwersji plików struktur chemicznych. Przeprowadź minimalizację energii przy użyciu Universal Force Field (UFF) przez 500 kroków, przypisz stan jonizacji dla pH 7,4, przypisz ładunki Gasteigera i zapisz plik w formacie PDBQT.
5. Zdefiniuj przestrzeń dokowania (docking box) o wymiarach 25 Å x 25 Å x 25 Å, wyśrodkowaną na centroidzie liganda kokrystalizowanego. Użyj centrów o współrzędnych (10, 19, 33) dla VDR, (137, 31, 78) dla FXR, (37, 34, 32) dla CASP3 oraz (12, −1, 61) dla HTR2A.
6. Przeprowadź dokowanie molekularne ligandów będących metabolitami do białek docelowych^26,27, przyjmując parametry: exhaustiveness = 8, seed = 42, num_modes = 9 oraz energy_range = 3 kcal/mol. Zapisz najwyżej oceniony wynik Vina (Vina score) oraz wartości średniokwadratowego odchylenia (RMSD) dla wszystkich póz.
7. Wybierz tryb 1 (mode 1) dla każdego priorytetowego kompleksu. Wygeneruj dwuwymiarowe schematy oddziaływań ligand-pozostałość w narzędziu do wizualizacji molekularnej i dwuwymiarowych diagramów oddziaływań oraz trójwymiarowe widoki receptor-ligand w programie do trójwymiarowej wizualizacji molekularnej.
8. Wykonaj kontrole redokowania dla VDR/1DB1 i FXR/3DCT. Zaakceptuj konfigurację dokowania receptora, gdy RMSD dla ciężkich atomów wynosi <2,0 Å w stosunku do pozy krystalograficznej.
9. Przeprowadź kontrole dokowania krzyżowego (cross-docking), dokując kwas litocholowy (LCA) do CASP3 oraz tryptaminę do VDR. Porównaj wyniki dla ligandów cognate i non-cognate oraz odnotuj przypadki, w których różnica wyników wynosi <1,0 kcal/mol.

8. Symulacja dynamiki molekularnej

1. Wygeneruj parametry ligandu za pomocą usługi parametryzacji ligandów kompatybilnej z CHARMM²⁸ (patrz Tabela Materiałów). Sprawdź wszystkie wyniki kar (penalty scores) i zaznacz każdy parametr z karą >50.
2. Konwertuj pliki strumieniowe ligandów na pliki .itp i .prm kompatybilne z silnikiem dynamiki molekularnej za pomocą skryptu do konwersji topologii pola siłowego. Połącz pliki topologii ligandu i białka dla każdego kompleksu.
3. Zastosuj ponowne rozdzielenie masy wodoru (hydrogen mass repartitioning) przy użyciu biblioteki konwersji topologii molekularnej i parametrów. Wygeneruj topologie wodne, wykorzystując pole siłowe białek stosowane w symulacjach dynamiki molekularnej²⁹ oraz jawny trójpunktowy model wody³⁰.
4. Solwatuj kompleksy wodne w pudełku dodekaedralnym z odstępem od krawędzi substancji rozpuszczonej wynoszącym co najmniej 1.2 nm. Zneutralizuj systemy i dodaj NaCl do stężenia 0.15 M.
5. Zbuduj system błonowy tryptamina-HTR2A za pomocą internetowego narzędzia do konstrukcji systemów błonowych^31,32,33 (patrz Tabela Materiałów). Użyj współrzędnych receptora wyrównanych z bazą danych orientacji białek błonowych³⁴ (patrz Tabela Materiałów), czystej dwuwarstwy POPC, warstw wody o grubości 22.5 Å oraz 0.15 M NaCl.
6. Zminimalizuj energię wszystkich systemów metodą najstromniejszego zejścia przez maksymalnie 50 000 kroków. Przed ekwilibracją potwierdź zbieżność przy Fmax <1000 kJmol^-1nm^-1.
7. Przeprowadź ekwilibrację systemów wodnych w fazach zespołu o stałej liczbie cząsteczek, objętości i temperaturze (NVT) oraz zespołu o stałej liczbie cząsteczek, ciśnieniu i temperaturze (NPT). Ekwilibruj system błonowy, korzystając z sześciostopniowego internetowego przepływu pracy przygotowania i ekwilibracji systemów molekularnych z stopniowo zwalnianymi więzami.
8. Przeprowadź 200 ns produkcyjnej dynamiki molekularnej (MD) dla wszystkich pięciu kompleksów. Zastosuj krok czasowy 4 fs z ponownym rozdzieleniem masy wodoru (HMR), termostat V-rescale w 310 K, barostat Parrinello-Rahman przy 1 bar, elektrostatykę Particle mesh Ewald (PME)³⁵ oraz więzy LINCS³⁶.
9. Przeanalizuj końcowe trajektorie za pomocą narzędzi do analizy trajektorii dynamiki molekularnej. Oblicz RMSD szkieletu, fluktuację średniokwadratową (RMSF) atomów Cα, promień żyroskopowy, powierzchnię dostępną dla rozpuszczalnika (SASA) oraz wiązania wodorowe białko-ligand, wykorzystując ostatnie 150 ns jako główny okno analizy.

9. Szacowanie wolnej energii wiązania metodą MM-PBSA

1. Wyodrębnij migawki trajektorii do analizy MM-PBSA. Wykorzystaj 2001 klatek dla każdego kompleksu w środowisku wodnym oraz 201 przetworzonych klatek dla podsystemu HTR2A osadzonego w błonie.
2. Uruchom narzędzie do obliczeń energii wiązania metodą mechaniki molekularnej i rozpuszczalnika kontinuum³⁷ z solwatacją Poissona-Boltzmanna, stałą dielektryczną wewnętrzną = 1, stałą dielektryczną zewnętrzną = 80, solwatacją niepolarną opartą na SASA oraz bez korekty entropowej. Podaj średnią swobodną energię wiązania oraz odchylenie standardowe.
3. Przeprowadź dekompozycję na poszczególne reszty dla wszystkich pięciu kompleksów. Wskaż reszty stabilizujące i destabilizujące o wartościach bezwzględnych wkładów ≥0,5 kcalmol⁻¹.

    

Wyniki

Kandydackie cele związane z metabolitami

Dziewięć metabolitów wygenerowało heterogeniczne zestawy przewidywanych celów w programach do predykcji oddziaływań chemicznych z białkami oraz dokowania molekularnego. Propionian, tryptamina, kwasy żółciowe i urolityna A wskazały kilka celów o znanej roli w sygnalizacji żołądkowo-jelitowej. Przewidywany krajobraz celów obejmował kanoniczne receptory błonowe, receptory jądrowe, enzymy wewnątrzkomórkowe, białka sygnałowe oraz białka związane z hormonami peptydowymi. Wyniki poniższe są zatem opisane jako geny związane z metabolitami (MAGs), a nie jako odkrycia dotyczące wyłącznie receptorów (Tabela 1).

Weryfikacja w odniesieniu do zgłoszonych oddziaływań metabolit-białko

Aby porównać wyniki przewidywania celów z istniejącą wiedzą eksperymentalną, przewidywane relacje między metabolitami a celami sklasyfikowano w trzech poziomach dowodów: (i) bezpośrednie lub bliskie na poziomie klasy interakcje metabolit-białko potwierdzone eksperymentalnie, w których zgłoszono, że metabolit lub blisko spokrewniony metabolit endogenny wiąże, aktywuje, hamuje lub reguluje funkcjonalnie kodowane białko; (ii) interakcje wspierane przez szlaki metaboliczne lub klasę celów, w których przewidywany cel należy do ustalonego szlaku reagującego na metabolity lub rodziny receptorów, ale bezpośrednie dowody dla konkretnej pary metabolit-białko są ograniczone; oraz (iii) powiązania oparte wyłącznie na obliczeniach, dla których w przeglądanej literaturze nie zidentyfikowano bezpośredniej interakcji eksperymentalnej. To porównanie służyło do osadzenia przewidywanych MAG w kontekście, a nie do ich walidacji.

Kilka przewidywań powtórzyło wcześniej opisane zjawiska biologiczne. Interakcję Propionate-FFAR2 uznano za wspartą eksperymentalnie, ponieważ FFAR2/GPR43 jest kanonicznym receptorem krótkołańcuchowych kwasów tłuszczowych. Interakcję Butyrate-HDAC3 sklasyfikowano jako wspartą eksperymentalnie lub klasowo, ponieważ maślan jest uznanym inhibitorem deacetylazy histonowej, a przewidywane nakładanie się dotyczyło członka rodziny HDAC. Przewidywania związane z kwasami żółciowymi obejmujące NR1H4/FXR i VDR uznano za wsparte uznaną biologią jądrowych receptorów kwasów żółciowych, szczególnie w przypadku hydrofobowych kwasów żółciowych, takich jak LCA; przewidywania dotyczące FXR związane z kwasem ursodeoksycholowym (UDCA) interpretowano z ostrożnością, ponieważ UDCA jest zazwyczaj słabszym lub zależnym od kontekstu ligandem FXR. Przewidywania dotyczące HTR1B, HTR2A, HTR2B i HTR6 związane z tryptaminą sklasyfikowano jako wsparte szlakiem serotonergicznym, a nie jako potwierdzone bezpośrednie oddziaływania specyficzne dla receptora, ponieważ tryptamina jest mikrobiologiczną monoaminą pochodzącą z tryptofanu, a receptory serotoninowe są uznanymi regulatorami motoryki i wydzielnictwa żołądkowo-jelitowego. Interakcję Urolithin A-CASP3 uznano za wspartą szlakiem na podstawie opublikowanych powiązań między urolityną A a odpowiedziami apoptotycznymi/zależnymi od kaspaz, a nie na podstawie bezpośrednich dowodów wiązania CASP3. Interakcje Indole-3-lactic acid-KYAT1 oraz indole-3-propionic acid-KYAT1 pozostawiono jako hipotezy wyłącznie obliczeniowe, ponieważ szersza literatura wspiera sygnalizację gospodarza poprzez mikrobiologiczne pochodne indolu, ale nie bezpośrednie wiązanie KYAT1 przez te konkretne metabolity^7,8,38,39,40.

W związku z tym w Tabeli 1 odróżniono nominację celu obliczeniowego od poziomu wcześniejszego wsparcia eksperymentalnego lub ścieżkowego. Podano w niej również dla każdego celu źródło predykcji (predykcja celu interakcji chemiczno-białkowej, program dokowania molekularnego lub oba), łączny wynik interakcji dla predykcji celu interakcji chemiczno-białkowej oraz prawdopodobieństwo z programu dokowania molekularnego w przypadku, gdy cel został zidentyfikowany za pomocą programu dokowania molekularnego. Przewidziane cele bez bezpośrednich wcześniejszych dowodów eksperymentalnych opisano jako kandydackie geny związane z metabolitami, które wymagają niezależnej walidacji na poziomie białka i odpowiedzi na ligand.

Pokrycie między przewidywanymi celami a genami różnicowo wyrażanymi w IBS-C

Przecięcie połączonych list przewidywanych celów oraz wyników różnicowej ekspresji na poziomie genów pozwoliło zidentyfikować 17 unikalnych przewidywanych genów związanych z metabolitami, których ekspresja była istotnie różnicowa w porównaniu osób z IBS-C i zdrowych ochotników. Wszystkie 17 genów wykazywało obniżoną ekspresję. Zestaw ten obejmował receptory błonowe i jądrowe (CASR, FFAR2, GPR68, HTR1B, HTR2A, HTR2B, HTR6, NR1H4, TBXA2R, VDR) oraz białka niereceptorowe (CASP3, GCG, GNAQ, GPHN, HDAC3, KYAT1, MLN) (Tabela 1, Rysunek 2A,B).

Wszystkie 17 MAGs spełniło kryterium poziomu fałszywych odkryć (FDR) poniżej 0,05; 16 z 17 spełniło bardziej rygorystyczny warunek FDR < 0,001, natomiast pozostały gen (HTR1B) był istotny przy FDR < 0,05. Siedem z 17 celów (CASP3, GCG, GNAQ, GPHN, GPR68, HDAC3, TBXA2R) spełniło oba warunki: FDR < 0,001 oraz wartość bezwzględną log₂ krotności zmiany przekraczającą 1,0 (zakres logFC od −1,34 do −1,10), co wskazuje na silną i spójną downregulację dla tej podgrupy. Pozostałe cele wykazały umiarkowaną, ale istotną statystycznie downregulację (zakres |logFC| od 0,45 do 0,97). Ten jednolity wzorzec opisowy zinterpretowano z ostrożnością, biorąc pod uwagę charakterystykę ekspresji w całym genomie dla danego zbioru danych (patrz. poniższa ocena statystyczna).

Ocena statystyczna nakładania się genów docelowych i DEG

Aby formalnie ocenić istotność statystyczną nakładania się 17 genów, zastosowano jednostronny dokładny test Fishera, przyjmując 17 przewidywanych genów docelowych jako zbiór zapytywania oraz wszystkie 18 296 unikalnych wpisów zredukowanych do poziomu genów, wykrytych w GSE36701, jako tło genomiczne. Z tego tła 17 296 genów (94,5%) wykazywało różnicową ekspresję przy FDR < 0.05, co odzwierciedla niemal powszechną supresję transkrypcyjną w porównaniu błony śluzowej odbytnicy IBS-C. Wszystkie 17 przewidywanych genów docelowych znajdowało się wśród genów o różnicowej ekspresji (obserwowany stopień nakładania się 17/17, 100%). Biorąc pod uwagę 94,5-procentowy wskaźnik różnicowej ekspresji w tle, oczekiwana liczba nakładających się genów dla dowolnego losowo wybranego zbioru 17 genów wynosi 16,1. Dokładny test Fishera dał wynik p = 0.384 przy ilorazie szans z korektą na ciągłość wynoszącym 2.03 (95% przedział ufności 0.12–33.73), co nie było statystycznie istotne przy α = 0.05 (Rycina 3A–C).

Wynik ten wskazuje, że obserwowana pokrywność 17/17 nie przekracza pokrywności spodziewanej przypadkowo przy profilu ekspresji w całym genomie dla tego zbioru danych. W związku z tym ustalenia te należy interpretować jako opisowy wzorzec kierunkowy, w którym wszystkie 17 przewidywanych celów było spójnie i istotnie obniżonych w tkance błony śluzowej odbytnicy w IBS-C, a nie jako dowód na statystyczne wzbogacenie lub niezależną walidację w odniesieniu do tła genomicznego. Formalne testowanie wzbogacenia wymagałoby replikacji w zbiorach danych transkryptomicznych o bardziej selektywnych profilach ekspresji różnicowej, w których istotność osiąga znacznie mniej niż połowa wszystkich genów. Należy podkreślić, że jednolity spadek ekspresji wszystkich 17 nakładających się genów jest obserwacją opisową, a nie oddzielnie zwalidowanym wynikiem statystycznym, ponieważ tło ekspresji różnicowej w tym zbiorze danych jest samo w sobie w przeważającej mierze obniżone; wspólny kierunek spadkowy wśród nakładających się genów był zatem oczekiwany i nie poddano go formalnemu testowi kierunkowości. Ten jednolity kierunek nie powinien zatem być interpretowany jako niezależny statystyczny dowód na skoordynowaną regulację specyficzną dla metabolitów.

Wzorce specyficzne dla metabolitów

Propionian posiadał największą liczbę nakładających się genów, w tym CASR, FFAR2, GCG, GNAQ, GPHN, GPR68, MLN i TBXA2R, co sugeruje możliwy udział sygnalizacji zależnej od krótkołańcuchowych kwasów tłuszczowych oraz związanej z białkiem Gq. Maślan wykazał nakładanie się z HDAC3, co jest zgodne z biologią deacetylazy histonowej związaną z maślanem, choć sama regulacja w dół mRNA nie potwierdza zmiany odpowiedzi na maślan. Nakładania związane z kwasami żółciowymi obejmowały receptory jądrowe VDR i NR1H4, które są uznanymi efektorami sygnalizacji kwasów żółciowych w jelitach^38,39. Tryptamina wykazała nakładanie się z HTR1B, HTR2A, HTR2B i HTR6, co wskazuje na sygnalizację serotonergiczną jako potencjalny moduł — system o dobrze udokumentowanej roli w motoryce i wydzielniczości przewodu pokarmowego⁴⁰. Kwas indole-3-mlekowy i kwas indole-3-propionowy nakładały się z KYAT1, a urolityna A nakładała się z CASP3.

Wzbogacanie szlaków

Analiza wzbogacenia funkcjonalnego 17 wspólnych genów zidentyfikowała szlaki związane z sygnalizacją downstream receptorów sprzężonych z białkiem G (GPCR), sygnalizacją Gαq, wiązaniem ligandów GPCR, synapsami serotoninergicznymi, interakcjami neuroaktywnych ligandów z receptorami, transdukcją sygnałów wapniowych, sygnalizacją cAMP oraz sekrecją hormonów peptydowych. Wyniki te są zgodne ze składem zestawu genów i potwierdzają jego spójność biologiczną, jednak odzwierciedlają one adnotacje funkcjonalne zgłoszonych genów, a nie stanowią niezależnego dowodu na aktywność na poziomie szlaków.

Struktura sieci oddziaływań białko-białko

Konstrukcję sieci interakcji białko-białko oraz analizę wzbogacenia ścieżek zinterpretowano w ramach trzech uzupełniających się sieci. W połączonej meta-sieci obejmującej 17 genów (Sieć 1) najbardziej wyraźną strukturą wspieraną przez adnotacje był komponent sygnalizacyjny GPCR/Gαq zorientowany wokół GNAQ, łączący GNAQ z genami związanymi z receptorami, takimi jak TBXA2R, CASR, HTR2A i HTR2B. Zachowano również ograniczoną łączność receptorów serotoninowych, najsilniej pomiędzy HTR2A i HTR2B, podczas gdy kilka innych genów pozostało izolowanych lub słabo połączonych przy wybranym progu ufności. Sieć specyficzna dla propionianu (Sieć 2) wykazała bardziej ograniczoną topologię, w której GNAQ zachował wspierane adnotacjami połączenia z CASR i TBXA2R, natomiast FFAR2, GPR68, GCG, GPHN i MLN były izolowane lub słabo połączone. Sieć tryptaminy/serotoniny (Sieć 3) obejmowała HTR1B, HTR2A, HTR2B i HTR6; w obrębie tego podzbioru HTR2A i HTR2B wykazały główne połączenie wspierane przez adnotacje, podczas gdy HTR1B i HTR6 nie były bezpośrednio połączone przy wybranym progu (Rysunek 4A–C).

Dokowanie molekularne

Dokowanie molekularne przeprowadzono dla pięciu wybranych kompleksów metabolit-białko. Pary kwas żółciowy-receptor jądrowy wykazały korzystniejsze wyniki Vina niż urolityna A-CASP3 i tryptamina-HTR2A. Najlepszy wynik uzyskał kompleks LCA-VDR (−10.0 kcal/mol), a następnie LCA-NR1H4/FXR (−9.9 kcal/mol) oraz UDCA-NR1H4/FXR (−9.4 kcal/mol). Urolityna A-CASP3 i tryptamina-HTR2A uzyskały niższe, lecz wciąż dopuszczalne wyniki na poziomie −7.1 kcal/mol (Tabela 2).

W przypadku kompleksu LCA-VDR (PDB ID: 1DB1) przewidywana poza była stabilizowana przez konwencjonalne wiązanie wodorowe między tlenem grupy karboksylowej LCA a Ser278 (4,29 Å), wraz z rozległymi kontaktami hydrofobowymi z udziałem Leu230, Val234, Trp286, Val300, His305, Tyr295, Leu233 oraz His397, a także dodatkowymi kontaktami van der Waalsa z Met272, Leu313, Ile271, Ile268, Leu309, Phe422, Val418, Ala231, Ala303, Cys288, Ser275 i Phe150. Poza o najwyższym rankingu uzyskała wynik Vina wynoszący −10,0 kcal/mol, rozmiar kawerny 2055 ų oraz środek siatki (10, 19, 33) (Tabela 3, Rycina 5A,B).

W przypadku kompleksu LCA-NR1H4/FXR (PDB ID: 3DCT), wynik dokowania wynoszący −9.9 kcal/mol był powiązany z przewidywanymi wiązaniami wodorowymi z udziałem His294 i Ile335, oddziaływaniem π-Sigma z His294 oraz hydrofobowymi kontaktami typu Alkyl lub π-Alkyl z udziałem Met290, Met328, Ala291, Leu287, Ile352 i His447, przy czym dodatkowe oddziaływania van der Waalsa wspierały dopasowanie szkieletu steroidowego w kieszeni FXR (Tabela 4, Rycina 6A,B).

Przewidywana poza kompleksu UDCA-NR1H4/FXR (PDB ID: 3DCT) wykazała konwencjonalne wiązanie wodorowe z His447 (3.66 Å), kolejne wiązanie wodorowe z Gly322 (3.46 Å), oddziaływanie π-anionowe z Val325 (4.96 Å) oraz wiązanie węgiel-wodór z Trp469 (4.51 Å). Mapa oddziaływań zidentyfikowała również niekorzystne kontakty donor-donor z Arg395 (3.89 Å) oraz Gln396 (3.40 Å), co sugeruje, że niższa wartość Vina score dla UDCA w porównaniu z LCA w tej samej kieszeni receptora może wynikać z mniej korzystnej lokalnej geometrii lub elektrostatyki (Tabela 5, Rysunek 7A,B).

W kompleksie urolityna A-CASP3 (PDB ID: 2DKO) przewidywany tryb wiązania charakteryzował się konwencjonalnymi wiązaniami wodorowymi z Gln161 (3,78 i 4,19 Å), Ser120 (3,95 Å) oraz Arg207 (3,05 i 3,77 Å), a dodatkowo był stabilizowany przez oddziaływania π-kation z Arg207, wiązanie wodorowe typu π-donora z Cys163 oraz dodatkowe kontakty π-alkilowe i siły van der Waalsa z udziałem Arg64, Ala162, His121, Ser205 i Trp206 (Tabela 6, Rysunek 8A,B).

W przypadku kompleksu tryptamina-HTR2A (PDB ID: 6A93), przewidywana poza została ustabilizowana przez elektrostatyczny mostek solny między protonowaną aminą tryptaminy a Asp155, konserwowanym asparaginianem helisy transbłonowej 3 (D3.32 w numeracji Ballesteros-Weinstein), który kotwiczy protonowaną aminę ligandów aminergicznych w receptorach serotoninowych i pokrewnych^41,42,43, wraz z wiązaniami wodorowymi z Thr160 i Ser159, kontaktami aromatycznymi z Phe340 i Trp336 oraz oddziaływaniami π-Alkyl z Val156 i Ile163. Dodatkowe kontakty van der Waalsa z Tyr370, Phe339, Ser242, Phe243, Phe332 i Leu123 potwierdziły wzorzec wiązania w kieszeni ortosterycznej (Tabela 7, Rycina 9A,B).

Walidacja protokołu dokowania

Aby ocenić niezawodność protokołu dokowania, przeprowadzono dwa komplementarne eksperymenty kontrolne. W przypadku kontroli redokowania (pozytywnych) ligandy kokrystaliczne zostały wyekstrahowane z ich referencyjnych struktur rentgenowskich i ponownie zadokowane do ich rodzimych miejsc wiązania. Najwyżej oceniona przewidywana poza dla analogu witaminy D VDX w VDR/1DB1 odbiegała o 0,87 Å od pozycji krystalograficznej, a ligand kokrystaliczny WAY-362450 w FXR/3DCT o 1,79 Å; obie wartości znajdowały się poniżej konwencjonalnego progu akceptacji 2,0 Å, co potwierdza poprawność geometryczną protokołu dokowania dla tych układów receptorowych (Rysunek 10A,B). W przypadku kontroli cross-dokowania (negatywnych) kwas litocholowy zadokowano do kaspazy-3 (2DKO), która jest cysteinową proteazą, dla której nie jest on znanym ligandem, co dało przewidywany wynik (−8,3 kcal/mol) słabszy o 1,7 kcal/mol niż w przypadku jego kognatywnego celu VDR (−10,0 kcal/mol), co jest zgodne z przewidywaną selektywnością miejsca wiązania. Zadokowanie tryptaminy do VDR dało przewidywany wynik −6,4 kcal/mol w porównaniu z −7,1 kcal/mol dla kognatywnego celu HTR2A; różnica 0,7 kcal/mol mieści się w raportowanej niepewności wyników dokowania molekularnego ligandów będących metabolitami do białek docelowych, co wskazuje zatem na jedynie umiarkowaną przewidywaną selektywność dla tego mniejszego ligandu (Rysunek 10C). Wspólnie kontrole te wskazują, że protokół dokowania odtwarza znane geometrie wiązania i rozróżnia pary kognatywne od niekognatywnych w badanych warunkach, pozostając przy tym przewidywaniami obliczeniowymi, które nie zastępują eksperymentalnych pomiarów powinowactwa (Tabela 8).

Symulacja dynamiki molekularnej

Symulacje dynamiki molekularnej dla pięciu priorytetowych kompleksów przeprowadzono w ramach trajektorii produkcyjnych o czasie trwania 200 ns. Cztery kompleksy rozpuszczalne i receptory jądrowe symulowano w jawnie określonym rozpuszczalniku wodnym, natomiast kompleks tryptaminy-HTR2A symulowano w jawnie określonej dwuwarstwie lipidowej POPC, aby zapewnić fizjologicznie odpowiednie środowisko błonowe dla tego receptora sprzężonego z białkiem G. Analizy pozwoliły na przetestowanie stabilności dynamicznej zadokowanych póz w warunkach zależnych od czasu oraz na porównanie względnego zachowania strukturalnego poszczególnych kompleksów (Tabela 9).

Profil RMSD kompleksu LCA-VDR/1DB1 wykazał krótki okres ekwilibracji w ciągu pierwszych 10 ns, po którym nastąpiło stabilne plateau z fluktuacjami głównie w zakresie 0,20–0,28 nm (Ryc. 11A). Wartości RMSF były niskie, a fluktuacje szkieletu wynosiły < 0,15 nm dla większości reszt (Ryc. 11B). Analiza wiązań wodorowych wykazała trwałą sieć od 2 do 5 wiązań wodorowych, z okazjonalnymi wzrostami do 7 (Ryc. 11C). Promień żyracji (Rg) utrzymywał się w zakresie 1,25–1,75 nm, a dostępna dla rozpuszczalnika powierzchnia cząsteczki (SASA) utrzymywała się na poziomie około 130 nm² (Ryc. 11D,E).

Kompleks urolithyna A-CASP3/2DKO wykazywał większą aktywność dynamiczną. Wartość RMSD początkowo wzrosła, a następnie oscylowała między 0,4 a 0,7 nm, z krótkotrwałym zdarzeniem o wysokim odchyleniu około 165 ns (Rysunek 12A). Analiza RMSF wykazała wysoką mobilność na poziomie reszt, przy czym największe fluktuacje wystąpiły w regionie elastycznej pętli wokół reszty 175 (Rysunek 12B). Analiza wiązań wodorowych ujawniła początkową rozległą sieć około 2–5 wiązań w ciągu pierwszych 30–40 ns, po której następowało w większości od 0 do 2 przerywanych wiązań (Rysunek 12C). Odpowiednie profile promienia żyracji i SASA przedstawiono na Rysunku 12D,E.

W przypadku systemów kwasów żółciowych NR1H4/FXR (3DCT), profil RMSD szkieletu białkowego utrzymywał się w stosunkowo wąskim zakresie przez większość trajektorii (Rycina 13A), podczas gdy profil RMSF wykazał niższą ruchliwość w regionach rdzeniowych i większe fluktuacje w regionach elastycznych (Rycina 13B). Kompleks LCA-3DCT utrzymywał około trzech do czterech trwałych wiązań wodorowych przez całą trajektorię, natomiast kompleks UDCA-3DCT wykazał większe fluktuacje wiązań wodorowych oraz zmniejszenie liczby wiązań po około 125 ns. Profile promienia żyroskopowego dla systemów związanych z LCA i UDCA przedstawiono odpowiednio na Rycini 13C,D, a odpowiadające im profile SASA przedstawiono na Rycini 13E,F.

Dynamika molekularna błony kompleksu tryptamina-HTR2A

Kompleks tryptamina-HTR2A/6A93 symulowano przez 200 ns w jawnej dwuwarstwie lipidowej POPC składającej się z 258 cząsteczek lipidów, przy zastosowaniu jawnego trzymiejscowego modelu wody i 0,15 M NaCl, co dało całkowitą wielkość układu wynoszącą około 100 925 atomów^33,44,45. Receptor pozostał stabilnie osadzony w dwuwarstwie przez cały czas trwania trajektorii (Rycina 14). RMSD szkieletu wzrosło z około 0,10 nm do stabilnego plateau w przedziale 0,15–0,20 nm w ciągu pierwszych 100 ns i pozostało stabilne w dalszym czasie, przy czym wszystkie wartości były niższe niż 0,25 nm, co wskazuje, że receptor zachował stabilną konformację w środowisku błonowym bez globalnego rozfałdowania (Rycina 15A). RMSF dla poszczególnych reszt wykazało niskie fluktuacje w rdzeniu helis transbłonowych z oczekiwaną wyższą mobilnością w regionach pętli i końców, co jest zgodne z typową elastycznością GPCR (Rycina 15B). Promień gyracji był ściśle ograniczony do przedziału od około 2,06 do 2,12 nm, a SASA fluktuowała w wąskim pasmie bez progresywnego dryfu, co w obu przypadkach potwierdziło zachowanie zwartego pęczka transbłonowego (Rycina 15C,D).

Wiązania wodorowe między białkiem a ligandem były utrzymywane przez cały przebieg trajektorii (Rysunek 15E), przy znaczących wahaniach liczby wiązań wodorowych w zakresie od 1 do 3. Aby specyficznie ocenić trwałość kluczowej interakcji jonowej, monitorowano przez całą trajektorię minimalną odległość między protonowanym azotem amoniowym tryptaminy a atomami tlenu grupy karboksylowej Asp155 (D3.32). Odległość ta pozostawała ściśle rozłożona wokół średniej wartości 0,270 nm (minimum 0,247 nm, maksimum 0,424 nm), a kontakt mostka solnego (< 0,4 nm) był utrzymywany przez 99,9% czasu symulacji, z jedynie dwoma krótkimi, przejściowymi odchyleniami i bez żadnego zdarzenia trwałej dysocjacji (Rysunek 16). Wyniki te sugerują, że konserwatywna interakcja jonowa Asp155 była wystarczająca, aby ustabilizować tryptaminę w kieszeni ortosterycznej HTR2A podczas symulacji w błonie.

Wolna energia wiązania MM-PBSA i dekompozycja na poszczególne reszty

Przeprowadzono analizę MM-PBSA, aby wprowadzić dodatkową warstwę priorytetyzacji energetycznej dla pięciu kompleksów (Tabela 10). W przypadku czterech kompleksów w środowisku wodnym, dekompozycja na poszczególne reszty pozwoliła zidentyfikować główne wkłady energetyczne dla każdego przewidywanego trybu wiązania. W kompleksie LCA-VDR/1DB1 ligand oraz Gln317 miały korzystny wkład, podczas gdy Trp286 wykazywał wkład niekorzystny. W kompleksie urolithin A-CASP3/2DKO reszty Arg64 i Arg207 wykazały silnie ujemne wkłady na poziomie poszczególnych reszt, co wskazuje na istotną stabilizację polarną lub elektrostatyczną; niemniej jednak odpowiadająca im trajektoria pozostała wysoce dynamiczna, co dowodzi, że same korzystne zależności energetyczne na poziomie reszt nie gwarantują trwałości stabilności kompleksu. W systemach 3DCT wiązanie LCA było napędzane głównie przez Arg331, natomiast wiązanie UDCA angażowało bardziej rozproszoną sieć energetyczną obejmującą Glu326, Asp394, Arg395, Arg441 i Asp470. W czterech systemach wodnych dekompozycja MM-PBSA potwierdziła względną priorytetyzację kompleksów opartych na LCA.

W przypadku kompleksu tryptaminy-HTR2A/6A93 osadzonego w błonie, analiza MM-PBSA została przeprowadzona na podsystemie białko–ligand wyekstrahowanym z trajektorii dwuwarstwy^46,47. Zaobserwowano korzystne wkłady dla liganda oraz Asp155 (D3.32), który był zdecydowanie dominującym stabilizującym wkładem na poziomie resydów, co jest zgodne z oddziaływaniem mostka solnego zidentyfikowanym zarówno w dokowaniu, jak i w analizach odległości w trajektorii. Trp137 wykazał największy niekorzystny wkład na poszczególną resydę spośród otaczających resydów kieszeni ortosterycznej (Ser86, Phe87, Phe133, Phe140, Phe141, Val156, Ser159, Thr160, Ile163, Val167, Tyr171), które wspólnie tworzą sieć kontaktów aromatycznych i polarnych wyściełających kieszeń wiążącą. Wartości te stanowią względne szacunki obliczeniowe służące do priorytetyzacji strukturalnej i nie są eksperymentalnymi powinowactwami wiązania.

Rycina 1: Schemat obliczeniowy priorytetyzacji genów gospodarza powiązanych z metabolitami w IBS-C. Schematyczna reprezentacja ośmioetapowego procesu integrującego wybór metabolitów, przewidywanie celów, analizę różnicowej ekspresji transkrypcyjnej, analizę nakładania się danych, analizę wzbogacenia sieci i szlaków, dokowanie molekularne, symulacje dynamiki molekularnej oraz analizę wolnej energii wiązania MM-PBSA. Kliknij tutaj, aby wyświetlić większą wersję tej ryciny.

Rysunek 2: Analiza różnicowej ekspresji i nakładania się celów metabolitów w błonie śluzowej IBS-C. (A) Wykres wulkaniczny (volcano plot) różnicowej ekspresji genów w zbiorze GSE36701. Niebieskie punkty: geny istotnie obniżone; czerwone punkty: geny istotnie podwyższone; szare punkty: geny nieistotne. Oznaczono wybrane nakładające się geny związane z metabolitami. (B) Diagram Venna przedstawiający nakładanie się 330 unikalnych przewidywanych celów metabolitów i genów o obniżonej ekspresji w GSE36701; wspólnych było 17 genów. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rycina 3: Ocena statystyczna 17 przewidywanych genów docelowych dla metabolitów w odniesieniu do GSE36701. (A) log2 fold change dla każdego z 17 genów, z kolorowaniem według poziomu istotności. (B) Stopień różnicowej ekspresji genów tła w porównaniu do przewidywanych celów, z zastosowaniem dokładnego testu Fishera. (C) Do dokładnego testu Fishera wykorzystano tabelę kontyngencji 2x2. Wszystkie 17 celów było istotnie obniżonych (downregulated); pokrycie interpretuje się jako opisowy wzorzec kierunkowy, a nie jako statystyczne wzbogacenie. Kliknij tutaj, aby wyświetlić większą wersję tej ryciny.

Rycina 4: Konstrukcja złożonej sieci oddziaływań białko-białko i wzbogacenie szlaków w sieciach oddziaływań białko-białko dla nakładających się genów związanych z metabolitami. (A) Sieć 1: połączona meta-sieć wszystkich 17 genów. (B) Sieć 2: specyficzna dla propionianu sieć ośmiu genów (CASR, FFAR2, GCG, GNAQ, GPHN, GPR68, MLN, TBXA2R). (C) Sieć 3: sieć tryptaminy/serotoniny czterech genów (HTR1B, HTR2A, HTR2B, HTR6). Sieci zostały wygenerowane dla Homo sapiens, stosując konstrukcję sieci oddziaływań białko-białko i wzbogacenie szlaków z poziomem ufności ≥ 0.700. Krawędzie reprezentują wsparte adnotacjami powiązania funkcjonalne Kliknij tutaj, aby wyświetlić większą wersję tej ryciny.

Rysunek 5: Trójwymiarowa i dwuwymiarowa reprezentacja strukturalna kwasu litocholowego w kompleksie z VDR (PDB ID: 1DB1). (A) Trójwymiarowa reprezentacja powierzchniowa i schematyczna (cartoon), na której kwas litocholowy przedstawiono w formie sfer. (B) Dwuwymiarowa mapa oddziaływań wykazująca wiązanie wodorowe Ser278 oraz otaczające kontakty hydrofobowe i siły van der Waalsa. Proszę kliknąć tutaj, aby wyświetlić większą wersję tego rysunku.

Rysunek 6. Trójwymiarowa i dwuwymiarowa reprezentacja strukturalna kwasu litocholowego w kompleksie z NR1H4/FXR (PDB ID: 3DCT). (A) Trójwymiarowa reprezentacja powierzchniowa i schematyczna. (B) Dwuwymiarowa mapa oddziaływań pokazująca wiązania wodorowe z His294 i Ile335, oddziaływanie π-Sigma oraz kontakty otaczające. Proszę kliknąć tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rycina 7: Trójwymiarowa i dwuwymiarowa reprezentacja strukturalna kwasu ursodeoksicholowego w kompleksie z NR1H4/FXR (PDB ID: 3DCT). (A) Trójwymiarowa reprezentacja powierzchniowa i schematyczna. (B) Dwuwymiarowa mapa oddziaływań wykazująca wiązania wodorowe z His447 i Gly322, oddziaływanie π-anionowe z Val325, wiązanie węgiel-wodór z Trp469 oraz niekorzystne kontakty donor-donor z Arg395 i Gln396. Kliknij tutaj, aby zobaczyć powiększoną wersję tej ryciny.

Rycina 8: Trójwymiarowa i dwuwymiarowa reprezentacja strukturalna urolityny A w kompleksie z CASP3 (PDB ID: 2DKO). (A) Trójwymiarowa reprezentacja powierzchniowa i schematyczna. (B) Dwuwymiarowa mapa oddziaływań wykazująca wiązania wodorowe z Gln161, Ser120 i Arg207, oddziaływania π-kationowe z Arg207, wiązanie wodorowe π-donorowe z Cys163 oraz kontakty otoczenia. Kliknij tutaj, aby wyświetlić większą wersję tej ryciny.

Rysunek 9: Trójwymiarowa i dwuwymiarowa reprezentacja strukturalna tryptaminy w kompleksie z HTR2A (PDB ID: 6A93). (A) Trójwymiarowa reprezentacja powierzchniowa i schematyczna (cartoon) wygenerowana w programie do trójwymiarowej wizualizacji cząsteczkowej. (B) Dwuwymiarowa mapa oddziaływań wygenerowana w narzędziu do wizualizacji cząsteczkowej i tworzenia dwuwymiarowych diagramów oddziaływań, ilustrująca mostek solny Asp155 oraz dodatkowe oddziaływania w miejscu wiązania. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rysunek 10: Walidacja protokołu dokowania. (A,B) Redokowanie ligandów kokrystalizowanych do VDR/1DB1 (RMSD 0,87 Å) oraz FXR/3DCT (RMSD 1,79 Å); nałożenie póz krystalograficznych i redokowanych, obie poniżej progu akceptacji 2,0 Å. (C) Selektywność dokowania krzyżowego: wyniki Vina dla ligandów kognatywnych i niekognatywnych w przypadku kwasu litocholowego i tryptaminy. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rycina 11. Analiza trajektorii dynamiki molekularnej kompleksu LCA-VDR/1DB1 w czasie 200 ns. (A) Profil RMSD. (B) Profil RMSF. (C) Liczba wiązań wodorowych. (D) Profil promienia żyroskopowego. (E) Profil SASA. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rysunek 12: Analiza trajektorii dynamiki molekularnej kompleksu urolityny A-CASP3/2DKO w czasie 200 ns. (A) Profil RMSD wykazujący szerokie fluktuacje konformacyjne i przejściowy zdarzenie o wysokim odchyleniu w okolicach 165 ns. (B) Profil RMSF wykazujący wyraźną elastyczność na poziomie reszt w pobliżu reszty 175. (C) Liczba wiązań wodorowych. (D) Profil promienia gyracji. (E) Profil SASA. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Rycina 13: Analiza trajektorii dynamiki molekularnej układów kwasów żółciowych NR1H4/FXR (3DCT) w czasie 200 ns. (A) Profil RMSD szkieletu białkowego dla kompleksu 3DCT. (B) Profil RMSF szkieletu białkowego. (C) Profil promienia żyroskopowego dla 3DCT-LCA. (D) Profil promienia żyroskopowego dla 3DCT-UDCA. (E) Profil SASA dla 3DCT-LCA. (F) Profil SASA dla 3DCT-UDCA. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rycina 14: Kompleks tryptaminy-HTR2A osadzony w jawnej dwuwarstwie lipidowej POPC. Receptor przedstawiono schematycznie jako strukturę przechodzącą przez dwuwarstwę, lipidy POPC jako linie z zaznaczonymi fosforanowymi grupami głowowymi, a tryptaminę w kieszeni ortosterycznej. Wodę przedstawiono powyżej i poniżej błony. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rycina 15: Analiza trajektorii dynamiki molekularnej kompleksu tryptaminy-HTR2A/6A93 przez 200 ns w jawnej warstwie lipidowej POPC. (A) Profil RMSD szkieletu białkowego. (B) Profil RMSF na poszczególne reszty. (C) Profil promienia gyracji. (D) Profil SASA. (E) Liczba wiązań wodorowych białko-ligand. Kliknij tutaj, aby zobaczyć powiększoną wersję tej ryciny.

Rysunek 16: Trwałość oddziaływania jonowego między tryptaminą a Asp155 (D3.32) w ciągu 200 ns trajektorii w błonie. Wykres przedstawia minimalną odległość między azotem amoniowym tryptaminy a atomami tlenu grupy karboksylowej Asp155 w funkcji czasu; linia przerywana oznacza próg kontaktu mostka solnego wynoszący 0.4 nm. Kontakt był utrzymywany przez 99.9% czasu symulacji. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Symbol genu	Metabolit(y) pochodzenia	Kategoria funkcjonalna	log2FC	FDR (skorygowana wartość P)	Poziom istotności
GCG	Propionian	Białko związane z hormonem peptydowym	−1.342	1.97e−7	FDR <0.001 & |logFC > 1
HDAC3	Maślan	Enzym	−1.234	2.44e−6	FDR <0.001 & |logFC| > 1
CASP3	Urolityna A	Enzym	−1.198	6.66e−7	FDR <0.001 & |logFC| > 1
GPR68	Propionian	Receptor błonowy	−1.137	4.35e−6	FDR <0.001 & |logFC| > 1
GNAQ	Propionian	Wewnątrzkomórkowe białko sygnałowe	−1.122	1.05e−6	FDR <0.001 & |logFC| > 1
GPHN	Propionian	Inne białko wewnątrzkomórkowe	−1.109	1.13e−6	FDR <0.001 & |logFC| > 1
TBXA2R	Propionian	Receptor błonowy	−1.104	4.04e−7	FDR <0.001 & |logFC| > 1
HTR6	Tryptamina	Receptor błonowy	−0.967	2.17e−5	FDR <0.001
VDR	Kwas litocholowy	Receptor jądrowy	−0.942	5.73e−7	FDR <0.001
HTR2A	Tryptamina	Receptor błonowy	−0.937	4.99e−6	FDR <0.001
FFAR2	Propionian	Receptor błonowy	−0.889	1.44e−4	FDR <0.001
NR1H4	Kwas litocholowy / kwas ursodeoksycholowy	Receptor jądrowy	−0.861	3.68e−6	FDR <0.001
HTR2B	Tryptamina	Receptor błonowy	−0.702	1.29e−4	FDR <0.001
MLN	Propionian	Białko związane z hormonem peptydowym	−0.605	7.39e−5	FDR <0.001
KYAT1	Kwas indolowy-3-mlekowy / kwas indolowy-3-propionowy	Enzym	−0.530	3.61e−4	FDR <0.001
CASR	Propionian	Receptor błonowy	−0.483	4.05e−4	FDR <0.001
HTR1B	Tryptamina	Receptor błonowy	−0.455	3.18e−2	FDR <0.05

Tabela 1: Przewidywane geny docelowe związane z metabolitami pokrywające się z genami o różnej ekspresji w zbiorze danych błony śluzowej odbytnicy IBS-C. Wszystkie wymienione pokrywające się geny były wyciszone. Tabela 1 została przesłana oddzielnie jako arkusz kalkulacyjny i wymienia dla każdego celu: metabolit(y) pochodzenia, kategorię funkcjonalną, źródło przewidywania celu (przewidywanie celu oddziaływania chemiczno-białkowego, program dokowania molekularnego lub oba), łączny wynik oddziaływania z przewidywania celu oddziaływania chemiczno-białkowego oraz prawdopodobieństwo z programu dokowania molekularnego, a jeśli dostępne – poziom przewidywania, log2 fold change oraz FDR z poziomem istotności ekspresji. Źródło: Wartości ekspresji genów uzyskano z tabeli różnicowej ekspresji GSE36701 z zapadniętymi genami (sonda o najniższym FDR na gen). Źródło przewidywania celu oraz wartości ufności opracowano na podstawie wyników przewidywania celu oddziaływania chemiczno-białkowego oraz programu dokowania molekularnego, stosując progi: łączny wynik oddziaływania z przewidywania celu oddziaływania chemiczno-białkowego ≥ 0,700 oraz prawdopodobieństwo z programu dokowania molekularnego ≥ 0,70. Wyniki przewidywania celu oddziaływania chemiczno-białkowego są wynikami łącznymi w skali 0–1; STP oznacza prawdopodobieństwo z programu dokowania molekularnego. Poziom 1 = ścisłe wsparcie przewidywania celu oddziaływania chemiczno-białkowego; Poziom 1+ = ścisłe wsparcie przewidywania celu oddziaływania chemiczno-białkowego wsparte dodatkowo przez program dokowania molekularnego.

Złożony	Białko (PDB ID)	Ligand	Wynik Vina (kcal/mol)	Rozmiar kawerny (A^3)	Centrum siatki X,Y,Z (A)	Pole wyszukiwania (A)
LCA-VDR	VDR (1DB1)	Kwas litocholowy	−10.0	2055	10, 19, 33	25 x 25 x 25
LCA-NR1H4/FXR	NR1H4/FXR (3DCT)	Kwas litocholowy	−9.9	3395	137, 31, 78	25 x 25 x 25
UDCA-NR1H4/FXR	NR1H4/FXR (3DCT)	Kwas ursodeoksycholowy	−9.4	3395	137, 31, 78	25 x 25 x 25
Urolityna A-CASP3	CASP3 (2DKO)	Urolityna A	−7.1	233	37, 34, 32	25 x 25 x 25
Tryptamina-HTR2A	HTR2A (6A93)	Tryptamina	−7.1	3238	12, −1, 61	25 x 25 x 25

Tabela 2: Wyniki dokowania molekularnego: najwyżej ocenione wyniki dokowania molekularnego ligandów-metabolitów do białek docelowych oraz parametry kawitacji dla pięciu priorytetowych kompleksów białko-ligand. Rozmiar kawitacji podano w ų. Źródło: Docking_Validation/Results/Docking_Validation_Results.xlsx, arkusz 'Original_Docking_Scores'. Dokowanie molekularne ligandów-metabolitów do białek docelowych; exhaustiveness = 8, seed = 42 (stałe), num_modes = 9 dla wszystkich kompleksów; przedstawiono najwyżej ocenioną pozę (mode 1).

Typ oddziaływania	Pozostałość (pozostałości)	Odległość (A)	Uwagi
Konwencjonalne wiązanie wodorowe	Ser278	4.29	Tlen karboksylowy LCA
Kontakt hydrofobowy / Pi-Alkyl	Leu230, Val234, Trp286, Val300, His305, Tyr295, Leu233, His397	-
Kontakt Van der Waalsa	Met272, Leu313, Ile271, Ile268, Leu309, Phe422, Val418, Ala231, Ala303, Cys288, Ser275, Phe150	-

Tabela 3: Sposoby wiązania wygenerowane dla dokowania kwasu litocholowego z VDR (PDB ID: 1DB1). Źródło: narzędzie do wizualizacji molekularnej i dwuwymiarowych diagramów interakcji – dwuwymiarowe diagramy interakcji ligand-pozostałość, zgodnie z opisem w sekcji Wyniki (Dokowanie molekularne) w manuskrypcie. „-” oznacza, że wartość odległości nie została indywidualnie podana dla danego kontaktu.

Typ oddziaływania	Pozostałość(i)	Odległość (A)	Uwagi
Wiązanie wodorowe	His294	-
Wiązanie wodorowe	Ile335	-
Oddziaływanie pi-sigma	His294	-
Alkilowy / Pi-alkilowy (hydrofobowy)	Met290, Met328, Ala291, Leu287, Ile352, His447	-
Kontakt van der Waalsa	Dodatkowe reszty kieszeni (nieokreślone indywidualnie w źródle)	-	Umożliwia dopasowanie szkieletu steroidowego

Tabela 4: Sposoby wiązania wygenerowane dla dokowania kwasu litocholowego z NR1H4/FXR (PDB ID: 3DCT).

Źródło: narzędzie do wizualizacji cząsteczkowej i dwuwymiarowych diagramów oddziaływań; dwuwymiarowe diagramy oddziaływań ligand-pozostałość, zgodnie z treścią sekcji Wyniki (Dokowanie molekularne) w manuskrypcie. Znak „-” oznacza, że dla danego kontaktu nie podano indywidualnej wartości odległości.

Typ oddziaływania	Pozostałość(ie)	Odległość (A)	Uwagi
Konwencjonalne wiązanie wodorowe	His447	3.66
Wiązanie wodorowe	Gly322	3.46
Oddziaływanie pi-anionowe	Val325	4.96
Wiązanie węgiel-wodór	Trp469	4.51
Niekorzystny kontakt donor-donor	Arg395	3.89
Niekorzystny kontakt donor-donor	Gln396	3.40

Tabela 5: Tryby wiązania wygenerowane dla dokowania kwasu ursodeoksycholowego z NR1H4/FXR (PDB ID: 3DCT). Źródło: narzędzie do wizualizacji molekularnej i dwuwymiarowych diagramów oddziaływań, dwuwymiarowe diagramy oddziaływań ligand-pozostałość, jak opisano w sekcji Wyniki (Dokowanie molekularne) manuskryptu. „-” oznacza, że dla danego kontaktu nie podano indywidualnej wartości odległości.

Typ oddziaływania	Pozostałość(y)	Odległość (A)	Uwagi
Konwencjonalne wiązanie wodorowe	Gln161	3.78
Konwencjonalne wiązanie wodorowe	Gln161	4.19	drugi kontakt
Konwencjonalne wiązanie wodorowe	Ser120	3.95
Konwencjonalne wiązanie wodorowe	Arg207	3.05
Konwencjonalne wiązanie wodorowe	Arg207	3.77	drugi kontakt
Oddziaływanie Pi-Kation	Arg207	-
Wiązanie wodorowe Pi-Donor	Cys163	-
Kontakt Pi-Alkyl / van der Waalsa	Arg64, Ala162, His121, Ser205, Trp206	-

Tabela 6: Tryby wiązania wygenerowane podczas dokowania urolityny A do CASP3 (PDB ID: 2DKO).Źródło: narzędzie do wizualizacji molekularnej i dwuwymiarowych diagramów interakcji; dwuwymiarowe diagramy interakcji ligand-pozostałość, zgodnie z opisem w sekcji Wyniki (Dokowanie molekularne) rękopisu. „-” oznacza, że wartość odległości nie została indywidualnie podana dla danego kontaktu.

Typ oddziaływania	Pozostałości	Odległość (A)	Uwagi
Elektrostatyczny mostek solny	Asp155 (D3.32)	-	protonowana amina tryptaminy
Wiązanie wodorowe	Thr160	-
Wiązanie wodorowe	Ser159	-
Kontakt aromatyczny	Phe340, Trp336	-
Oddziaływanie pi-alkilowe	Val156, Ile163	-
Kontakt van der Waalsa	Tyr370, Phe339, Ser242, Phe243, Phe332, Leu123	-

Tabela 7: Tryby wiązania wygenerowane dla dokowania tryptaminy do HTR2A (PDB ID: 6A93). Źródło: narzędzie do wizualizacji molekularnej i dwuwymiarowych diagramów oddziaływań, dwuwymiarowe diagramy oddziaływań ligand-pozostałość, przedstawione w sekcji Wyniki (Dokowanie molekularne) w manuskrypcie. Symbol „-” oznacza, że wartość odległości dla danego kontaktu nie została podana indywidualnie.

(A) Walidacja redockingu (kontrole pozytywne)
PDB ID	Białko	Ligand kokrystaliczny	Wynik Vina (kcal/mol)	RMSD (A)	Próg (A)	Wynik
1DB1	VDR	VDX (analog witaminy D)	−13.0	0.87	2.0	ZALICZONO
3DCT	FXR	WAY-362450 (064)	−11.9	1.79	2.0	ZALICZONO
(B) Walidacja cross-dockingu (kontrole negatywne)
Ligand	Cel kognatywny (PDB)	Wynik kognatywny (kcal/mol)	Cel niekognatywny (PDB)	Wynik niekognatywny (kcal/mol)	Delta (kcal/mol)	Selektywność
Kwas litocholowy	VDR (1DB1)	−10.0	CASP3 (2DKO)	−8.3	1.7	Potwierdzona
Tryptamina	HTR2A (6A93)	−7.1	VDR (1DB1)	−6.4	0.7	Umiarkowana (w granicach niepewności Vina +/−0.5–1.0)

Tabela 8: Wyniki walidacji protokołu dokowania: wartości RMSD dla redokowania (kontrole dodatnie) oraz wyniki cross-dockingu (kontrole ujemne).Źródło: Docking_Validation/Results/Docking_Validation_Results.xlsx oraz Docking_Validation/Logs/*.log (dokowanie molekularne ligandów metabolitów do białek docelowych, exhaustiveness = 8, seed = 42, pole 25 Å × 25 Å × 25 Å). RMSD obliczone na podstawie dopasowania ciężkich atomów i nazw atomów (bez superpozycji).

Kompleks	RMSD (nm), średnia + / –SD (zakres)	Rg (nm), średnia + / – SD (zakres)	SASA (nm^2), średnia + / – SD (zakres)	Wiązania wodorowe, średnia + / –SD (zakres)	RMSF (nm), średnia (max)
LCA-VDR/1DB1	0.230 + / – 0.025 (0.167–0.296)	1.889 + / − 0.009 (1.863–1.919)	130.4 + / − 2.3 (122.3–137.4)	1.9 + / − 0.9 (0–7)	0.093 (max 0.600 w pozostałości 120)
LCA-NR1H4/FXR/3DCT	0.190 + / – 0.020 (0.135–0.281)	1.824 + / − 0.008 (1.804–1.849)	129.7 + / − 2.3 (121.9–138.1)	3.8 + / − 0.7 (1–6)	0.113 (max 0.298)
UDCA-NR1H4/FXR/3DCT	0.190 + / – 0.020 (0.135–0.281)	1.834 + / − 0.013 (1.809–1.921)	131.0 + / −3.4 (121.6–143.5)	1.1 + / − 1.1 (0–5)	0.113 (max 0.298)
Urolithin A-CASP3/2DKO	0.521 + / – 0.058 (0.244–0.755)	1.892 + / − 0.024 (1.839–1.984)	134.9 + / − 3.0 (126.4–146.4)	0.6 + / − 0.7 (0–3)	1.172 (max 2.532 w pozostałości 175)
Tryptamine-HTR2A/6A93 (błona)	0.177 + / –0.017 (0.131–0.227)	2.089 + / − 0.007 (2.070–2.116)	165.1 + / − 2.7 (156.–172.7)	1.7 + / − 0.7 (0–4)	0.090 (max 0.319)

Tabela 9: Podsumowanie przebiegu symulacji dynamiki molekularnej trwającej 200 ns dla pięciu priorytetowych kompleksów białko-ligand, w tym układu tryptaminy-HTR2A osadzonego w błonie.Źródło: pliki narzędzi do analizy trajektorii dynamiki molekularnej (.xvg) — gmx rms, gmx gyrate, gmx sasa, gmx hbond, gmx rmsf — obliczone dla ostatnich 150 ns (50–200 ns) każdego 200 ns przebiegu produkcyjnego, zgodnie z krokiem 8.8 protokołu. RMSD/Rg dopasowane do szkieletu białkowego; promień sondy SASA 0,14 nm; odcięcie donor-akceptor wiązania wodorowego 0,35 nm / 30 °. LCA-3DCT i UDCA-3DCT współdzielą jedną trajektorię szkieletu białkowego (RMSD, RMSF) przy specyficznych dla ligandów wartościach Rg/SASA/wiązań wodorowych.

Tryptamina-HTR2A/6A93 (błona) — ilościowy rozkład na poszczególne reszty
Reszta	Całkowity wkład ddG (kcal/mol), średnia + / − SD	Kierunek
Asp155 (D3.32)	−89.94 + / − 6.81	Stabilizujący (dominujący)
Tryptamina (ligand)	−13.01 + / − 6.22	Stabilizujący
Tyr171	13.62 + / − 4.54	Destabilizujący
Val167	23.32 + / − 3.96	Destabilizujący
Val156	20.03 + / − 3.81	Destabilizujący
Thr160	4.86 + / − 3.64	Destabilizujący
Ser159	24.16 + / − 3.48	Destabilizujący
Ser86	24.48 + / − 3.65	Destabilizujący
Phe87	35.18 + / − 4.04	Destabilizujący
Phe133	32.80 + / −3.70	Destabilizujący
Phe140	30.63 + / − 3.84	Destabilizujący
Phe141	35.25 + / − 3.55	Destabilizujący
Ile163	27.64 + / − 3.71	Destabilizujący
Trp137	53.77 + / − 4.32	Destabilizujący (największy niekorzystny)
Pozostałe cztery kompleksy — reszty zidentyfikowane w rozkładzie na poszczególne reszty (jakościowy)
Kompleks	Reszta	Kierunek
LCA-VDR/1DB1	Ligand (LCA)	Korzystny
LCA-VDR/1DB1	Gln317	Korzystny
LCA-VDR/1DB1	Trp286	Niekorzystny
LCA-NR1H4/FXR/3DCT	Arg331	Korzystny (dominujący)
UDCA-NR1H4/FXR/3DCT	Glu326	Mieszana/rozproszona sieć
UDCA-NR1H4/FXR/3DCT	Asp394	Mieszana/rozproszona sieć
UDCA-NR1H4/FXR/3DCT	Arg395	Mieszana/rozproszona sieć
UDCA-NR1H4/FXR/3DCT	Arg441	Mieszana/rozproszona sieć
UDCA-NR1H4/FXR/3DCT	Asp470	Mieszana/rozproszona sieć
Urolityna A-CASP3/2DKO	Arg64	Silnie korzystny (polarny/elektrostatyczny)
Urolityna A-CASP3/2DKO	Arg207	Silnie korzystny (polarny/elektrostatyczny)

Tabela 10: Dekompozycja MM-PBSA na poszczególne reszty KRÓTKI STRESZCZENIE: reszty stabilizujące i destabilizujące (wkład bezwzględny ≥ 0.5 kcal mol⁻¹) dla każdego z pięciu priorytetowych kompleksów białko-ligand, w tym systemu tryptaminy-HTR2A osadzonego w błonie. Źródło: Membrane Simulation/03_MMPBSA/results/FINAL_DECOMP_MMPBSA.dat (narzędzie do obliczeń energii wiązania metodą mechaniki molekularnej/ciągłego rozpuszczalnika – dekompozycja na poszczególne reszty metodą Generalized Born (GB), „Complex: Total Energy Decomposition”). Numery reszt przeliczono z numeracji wewnętrznej systemu zbudowanego w CHARMM-GUI (przesunięcie +68) na oryginalną numerację PDB 6A93 stosowaną w pozostałych częściach niniejszego manuskryptu.

Źródło: Dane z poprzednich symulacji MD/1DB1,2KD0, LCA & UDCA_3DCT}/mmpbsa_*/Decomposition_NORMAL_GB_Complex_TDC*.svg oraz wyniki w manuskrypcie (wolna energia wiązania MM-PBSA i dekompozycja na poszczególne reszty). Dla tych czterech kompleksów w katalogu projektu nie ma numerycznych wyników dekompozycji na reszty w formacie .dat/.csv (dostępne są jedynie wyrenderowane wykresy SVG z tekstem w ścieżkach wektorowych, który nie jest możliwy do wyodrębnienia maszynowo); zgodnie z treścią manuskryptu zgłoszono jedynie tożsamość reszt oraz kierunek korzystny/niekorzystny. Dokładne wkłady w kcal/mol dla tych czterech kompleksów nie są dostępne w repozytorium źródłowym.

Dyskusja

Niniejsze eksploracyjne badanie obliczeniowe demonstruje zintegrowany, powtarzalny schemat postępowania służący do priorytetyzacji genów gospodarza związanych z metabolitami drobnoustrojów oraz kompleksów białko-ligand, zastosowany w tym przypadku do publicznego zbioru danych transkrypomicznych błony śluzowej odbytnicy w IBS-C. Przy użyciu tego schematu podzbiór przewidywanych genów związanych z metabolitami drobnoustrojów pokrywał się z genami konsekwentnie obniżonymi w zbiorze danych, grupując się w obrębie szlaków związanych z GPCR, układem serotoninergicznym, sygnalizacją wapniową, neuroaktywnymi receptorami ligandów oraz receptorami jądrowymi – systemami coraz częściej wiązanymi z komunikacją między mikrobiotą a gospodarzem^48,49. Wyniki te należy interpretować ściśle jako generujące hipotezy: analiza nie mierzy stężeń metabolitów drobnoustrojów, ilości białek receptorowych, wiązania ligandów, aktywacji receptorów, sygnalizacji w dół, motoryki, sekrecji, reakcji bólowych ani wyników klinicznych. Najsilniejszym uzasadnionym wnioskiem jest stwierdzenie, że zidentyfikowane geny i szlaki są kandydatami do walidacji eksperymentalnej, a nie potwierdzonymi mechanizmami choroby.

Kluczowe znaczenie niniejszego protokołu, w odniesieniu do wcześniejszych prac badających pojedyncze pary metabolit-receptor w izolacji, polega na integracji przewidywania celu, publicznych danych transkrypcyjnych, analizy sieci, dokowania z kontrolami walidacyjnymi, dynamiki molekularnej oraz MM-PBSA w jeden sekwencyjny potok priorytetyzacji. Każdy etap zawęża i kontekstualizuje zestaw kandydatów wygenerowany na poprzednim etapie, a właśnie to sekwencyjne filtrowanie sprawia, że końcowa lista kandydatów jest możliwa do zweryfikowania eksperymentalnie. Moduł GPCR skoncentrowany na GNAQ oraz moduł związany z receptorem serotoninowym zidentyfikowane w tej pracy są biologicznie prawdopodobne, biorąc pod uwagę rolę sygnalizacji Gq w aktywacji fosfolipazy C, produkcji inozytolo-1,4,5-trisfosforanu, mobilizacji wapnia, sekrecji i funkcji enteroendokrynnych, a także uznaną rolę sygnalizacji krótkołańcuchowych kwasów tłuszczowych i pochodnych tryptofanu w homeostazie śluzówki oraz sygnalizacji serotonergicznej w motoryce przewodu pokarmowego, sekrecji, wrażliwości trzewnej i komunikacji oś gut-brain^10,18,50,51.

Kluczową cechą metodologiczną niniejszego badania jest traktowanie receptora HTR2A wbudowanego w błonę. Ponieważ symulacja w fazie rozpuszczalnej nie może odtworzyć środowiska lipidowego, które determinuje zachowanie konformacyjne receptora sprzężonego z białkiem G, kompleks tryptaminy z HTR2A symulowano w jawnej dwuwarstwie POPC. W tym środowisku błonowym receptor pozostał stabilny strukturalnie przez całą 200 ns trajektorię, a mostek solny między grupą amoniową tryptaminy a Asp155 (D3.32) był utrzymywany przez niemal cały czas symulacji. Fakt, że trzy niezależne linie dowodowe – poza dokowania, trwały dystans kontaktu w trakcie trajektorii oraz dominujący wkład MM-PBSA na poszczególne reszty – wskazują na tę samą konserwowaną oddziaływanie z D3.32, zapewnia spójność wewnętrzną przewidzianego sposobu wiązania tryptaminy, który odzwierciedla kanoniczną geometrię wiązania ligandów aminergicznych w receptorach serotoninowych.

Podczas odtwarzania tego schematu pracy należy wziąć pod uwagę pewne kwestie metodologiczne. Błędy w strukturze kanonicznej lub związki oznaczone jako PAINS (Pan-assay interference compounds) przenoszą się na etap przewidywania celów i dokowania, co wymaga dokładnego wyboru metabolitów i kurateli chemoinformatycznej. Należy zminimalizować zestawy celów wynikające z szumu, konsekwentnie stosując kryteria ufności (przewidywanie celów oddziaływania chemiczno-białkowego: ≥0.700; program dokowania molekularnego: ≥0.70; konstrukcja sieci oddziaływań białko-białko i wzbogacenie szlaków: ≥0.700). Przewidziane cele powinny być grupowane według kategorii funkcjonalnych, aby uniknąć błędnego zakwalifikowania wszystkich genów powiązanych z metabolitami jako receptorów. Precyzyjne przetwarzanie wstępne struktur PDB, minimalizacja energii ligandów oraz rozmieszczenie siatki wokół znanych reszt wiążących są kluczowymi aspektami dokowania, a wprowadzone tutaj kontrole redokowania i dokowania krzyżowego zapewniają obiektywną miarę poprawności metodologii dokowania. Zakres powtarzalności w dynamice molekularnej jest definiowany przez kombinację parametryzacji pola siłowego, odpowiedniej solwatacji lub budowy błony, etapowej ekwilibracji oraz właściwego próbkowania produkcyjnego.

Typowe adaptacje i kroki rozwiązywania problemów obejmują poluzowanie progów, jeśli przewidywanie celu nie zwraca żadnych trafień, sprawdzanie spójności kierunkowej na poziomie sond dla genów z wieloma sondami oraz interpretowanie odizolowanych węzłów w konstrukcji sieci oddziaływań białko-białko i wzbogacaniu szlaków jako zależnych od progu, a nie jako biologicznie nieistotnych. W przypadku receptorów błonowych należy stosować jawne symulacje dwuwarstwy lipidowej zamiast symulacji w środowisku wodnym, co zilustrowano w opisanym tutaj метоdzie HTR2A. W sytuacjach, gdy konieczna jest dekompozycja energii na poszczególne reszty, obliczenia muszą być wykonane przy użyciu silnika umożliwiającego dekompozycję, a raportowana numeracja reszt powinna zostać uzgodniona z natywną numeracją receptora, aby uniknąć niejednoznaczności. Sugerujemy, aby wyniki wzbogacania szlaków traktować przede wszystkim jako kontekst organizacyjny dla listy kandydatów, a nie jako walidację na poziomie szlaków. Mechanicznie wzbogacenie terminów dotyczących GPCR, układu serotonergicznego lub sygnalizacji wapniowej będzie występować zawsze, gdy lista genów zawiera wiele genów receptorów serotoninowych, niezależnie od koregulacji na poziomie białek. Wartości RMSD, Rg i RMSF dla systemu HTR2A osadzonego w błonie powinny być interpretowane z uwzględnieniem dwuwarstwy lipidowej: spadek Rg w późniejszej trajektorii może odzwierciedlać indukowaną przez dwuwarstwę adaptację konformacyjną wiązki transbłonowej, a nie globalne rozfałdowanie, natomiast trwałe wiązania wodorowe ligand-białko należy interpretować łącznie z ogólną stabilnością RMSD.

Ograniczenia niniejszego badania są znaczące i ograniczają interpretację wyników. Badanie oparto na jednym, stosunkowo niewielkim publicznym zbiorze danych, a przeszukiwanie głównych publicznych repozytoriów transkrypcyjnych (internetowe narzędzie do analizy różnicowej ekspresji genów oraz ArrayExpress) nie pozwoliło w momencie analizy na zidentyfikowanie niezależnego zbioru danych transkrypcyjnych błony śluzowej odbytnicy w przebiegu IBS-C o porównywalnym projekcie i platformie, który mógłby służyć jako kohorta replikacyjna. Brak niezależnej replikacji transkrypcyjnej jest istotnym ograniczeniem i żadne stwierdzenie w niniejszym manuskrypcie nie powinno być interpretowane jako zewnętrzna walidacja wyników uzyskanych z pojedynczego zbioru danych. Zbiór danych wykazuje niemal powszechną różnicową ekspresję (około 94,5% genów jest istotnych, z czego ogromna większość jest wyciszona), co sprawia, że konwencjonalne statystyki wzbogacenia są nieinformacyjne i uniemożliwia wyciągnięcie wniosków na temat specyfiki wyciszenia genów docelowych w stosunku do tła genomicznego; w związku z tym nakładka jest raportowana jako opisowy wzorzec kierunkowy, a nie statystyczne wzbogacenie. Masowa transkryptomika błony śluzowej nie pozwala odróżnić rzeczywistej regulacji genów od zmian w składzie komórkowym. Ekspresja mRNA nie determinuje zawartości białka ani odpowiedzi funkcjonalnej. Bazy danych do przewidywania celów obciążone są błędem adnotacji, a wyniki dokowania, MD oraz MM-PBSA zależą od wyboru pola siłowego, parametryzacji ligandu, pozycji startowej, czasu symulacji oraz adekwatności próbkowania. Dokładne identyfikatory wersji/budowy dla niektórych komponentów serwerów WWW i pakietów, w tym usługi parametryzacji ligandów kompatybilnej z CHARMM, CHARMM-GUI oraz środowiska obliczeń statystycznych

Wersje kompilacji pakietów oraz podwersje narzędzia do obliczeń energii wiązania metodą mechaniki molekularnej/ciągłego rozpuszczalnika nie zostały w pełni odzyskanymi z archiwalnego zapisu projektu i powinny być raportowane jako dostępne w oddzielnej Tabeli materiałów. Wartości MM-PBSA są szacunkami względnymi, nie zawierają jawnego składnika entropii konfiguracyjnej i nie powinny być interpretowane jako powinowactwa eksperymentalne. Badaniu brakuje danych metabolomicznych, w związku z czym nie można określić, czy dostępność ligandów jest zmieniona w IBS-C, ani czy zaobserwowane zmiany w ekspresji są przyczynami, konsekwencjami, odpowiedziami kompensacyjnymi czy nieistotnymi korelacjami.

Przyszłe zastosowania tej metody powinny obejmować niezależną replikację transkrypcyjną, ilościową reakcję łańcuchową polimerazy (qPCR) i walidację na poziomie białka, lokalizację typów komórek za pomocą transkryptomiki jednokomórkowej lub przestrzennej, profilowanie metabolomiczne odpowiednich klas metabolitów oraz funkcjonalne testy odpowiedzi na ligand w kolonoidach pochodzących od pacjentów, eksplantach błony śluzowej lub porównywalnych modelach. Porównania z kohortami pacjentów z IBS z przewagą biegunek, IBS mieszanym, nieswoistymi zapaleniami jelit oraz zaparciami niepowiązanymi z IBS^1,2 pomogłyby ustalić specyficzność choroby. W odniesieniu do komponentu strukturalnego, replikacja trajektorii MD, przeprowadzenie analiz wrażliwości z alternatywnymi pozami startowymi oraz pełna dokumentacja depozycji topologii, trajektorii oraz plików wejściowych i wyjściowych MM-PBSA dodatkowo wzmocniłyby powtarzalność. Eksperymentalne testy odpowiedzi na ligand pozostają niezbędne do ustalenia, czy priorytetowe kompleksy są funkcjonalnie istotne; obecne wyniki nie stanowią podstawy do twierdzeń klinicznych lub terapeutycznych.

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
AutoDock Vina	Scripps Research / open source	v1.2.7; https://vina.scripps.edu/	Dokowanie molekularne ligandów metabolitów do białek docelowych.
CGenFF/ParamChem	SilcsBio / University of Maryland	v4.6; https://cgenff.com/	Parametryzacja pola siłowego ligandów dla dynamiki molekularnej.
Pole siłowe CHARMM36m	Twórcy CHARMM / open source	CHARMM36m; https://www.charmm.org/charmm/resources/charmm-force-fields/	Pole siłowe białek stosowane w symulacjach dynamiki molekularnej.
CHARMM-GUI Membrane Builder	CHARMM-GUI / Lehigh University	Serwer WWW; dokładna wersja nieodzyskiwalna; https://www.charmm-gui.org/?doc=input/membrane	Konstrukcja i konfiguracja ekwilibracji jawnego systemu błony POPC.
Discovery Studio Visualizer	BIOVIA (Dassault Systèmes)	2021; https://discover.3ds.com/discovery-studio-visualizer-download	Dwuwymiarowa analiza oddziaływań ligand-pozostałość.
GEO2R	NCBI Gene Expression Omnibus	Narzędzie WWW; dostęp w okresie styczeń-maj 2026; https://www.ncbi.nlm.nih.gov/geo/geo2r/	Analiza różnicowej ekspresji GSE36701.
GeneCards	Instytut Nauki Weizmanna	Baza danych WWW; dostęp w okresie styczeń-maj 2026; https://www.genecards.org/	Weryfikacja symboli i informacji o genach podczas standaryzacji celów.
gmx_MMPBSA	Open source (Valdés-Tresanco i in.)	1.5.x; https://valdes-tresanco-ms.github.io/gmx_MMPBSA/	Szacowanie wolnej energii wiązania MM-PBSA i dekompozycja na poszczególne pozostałości.
GROMACS	Zespół programistów GROMACS / open source	2024.2; https://www.gromacs.org/	Silnik symulacji dynamiki molekularnej.
Zbiór danych transkrypcyjnych GSE36701	NCBI Gene Expression Omnibus	GSE36701; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36701	Publiczny zbiór danych ekspresji śluzówki odbytnicy IBS-C.
Open Babel	Open source	3.2.0; https://openbabel.org/	Konwersja formatów plików chemicznych, generowanie trójwymiarowych ligandów i przygotowanie ligandów.
Baza danych OPM	University of Michigan	Baza danych WWW; dostęp w okresie styczeń-maj 2026; https://opm.phar.umich.edu/	Współrzędne orientacji białek w błonach (Orientation of Proteins in Membranes) użyte do wyrównania HTR2A.
ParmEd	Twórcy ParmEd / open source	4.x; https://parmed.github.io/ParmEd/html/index.html	Repartycja masy wodoru i przetwarzanie topologii symulacji molekularnych.
PyMOL	Schrödinger / open source	2.x; https://www.pymol.org/	Trójwymiarowa wizualizacja strukturalna i ryciny receptor-ligand.
RCSB Protein Data Bank	RCSB PDB	Baza danych WWW; dostęp w okresie styczeń-maj 2026; https://www.rcsb.org/	Źródło eksperymentalnych struktur białek i współrzędnych PDB.
STITCH	Konsorcjum STITCH (EMBL)	v5.0; https://stitch.embl.de/	Przewidywanie celów oddziaływań chemicznych z białkami.
STRING	Konsorcjum STRING / ELIXIR	v12.0; https://version-12-0.string-db.org/	Konstrukcja sieci oddziaływań białko-białko i wzbogacenie szlaków.
SwissADME	SIB Swiss Institute of Bioinformatics / University of Lausanne	Narzędzie WWW; dostęp w okresie styczeń-maj 2026; https://www.swissadme.ch/	Deskryptory chemoinformatyczne, przewidywania farmakokinetyczne i ocena PAINS.
SwissTargetPrediction	SIB Swiss Institute of Bioinformatics / University of Lausanne	Narzędzie WWW; dostęp w okresie styczeń-maj 2026; https://www.swisstargetprediction.ch/	Przewidywanie ludzkich białek docelowych na podstawie ligandów.
UniProt ID Mapping	Konsorcjum UniProt	Usługa WWW; dostęp w okresie styczeń-maj 2026; https://www.uniprot.org/id-mapping	Mapowanie identyfikatorów białek na standaryzowane symbole genów zatwierdzone przez HGNC.
NVIDIA RTX 3080	NVIDIA Corporation	RTX 3080; ≥8 GB VRAM; wersja CUDA/sterownika nie określona w manuskrypcie	Jednostka przetwarzania graficznego obsługująca CUDA, użyta do symulacji dynamiki molekularnej.
GPU kompatybilne z CUDA	NVIDIA Corporation	Wersja zestawu narzędzi CUDA nie określona w manuskrypcie; ≥8 GB VRAM	GPU obsługujące CUDA z ≥8 GB VRAM; stacja robocza wymagała również ≥32 GB RAM i 6-rdzeniowego procesora CPU.
Ubuntu Linux	Canonical Ltd. / open source	22.04 LTS	64-bitowy system operacyjny Linux.
Python 3.9	Python Software Foundation	3.9	Ogólne środowisko programistyczne użyte do skryptowania i analizy przepływu pracy.
Gene Expression Omnibus (GEO)	NCBI / U.S. National Library of Medicine	Publiczne repozytorium WWW; wersja oprogramowania nie określona w manuskrypcie	Publiczne repozytorium danych genomiki funkcjonalnej.
AutoDockTools/MGLTools	Molecular Graphics Laboratory, Scripps Research	1.5.7	Zestaw narzędzi do przygotowania struktur molekularnych i danych wejściowych do dokowania.
Narzędzia analityczne GROMACS	Zespół programistów GROMACS / open source	2024.2	Narzędzia do analizy trajektorii dynamiki molekularnej.
Sześciostopniowy protokół CHARMM-GUI	CHARMM-GUI / Lehigh University	Protokół WWW; dokładna wersja nieodzyskiwalna	Wieloetapowy, internetowy przepływ pracy przygotowania i ekwilibracji systemu molekularnego.
cgenff_charmm2gmx_py3.py	Skrypt konwersyjny open-source; źródło nie określone w manuskrypcie	Wersja nie określona w manuskrypcie	Skrypt konwersyjny topologii pola siłowego.
Python	Python Software Foundation	3.9	Ogólne środowisko programistyczne.
SciPy	Społeczność SciPy / open source	Wersja nie określona w manuskrypcie	Biblioteka obliczeń naukowych.
scipy.stats.fisher_exact	Społeczność SciPy / open source	Wersja SciPy nie określona w manuskrypcie	Implementacja dokładnego testu Fishera’a.
R	R Foundation for Statistical Computing	4.3.x	Środowisko obliczeń statystycznych.
Bioconductor	Projekt Bioconductor / open source	3.18	Framework oprogramowania bioinformatycznego.
limma	Projekt Bioconductor / open source	Wersja nie określona w manuskrypcie	Pakiet do analizy różnicowej ekspresji genów.
Procedura Benjamini–Hochberg	Metoda statystyczna	Nie dotyczy (procedura statystyczna)	Metoda korekty współczynnika fałszywych odkryć (FDR).
NVIDIA RTX 3080	NVIDIA Corporation	RTX 3080; ≥8 GB VRAM; wersja CUDA/sterownika nie określona w manuskrypcie	Jednostka przetwarzania graficznego z co najmniej 8 GB pamięci wideo.
GPU kompatybilne z CUDA	NVIDIA Corporation	Wersja zestawu narzędzi CUDA nie określona w manuskrypcie; ≥8 GB VRAM	Jednostka przetwarzania graficznego wspierająca ogólne obliczenia równoległe.
Ubuntu Linux 22.04 LTS	Canonical Ltd. / open source	22.04 LTS	64-bitowy system operacyjny Linux.
TIP3P	Twórcy pola siłowego CHARMM / open source	TIP3P; brak wersji oprogramowania	Jawny trójpunktowy model wody.
MM/PBSA	Twórcy gmx_MMPBSA / open source	gmx_MMPBSA 1.5.x	Metoda obliczania energii wiązania mechanika molekularnego/powierzchni Poisson–Boltzmana.