Protokół do ujawnienia neuronalnych mechanizmów działania neuroterapeutycznego na elektroencefalografii z wykorzystaniem ludzkiego neurosolvera kory nowej mózgu Wprowadzenie W ostatnich latach neuroterapia stała się coraz bardziej popularna jako metoda leczenia różnych zaburzeń psychicznych i neurologicznych. Jednak mechanizmy działania neuroterapii na poziomie neuronalnym pozostają częściowo niejasne. Elektroencefalografia (EEG) jest jednym z najczęściej stosowanych narządów diagnostycznych w neuroterapii, ponieważ umożliwia rejestrację aktywności elektrycznej mózgu w czasie rzeczywistym. W tym artykule przedstawiamy protokół do badania neuronalnych mechanizmów działania neuroterapii na EEG za pomocą ludzkiego neurosolvera kory nowej mózgu. Materiał i metody 1. Przygotowanie danych Aby rozpocząć, zebrane zostaną dane EEG od pacjentów, którzy przechodzą neuroterapię. Dane te będą przetwarzane w celu usunięcia artefaktów i zakłóceń, a następnie zostaną przefiltrowane w celu uzyskania odpowiednich zakresów częstotliwości. 2. Modelowanie matematyczne Następnie, wykorzystując ludzki neurosolver kory nowej mózgu, dane EEG zostaną przekształcone w model matematyczny opisujący aktywność neuronalną w korze nowej mózgu. Neurosolver kory nowej mózgu jest zaawansowanym narzędziem komputerowym, które pozwala na symulację działania neuronów i sieci neuronowych w korze mózgu. 3. Analiza wrażliwości Po utworzeniu modelu matematycznego, przeprowadzona zostanie analiza wrażliwości w celu zidentyfikowania kluczowych parametrów i zmiennych, które mają wpływ na działanie neuroterapii. Analiza ta pomoże zrozumieć, jak neuroterapia modyfikuje aktywność neuronalną w korze mózgu i jak te zmiany wpływają na działanie terapeutyczne. 4. Weryfikacja modelu Na koniec, model matematyczny zostanie zweryfikowany poprzez porównanie jego prognoz z rzeczywistymi danymi EEG pacjentów, którzy przeszli neuroterapię. Ta weryfikacja pomoże określić dokładność i przydatność modelu w przewidywaniu działania neuroterapii. Wyniki Protokół ten pozwoli na lepsze zrozumienie mechanizmów działania neuroterapii na poziomie neuronalnym. Dzięki zastosowaniu ludzkiego neurosolvera kory nowej mózgu, będziemy w stanie zidentyfikować kluczowe elementy aktywności neuronalnej, które są modyfikowane przez neuroterapię i które mają bezpośredni wpływ na efekty terapeutyczne. Podsumowanie Przedstawiony protokół do badania neuronalnych mechanizmów działania neuroterapii na EEG z wykorzystaniem ludzkiego neurosolvera kory nowej mózgu zapewnia nowatorskie podejście do zrozumienia i ulepszania neuroterapii. Dzięki temu protokołowi będziemy w stanie lepiej zrozumieć, jak neuroterapia działa na poziomie neuronalnym i jak te informacje mogą być wykorzystane do ulepszania metod leczenia różnych zaburzeń psychicznych i neurologicznych.

Podsumowanie

Ten protokół pokazuje, jak symulacje neuronalne oparte na fizyce mogą być wykorzystywane do interpretacji biomarkerów elektro fizjologicznych neuroterapeutyków i ujawniania ich wpływu na obwody neuronalne, zapewniając uzasadnione mechanistycznie podejście do rozwoju neuroterapeutyki.

Streszczenie

Elektroencefalografia (EEG) i metody elektrofizjologiczne dostarczają biomarkerów o rozdzielczości milisekundowej w przypadku zaburzeń ośrodkowego układu nerwowego i są szeroko stosowane do oceny efektów związanych z leczeniem. Jednak ograniczona wiedza na temat mechanizmów neuronowych generowanych przez te biomarkery utrudnia rozwój diagnostyk i terapii opartych na tych sygnałach. Human Neocortical Neurosolver (HNN) to oprogramowanie do modelowania biofizycznego o otwartym kodzie źródłowym, które łączy zlokalizowany biomarker EEG z ich wieloskalowymi generatorami neuronowymi. Ten protokół przedstawia przepływ pracy oparty na hipotezie, wykorzystujący HNN do testowania mechanizmów neuronowych wywołanych biomarkerami EEG przez optymalizację parametrów modelu, aby osiągnąć dopasowanie między symulowanymi a empirycznymi falami źródeł prądowych. Można wtedy wizualizować i ilościowo określać odpowiednią wieloskalową aktywność komórkową i obwodową, zapewniając cele walidacji dla przewidywań modelu w następczych badaniach empirycznych. Podano przykład, który pokazuje, jak badać mechanizmy neuronowe leżące u podstaw wczesnych składowych potencjału związanego z wydarzeniem w słyszalnej odpowiedzi wywołanej (P1, N1 i P2), oraz ocenić zmiany po modyfikacjach wywołanych neuroterapeutycznie w aktywności obwodów neuronowych. Ten protokół umożliwia projektowanie eksperymentów symulacyjnych w celu wygenerowania przewidywań, które można przetestować, łącząc biomarkery EEG z podstawowymi mechanizmami obwodów neuronowych. Podobny przepływ pracy można zastosować do badania mechanizmów chorobowych lub innych interwencji terapeutycznych.

Wprowadzenie

Opracowywanie leków ukierunkowanych na ośrodkowy układ nerwowy (OUN) wiąże się z wyjątkowymi wyzwaniami, a wskaźniki zatwierdzeń są niższe niż w przypadku innych obszarów terapeutycznych, co podkreśla potrzebę stosowania innowacyjnych podejść metodologicznych¹, szczególnie tych, które pozwalają na wykrycie efektów leczenia w obrębie dynamiki mózgu. Ugruntowaną metodą badania wpływu środków terapeutycznych na aktywność neuronalną jest elektroencefalografia (EEG)^2,3. EEG dostarcza charakterystycznego obrazu dynamiki mózgu na poziomie obwodów in vivo i oferuje silny potencjał translacyjny z modeli gryzoni na badania z udziałem ludzi, ponieważ obwody neuronalne generujące sygnały EEG wykazują homologię międzygatunkową^4,5,6,7,8. W rozwoju farmaceutycznym EEG może pełnić wiele funkcji, w tym służyć jako translacyjny wskaźnik między badaniami na zwierzętach a badaniami z udziałem ludzi, umożliwiać ocenę bezpieczeństwa leku, kierować wyborem związku, dostarczać informacji o zależnościach dawka-odpowiedź, pozwalać na ocenę dowodu mechanizmu (proof-of-mechanism) we wczesnych fazach klinicznych oraz umożliwiać stratyfikację prób klinicznych i wzbogacanie kohort^{9,10,11,12,13,14}. Pomimo tych zalet, interpretacja sygnałów EEG pozostaje głównym wyzwaniem, szczególnie w przypadku prób powiązania zaobserwowanych zmian z leżącymi u ich podstaw mechanizmami neuronalnymi.

Solidnym biomarkerem EEG stosowanym w poszukiwaniu leków oddziałujących na OUN jest potencjał wywołany (ERP). ERP odzwierciedlają zsynchronizowaną czasowo aktywność mózgu wywołaną bodźcami sensorycznymi i są szeroko wykorzystywane do badania zaburzeń neurorozwojowych i neuropsychiatrycznych, w tym depresji^15,16, schizofrenii^17,18, zaburzeń spektrum autyzmu^19,20 oraz choroby Alzheimera²¹. ERP są również używane do oceny efektów terapeutycznych i zakresów dawek w odniesieniu do obwodów mózgowych^22,23,24,25, gdzie normalizacja w kierunku odpowiedzi zdrowych osobników może wskazywać na skuteczność terapeutyczną²⁶. Jednak kluczowym ograniczeniem ERP i innych biomarkerów EEG (np. oscylacji mózgowych) jest fakt, że ich powiązania ze stanami chorobowymi lub efektami leków mają w dużej mierze charakter korelacyjny. Choć analizy statystyczne pozwalają zidentyfikować zależności między biomarkerami a wynikami, nie dostarczają one wglądu w mechanizmy wyjaśniające, w jaki sposób konkretne elementy obwodów neuronalnych generują te sygnały. W związku z tym przyczynienie się konkretnych typów komórek i mechanizmów obwodowych do tych zjawisk pozostaje niejasne. Zrozumienie komórkowego i obwodowego pochodzenia sygnałów EEG mogłoby znacznie zwiększyć ich wartość poprzez powiązanie obserwowanych sygnatur z leżącą u ich podstaw fizjologią^27,28. W niniejszym manuskrypcie termin „biomarker” EEG odnosi się do mierzalnych zmian w sygnałach EEG po interwencji terapeutycznej, zgodnie z definicją w ramach FDA–NIH BEST (Biomarkers, EndpointS, and other Tools)²⁹, a nie oznacza formalnej kwalifikacji do konkretnego zastosowania klinicznego³⁰.

Podczas gdy inwazyjne rejestracje elektrofizjologiczne mogą dostarczać szczegółowych informacji na poziomie komórkowym i obwodowym, podejścia te są w dużej mierze ograniczone do modeli zwierzęcych i trudne do bezpośredniego przeniesienia na badania z udziałem ludzi. Alternatywne podejścia, takie jak techniki modelowania odwrotnego, pozwalają na szacowanie aktywności źródła na podstawie sygnałów EEG, jednak często brakuje im jawnych reprezentacji mechanistycznych leżących u podstaw obwodów neuronalnych. Symulacje biofizyczne oferują komplementarne ramy poprzez modelowanie procesów fizycznych, dzięki którym obwody neuronalne generują mierzalne sygnały EEG^31,32,33,34 (Rycina 1). W porównaniu do czysto statystycznych analiz biomarkerów lub metod odwrotnych bez ugruntowania mechanistycznego, modelowanie biofizyczne umożliwia bezpośrednie testowanie hipotez łączących dynamikę obwodów neuronalnych z obserwowanymi sygnałami elektrofizjologicznymi.

Rycina 1. Modelowanie biofizyczne w celu opracowania i przetestowania hipotez mechanistycznych leżących u podstaw farmakologicznych biomarkerów elektroencefalograficznych (EEG). (A) Identyfikacja biomarkera EEG w oparciu o różnice w sygnałach mózgowych między warunkami. Przykładem jest potencjał wywołany (ERP) odpowiedzi słuchowej, który jest zredukowany w stanie po leczeniu (kolor czerwony) w stosunku do stanu przed leczeniem (kolor niebieski). (B) Modelowanie biofizyczne umożliwia testowanie hipotez mechanistycznych wyjaśniających, w jaki sposób powstają biomarkery EEG i jak zmieniają się pod wpływem interwencji farmakologicznej. Formułowane są hipotezy dotyczące zmian w aktywności neuronalnej indukowanych przez lek, a następnie identyfikuje się odpowiadające im parametry modelu. (C) Domyślny model Human Neocortical Neurosolver (HNN) służy jako punkt wyjścia do testowania hipotez poprzez ręczną modyfikację parametrów modelu lub stosowanie zautomatyzowanych algorytmów optymalizacji i wnioskowania. Różnice w wartościach parametrów między stanem przed leczeniem a po leczeniu odpowiadają przewidywaniom opartym na modelu. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Niniejszy protokół wykorzystuje Human Neocortical Neurosolver (HNN), otwartoźródłowe środowisko do modelowania biofizycznego, w celu powiązania biomarkerów ERP efektów związanych z leczeniem z ich podstawowymi mechanizmami na poziomie komórkowym i obwodowym³³ (Rycina 2). HNN opiera się na zasadzie, że synchroniczny przepływ prądu wewnątrzkomórkowego w uporządkowanych dendrytach neuronów piramidalnych generuje główne dipole prądowe leżące u podstaw sygnałów EEG^6,35,36,37. Model przedstawia kanoniczną kolumnę neokortykalną składającą się z pobudzających neuronów piramidalnych i hamujących interneuronów rozmieszczonych w warstwach korowych 2/3 i 5. Domyślna sieć HNN zawiera 100 neuronów piramidalnych i 33 neurony hamujące na warstwę, tworząc uproszczoną, lecz ugruntowaną biologicznie reprezentację obwodów korowych. Neurony piramidalne są modelowane za pomocą wielokomorowych struktur dendrytycznych, aby uchwycić kluczowe cechy morfologiczne³⁸, natomiast neurony hamujące są reprezentowane jako pojedyncze komory ze względu na ich ograniczony wkład w prądy zewnątrzkomórkowe³³. Interakcje synaptyczne obejmują pobudzające receptory kwasu α-amino-3-hydroksy-5-metylo-4-izoksazolopropionowego (AMPA) i N-metylo-D-asparaginianu (NMDA) oraz hamujące receptory kwasu gamma-aminomasłowego typu A i kwasu gamma-aminomasłowego typu B (GABA_B), przy czym wszystkie neurony zawierają aktywne przewodności jonowe sterowane dynamiką Hodgkina-Huxleya.

Rysunek 2. Schemat modelu HNN. Wizualizacja głównych komponentów modelu HNN, w tym lokalnych połączeń sieciowych między neuronami pobudzającymi a hamującymi oraz egzogennych dróg wejściowych określanych jako „napęd proksymalny” (proximal drive) i „napęd dystalny” (distal drive). Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.

Aktywność neuronalna w HNN jest napędzana przez wejścia egzogenne reprezentujące ścieżki sprzężenia w przód oraz sprzężenia zwrotne. „Proksymalne” napędy sprzężenia w przód odpowiadają wejściom z wzgórza lemniskalnego, które celują w dendryty proksymalne, podczas gdy „dystalne” napędy sprzężenia zwrotnego reprezentują wejścia kortyko-kortykalne i wzgórzowe nielemniskalne celujące w dendryty dystalne. Wejścia te są modelowane jako pociągi potencjałów czynnościowych, które wywołują prądy synaptyczne i generują wewnątrzkomórkowy przepływ prądu wzdłuż dendrytów neuronów piramidowych. Wynikający z tego dipol prądowy na poziomie populacji jest wyrażony w nanoamperometrach, co umożliwia bezpośrednie porównanie z danymi EEG lub magnetoencefalografii (MEG) z lokalizacją źródeł przy ograniczonej orientacji. Domyślna parametryzacja HNN opiera się na danych empirycznych z badań nad korą somatosensoryczną^39,40,41 i została z sukcesem zastosowana do sygnałów z kory słuchowej^42,43,44, wzrokowej⁴⁵ oraz czołowej⁴⁶, a przewidywania wynikające z modelu zostały zwalidowane w późniejszych badaniach eksperymentalnych^7,41,47.

Symulacje HNN mogą być stosowane na wielu etapach badań i rozwoju farmaceutycznego, w tym podczas walidacji celu, porównywania mechanizmów działania leków, optymalizacji dawkowania oraz generowania hipotez dla kolejnych eksperymentów^14,48,49,50. Pozwala to użytkownikom na włączenie modelowania mechanistycznego do praktycznych procesów badawczych, wspierając generowanie i testowanie hipotez dotyczących wpływu neuroterapiutyków na obwody neuronalne. W niniejszym protokole skupiamy się na wczesnych komponentach P1, N1 i P2 słuchowych ERP, ponieważ cechy te są dobrze scharakteryzowane i stanowią ograniczenia dla modelowania opartego na hipotezach⁵¹. Choć głównym punktem zainteresowania są zmiany indukowane lekami, podejście to można rozszerzyć na inne interwencje neuroterapeutyczne, takie jak stymulacja mózgu lub trening behawioralny, a także na badania nad zaburzeniami OUN.

Stosowanie HNN opiera się na iteracyjnym schemacie modelowania, w którym struktura i parametry modelu są początkowo ograniczone przez istniejące dane, a następnie doprecyzowane poprzez porównanie z obserwacjami empirycznymi. Wielkoskalowe modele neuronowe zawierają wiele parametrów, lecz jedynie podzbiór – określany jako parametry interesujące – jest dostosowywany w celu przetestowania konkretnych hipotez. Parametry te nie są wybierane w sposób dowolny; są one dobierane na podstawie wcześniejszych dowodów eksperymentalnych oraz literatury opisującej potencjalne mechanizmy działania neuroterapeutyku. W niniejszym protokole parametry związane z czasem i siłą wejścia egzogennego, lokalną łącznością hamującą oraz przewodnictwem kanałów jonowych w dendrytach zostały wybrane jako przykłady biologicznie interpretowalnych zmiennych, na które mogą wpływać neuroterapeutyki.

Zaczynając od modelu domyślnego, użytkownicy w pierwszej kolejności dopasowują parametry do danych ERP sprzed leczenia, stosując kombinację ręcznego strojenia i automatycznej optymalizacji. Ręczne strojenie pozwala na regulację globalnego skalowania oraz parametrów wejściowych w celu przybliżenia empirycznego przebiegu fali, co zapewnia intuicyjne zrozumienie wpływu zmian parametrów na wynik modelu. Następnie stosuje się metody automatyczne, takie jak strategia ewolucyjna adaptacji macierzy kowariancji (CMA-ES), optymalizacja bayesowska oraz optymalizacja z ograniczeniami metodą aproksymacji liniowej, aby doprecyzować wartości parametrów i poprawić dopasowanie. Po ustanowieniu modelu sprzed leczenia, parametry, które według hipotezy odpowiadają za zmiany po leczeniu, są dostosowywane do danych ERP uzyskanych po przeprowadzeniu leczenia.

Aby rozwiązać problem niepewności w szacowaniu parametrów, stosuje się wnioskowanie oparte na symulacjach (SBI), w celu oszacowania rozkładów wartości parametrów, które odtwarzają zaobserwowane dane^52,53. SBI uwzględnia możliwość, że wiele kombinacji parametrów może generować podobne wyniki, i umożliwia kwantyfikację niepewności parametrów. Różnice między rozkładami parametrów przed i po zastosowaniu leczenia można ocenić za pomocą indeksu nakładania się (OVL)^54,55, co pozwala na uzyskanie informacji o potencjalnych mechanizmach działania.

Kluczową zaletą tego podejścia jest to, że dopasowanie modelu do konkretnej modalności danych generuje przewidywania na wielu poziomach aktywności neuronalnej, w tym w zakresie impulsowania komórek, lokalnych potencjałów polowych (LFPs) specyficznych dla warstw oraz gęstości prądów źródłowych (CSD). Przewidywania te stanowią cele dla walidacji eksperymentalnej z wykorzystaniem technik komplementarnych. Jeśli przewidywania nie zostaną potwierdzone danymi empirycznymi, model można zaktualizować poprzez wprowadzenie nowych ograniczeń, co tworzy iteracyjny cykl generowania hipotez, ich testowania i udoskonalania (Rysunek 3).

Rysunek 3. Iteracyjny schemat opracowywania i testowania predykcji biomarkerów ERP za pomocą HNN. Schemat odpowiada poszczególnym etapom protokołu. Identyfikacja biomarkera EEG oraz inicjalizacja domyślnego modelu HNN są zaznaczone na czerwono (kroki 1–2). Ręczne dostrajanie i optymalizacja służą do dopasowania parametrów modelu do sygnałów ERP przed i po leczeniu (kolor fioletowy; kroki 3–5). Kwantyfikacja niepewności z wykorzystaniem wnioskowania opartego na symulacjach (SBI) jest zaznaczona na zielono (krok 6). Następnie predykcje modelu są analizowane i porównywane z danymi eksperymentalnymi w celu walidacji lub dalszego ograniczenia modelu (kolor pomarańczowy; krok 7). Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Niniejszy protokół został opracowany do użytku z danymi EEG lub MEG o ograniczonej orientacji i zlokalizowanym źródle, zebranymi podczas paradygmatów odpowiedzi wywołanych. Do wygenerowania wymaganych danych wejściowych można wykorzystać standardowe metody preprocessingu i lokalizacji źródła (np. minimum norm estimation [MNE]-Python⁵⁶). Sygnały na poziomie źródła wyrażone w nanoamperometrach są bezpośrednio porównywalne z wynikami HNN. W przypadku szybkich odpowiedzi sensorycznych sygnały na poziomie źródła i sensora są często bardzo podobne, co pozwala na wykorzystanie wniosków z modelowania z lokalizacją źródła do interpretacji danych EEG na poziomie sensora^57,58.

Protokół

Wszystkie procedury obejmujące dane ludzkie zostały przeprowadzone zgodnie z odpowiednimi wytycznymi i regulacjami instytucjonalnymi. Zbiór danych wykorzystany w niniejszym badaniu pochodzi z wcześniej opublikowanej pracy⁴³ i nie wymagał dodatkowej zgody komisji etycznej. W niniejszym protokole nie wykorzystano materiałów ani procedur niebezpiecznych.

1. Identyfikacja biomarkera potencjału wywołanego w EEG indukowanego terapią i zdefiniowanie hipotez modelu

1. Zebrać lub zidentyfikować zestaw danych zawierający eksperymentalnie zarejestrowane sygnały EEG od badanych osób (np. przed i po leczeniu w kontekście neuroterapii). Rejestrować pomiary EEG podczas prezentacji bodźca sensorycznego i jednocześnie zapisywać znaczniki czasu wystąpienia bodźca sensorycznego wraz z danymi EEG, aby umożliwić segmentację na próby. Upewnić się, że dane EEG są przechowywane w formacie kompatybilnym z oprogramowaniem do wstępnego przetwarzania (np. .fif, .set lub .edf).
   UWAGA: Powiązane repozytorium kodu (https://github.com/ntolley/hnn_jove) udostępnia pliki danych wykorzystane do wygenerowania reprezentatywnych wyników. Repozytorium zawiera wstępnie przetworzone ERP z MEG słuchowego według Kohl i wsp. (2022), które służy jako ERP przed leczeniem (oryginalne dane dostępne pod adresem: https://github.com/kohl-carmen/HNN-AEF). Hipotetyczny ERP po leczeniu jest generowany poprzez skalowanie przebiegu przed leczeniem przy użyciu okna z wygaszaniem gaussowskim. Odpowiadające im pliki danych znajdują się w repozytorium w lokalizacji data/pre-treatment.txt i data/post-treatment.txt. Ponieważ sygnały MEG i EEG odzwierciedlają podobne podstawowe generatory neuronalne, protokół ten ma zastosowanie do obu modalności.
2. Zidentyfikować zestaw kandydujących cech biomarkerów ERP, które w hipotezie mają odróżniać efekty związane z leczeniem (np. czas wystąpienia i amplitudy pików ERP).
   UWAGA: W tym przykładowym protokole amplitudy pików są wykorzystywane jako biomarkery zainteresowania.
3. Przeprowadzić wstępne przetwarzanie danych EEG i wyekstrahować interesujące cechy biomarkerów.
   UWAGA: Wstępne przetwarzanie i analizę ERP wspiera kilka pakietów oprogramowania, w tym MNE-Python⁵⁶, EEGLAB⁵⁹ oraz FieldTrip⁶⁰. Do modelowania sygnałów ERP zaleca się lokalizację źródeł, choć nie jest ona wymagana. Przykładowy przepływ pracy jest dostępny pod adresem https://jonescompneurolab.github.io/hnn-core/stable/auto_examples/workflows/plot_simulate_somato.html. W kilku wcześniejszych pracach opisano szczegółowo wstępne przetwarzanie i analizę sygnałów EEG; czytelników zachęca się szczególnie do zapoznania się z publikacjami^56,61 w celu uzyskania pełniejszego wprowadzenia.
   1. Wykonać lokalizację źródeł przy użyciu sygnałów na poziomie czujników ze wszystkich kanałów lub wybrać czujniki EEG do analizy. Użyć danych z lokalizacją źródeł do bezpośredniego porównania z wynikami modelu; dane na poziomie czujników nie będą miały odpowiednio dopasowanych jednostek.
      UWAGA: Jednostkowa odpowiedniość jeden-do-jednego opisana poniżej nie będzie obowiązywać dla sygnałów na poziomie czujników.
   2. Posegmentować zarejestrowane dane EEG na próby, korzystając ze znaczników czasu wystąpienia bodźca sensorycznego.
   3. Obliczyć uśrednione dla prób przebiegi ERP dla warunków przed leczeniem i po leczeniu.
   4. Wyekstrahować kandydujące biomarkery ERP z uśrednionych przebiegów (np. obliczyć amplitudy piku N1). Przed ekstrakcją zdefiniować kryteria detekcji pików (np. okno czasowe i polaryzację).
4. Przeprowadzić testy statystyczne w celu ustalenia, które cechy ERP różnią się istotnie pomiędzy warunkami (np. przed leczeniem a po leczeniu). Wybrać odpowiednie testy statystyczne w zależności od projektu badania i w razie potrzeby zastosować korektę dla wielokrotnych porównań (np. ANOVA z powtarzanymi pomiarami, a następnie test post-hoc Tukeya HSD dla wielokrotnych porównań).
   UWAGA: Przykład kodu dotyczący testowania statystycznego jest dostępny pod adresem https://mne.tools/stable/auto_tutorials/stats-sensor-space/20_erp_stats.html.
5. Wyznaczyć konkretne, statystycznie istotne, rozróżniające cechy biomarkerów EEG (np. różnice w amplitudach N1). Zapisać wyniki do wykorzystania w kolejnych krokach.
6. Sformułować w oparciu o literaturę hipotezy dotyczące mechanizmów działania leku i związanych z nimi parametrów modelu, które są przedmiotem zainteresowania. Przeanalizować wcześniejszą literaturę i dane eksperymentalne, aby zidentyfikować właściwości biofizyczne zmienione przez neuroterapię, które mogą tłumaczyć różnice w cechach.
7. Zidentyfikować, które parametry biofizycznego modelu neuronalnego (HNN) są bezpośrednio reprezentowane lub pośrednio powiązane z właściwościami biologicznymi zidentyfikowanymi w kroku 1.6. Zdefiniować je jako parametry zainteresowania. Powiązać mechanizmy biologiczne z parametrami modelu, korzystając z literatury i dokumentacji HNN.
8. Wyznaczyć zidentyfikowany zestaw parametrów modelu zainteresowania odpowiadający właściwościom biofizycznym, które w hipotezie generują zidentyfikowane różnice w cechach EEG. Użyć domyślnego modelu HNN (inicjalizowanego w kroku 2) jako punktu wyjścia dla wszystkich wartości parametrów i zapisać wyniki do kolejnych kroków.

2. Inicjalizacja domyślnego modelu HNN: Instalacja oprogramowania do modelowania i konfiguracja folderu projektu

UWAGA: Wersje oprogramowania wykorzystane w niniejszym badaniu są określone w Tabeli Materiałów, wraz z minimalnymi wymaganiami systemowymi. Dostępne są różne opcje instalacji (tzn. pip, conda oraz instalacja ze źródeł) dla systemów Linux, macOS i Windows.

1. Pobrać i zainstalować działającą wersję Anaconda Python. Utworzyć i aktywować nowe środowisko Python w celu instalacji wymaganych pakietów oprogramowania.
2. Zainstalować oprogramowanie do biofizycznego modelowania sieci neuronowych HNN-core, korzystając z instrukcji instalacji dla konkretnego systemu operacyjnego dostępnych pod adresem https://jonescompneurolab.github.io/textbook/content/01_getting_started/installation.html.
   UWAGA: Aby efektywnie zainstalować zależności oprogramowania wykorzystane w niniejszym badaniu, powiązane repozytorium kodu (https://github.com/ntolley/hnn_jove) wykorzystuje pixi (https://pixi.prefix.dev/latest/). Należy postępować zgodnie z instrukcjami w pliku README repozytorium, aby zainstalować pixi i skonfigurować lokalną wersję repozytorium kodu.
3. Zweryfikować, czy zainstalowana wersja oprogramowania do biofizycznego modelowania sieci neuronowych to 0.6.0 lub nowsza, wpisując w terminalu następującą komendę: pip show hnn_core
4. Upewnić się, że środowisko Python jest aktywne i że instalacja zakończyła się pomyślnie. Uruchomić graficzny interfejs użytkownika (GUI), wpisując w terminalu hnn-gui i naciskając Enter.
5. Utworzyć w systemie plików komputera nowy folder projektu do przechowywania wszystkich plików danych wygenerowanych w ramach niniejszego protokołu. Utworzyć folder w dostępnym katalogu (np. w katalogu domowym lub roboczym katalogu projektu).

3. Ustalenie dopasowania modelu przed obróbką za pomocą strojenia ręcznego

1. Rozpocznij od kanonicznej symulacji ERP w HNN i jej domyślnych parametrów. Ręcznie dostrój współczynnik skalowania oraz parametry wzbudzenia egzogennego, aby dopasować je do ERP przed leczeniem (np. ERP przed leczeniem).
   UWAGA: Interfejs GUI HNN automatycznie ładuje parametry modelu dopasowane do ERP somatosensorycznego⁴⁰, który w licznych badaniach okazał się dobrym punktem wyjścia jako „kanoniczny ERP”. Niniejszy samouczek koncentruje się na modyfikacji współczynnika skalowania i parametrów wejścia egzogennego w odniesieniu do tego punktu wyjścia.
2. Załaduj empiryczny przebieg fali ERP przed leczeniem z Kroku 1 do interfejsu GUI HNN (Rysunek 4A–4F)
   1. Kliknij przycisk Load data na pasku menu znajdującym się w dolnej lewej części okna GUI (Rysunek 4D).
      UWAGA: Nomenklatura dotycząca nazewnictwa pików ERP znacznie różni się w literaturze; etykiety P1/N1/P2 na Rysunku 4F służą wyłącznie celom ilustracyjnym i mogą nie odpowiadać konwencjom nazewniczym stosowanym w innych badaniach.
   2. W oknie przeglądarki plików wybierz plik .csv lub .txt zawierający przebieg fali ERP do zamodelowania (tzn. przebieg docelowy). Upewnij się, że plik jest rozdzielony przecinkami i sformatowany w dwóch kolumnach: pierwsza kolumna zawiera czas (ms), a druga kolumna zawiera zlokalizowany źródłowo empiryczny przebieg dipola (nAm). Pierwszy wiersz jest traktowany jako nagłówek i nie powinien zawierać wartości danych. Opcjonalnie można dodać informacyjne etykiety kolumn (np. „Time (ms)” i „Dipole (nAm)”).}
      UWAGA: Empiryczny plik danych w powiązanym repozytorium kodu nazywa się pre-treatment.txt.
   3. Sprawdź przebieg fali, który zostanie automatycznie wykreślony w panelu rysunków (Rysunek 4F).
3. Uruchom domyślną symulację kanonicznego ERP
   1. Ustaw wartości parametrów tstop, dt, Trials, Backend oraz Cores w panelu Simulation Parameters (Rysunek 4B) na pożądane wartości. Użyj tstop do kontroli długości symulacji, dt do kontroli kroku czasowego całkowania oraz Trials do kontroli liczby powtórzeń symulacji uruchamianych z tymi samymi wartościami parametrów modelu. Wybierz Backend jako sekwencyjny (Joblib) lub równoległy (MPI) i określ liczbę rdzeni komputera Cores.
      UWAGA: Zmienność między próbami wynika z odchylenia standardowego czasu wzbudzenia egzogennego opisanego w Kroku 3.5 poniżej.
   2. Kliknij przycisk Run (Rysunek 4D), aby rozpocząć domyślną symulację kanonicznego ERP.
4. Utwórz wykres porównujący symulowany ERP z empirycznym ERP
   1. Kliknij kartę Visualization w górnej lewej części okna GUI (Rysunek 4A).
   2. Kliknij menu rozwijane oznaczone jako Data to compare (niewidoczne) i wybierz załadowany przebieg docelowy z Kroku 3.2.
   3. Kliknij Clear axis, aby zresetować wykres.
   4. Kliknij Add plot, aby wygenerować nowy wykres z nałożonymi symulowanym początkowym przebiegiem fali ERP (niebieski) i przebiegiem docelowym (pomarańczowy), wraz z tekstem wskazującym automatycznie obliczony współczynnik korelacji (Corr) i błąd średnio kwadratowy (RMSE) między dwoma przebiegami (Rysunek 4F).
      UWAGA: HNN-GUI oferuje możliwość obliczenia dwóch miar dopasowania: Corr i RMSE. Miary te są wykorzystywane do ręcznego strojenia i optymalizacji (Krok 4).
5. Zmodyfikuj współczynnik skalowania
   1. Zmodyfikuj współczynnik skalowania poprzez ręczne strojenie, aby w przybliżeniu dopasować amplitudy symulowanych i empirycznych przebiegów dipolowych. Ustaw domyślny parametr Dipole scaling (Rysunek 4C) w karcie Simulation (Rysunek 4A) na 3000.
      UWAGA: Współczynnik skalowania odpowiada przewidywanej szacowanej liczbie neuronów leżących u podstaw generowania sygnału EEG. Domyślna wartość 3000 wskazuje, że do wygenerowania odpowiedzi wywołanej o amplitudzie w nAm wskazanej na osi y Rysunku 4F konieczne jest 200 neuronów piramidowych (rozmiar modelu HNN) × 3000 = 600 000 neuronów.
6. Zmodyfikuj czas wzbudzeń egzogennych
   1. Zmodyfikuj średnią i odchylenie standardowe wzbudzeń egzogennych poprzez ręczne strojenie, aby uzyskać dokładniejsze dopasowanie do czasu empirycznie zarejestrowanych pików ERP przed leczeniem (tzn. P1/N1/P2) (Rysunek 5A–5D).
      UWAGA: Domyślne parametry lokalnej łączności i komórek dostarczone z HNN zostały dostrojone w celu odtworzenia zdrowych wzorców aktywności pojedynczych komórek i sieci. Chociaż parametry lokalnej sieci można dostosować, zaleca się początkowo pozostawienie stałych parametrów wstępnie dostrojonego lokalnego modelu szablonu neokortykalnego HNN i sprawdzenie, czy wiarygodne dopasowanie można osiągnąć, regulując jedynie wzbudzenia egzogenne.
   2. Zidentyfikuj, które symulowane piki ERP są przesunięte w czasie względem empirycznego przebiegu fali ERP (Rysunek 4).
      UWAGA: Ten przykład zakłada trzy wczesne piki w empirycznym ERP, podobnie jak w domyślnej symulacji kanonicznego ERP. Aby dodać piki, symuluj dodatkowe wzbudzenia zewnętrzne.
   3. Kliknij kartę External drives w górnej lewej części okna GUI (Rysunek 4A oraz Rysunek 5A).
      UWAGA: Widoczne są parametry dla trzech zdefiniowanych wzbudzeń egzogennych, reprezentujących proksymalne sprzężenie w przód (evprox1), dystalne sprzężenie zwrotne (evdist1) oraz ponownie pojawiające się proksymalne sprzężenie w przód (evprox2), które generują domyślne symulacje kanonicznego ERP (szczegóły modelu HNN i struktury wzbudzeń egzogennych znajdują się we Wstępie). Histogramy przedstawiające czasy i liczbę impulsów pokazano na Rysunku 4E.
   4. Kliknij menu rozwijane wzbudzenia egzogennego, którego Mean time jest najbliższe przesuniętemu piku.
   5. Zmodyfikuj wartości w polach tekstowych dla Mean time i Std dev time, aby lepiej dopasować czas i szerokość pików w przebiegu docelowym (Rysunek 5B–5D). Dostosuj Mean time, aby przesunąć czas piku, oraz Std dev time, aby zmienić szerokość piku.
      UWAGA: Mean time i Std dev time kontrolują średnią i wariancję impulsów egzogennych, które aktywują lokalną sieć w proksymalnych lub dystalnych wzorcach projekcji (patrz histogramy na Rysunku 4E). Parametry te nie determinują w pełni czasu ani szerokości piku ERP. Dokładny czas i szerokość zależą zarówno od wzbudzeń egzogennych, jak i od wewnętrznej aktywności sieci.
      1. Ustaw Mean time dla wzbudzenia zewnętrznego evprox1 na 60 ms.
      2. Ustaw Mean time dla wzbudzenia zewnętrznego evdist1 na 100 ms.
      3. Ustaw Mean time dla wzbudzenia zewnętrznego evprox2 na 150 ms.
7. Zmodyfikuj amplitudę wzbudzeń egzogennych
   1. Zmodyfikuj wagi synaptyczne (przewodność postsynaptyczną) wzbudzeń egzogennych poprzez ręczne strojenie, aby uzyskać dokładniejsze dopasowanie do amplitudy empirycznie zarejestrowanych pików ERP (tzn. P1/N1/P2) (Rysunek 6A oraz Rysunek 6B).
   2. Zidentyfikuj, które symulowane piki ERP są niedopasowane pod względem amplitudy do empirycznego przebiegu fali ERP.
   3. Kliknij kartę External drives w górnej lewej części okna GUI (Rysunek 4A).
   4. Kliknij menu rozwijane wzbudzenia egzogennego, którego Mean time jest najbliższe niedopasowanemu piku.
   5. Zmodyfikuj wartości w polach tekstowych pod AMPA weights i NMDA weights, aby dostosować przewodnictwa synaptyczne. Zwiększenie siły wzbudzenia proksymalnego neuronów piramidowych L5 i L2/3 zazwyczaj generuje bardziej dodatnie piki, natomiast zwiększenie siły wzbudzenia dystalnego zazwyczaj generuje bardziej ujemne piki.
      UWAGA: Podobnie jak w przypadku czasu wzbudzeń egzogennych, amplituda piku ERP nie jest w pełni determinowana przez siłę wzbudzenia. Dynamika impulsowania może wywoływać nieintuicyjne efekty. Testuj zmiany w zakresie jednego rzędu wielkości (np. AMPA L5_pyramidal od 0.014 do 0.14) i doprecyzuj je iteracyjnie. Rysunek 6 pokazuje wartości ustawione na 10× mniejsze niż w domyślnej symulacji.
      1. Ustaw AMPA weights wzbudzenia evdist1 na L5_pyramidal = 0.014243 i L2_pyramidal = 0.0000007.
      2. Ustaw NMDA weights wzbudzenia evdist1 na L5_pyramidal = 0.0080074 i L2_pyramidal = 0.0004317.
      3. Ustaw AMPA weights wzbudzenia evprox2 na L5_pyramidal = 0.0684013 i L2_pyramidal = 0.143884.
         UWAGA: Pełny zestaw parametrów użytych do wygenerowania reprezentatywnych wyników jest dostępny w powiązanym repozytorium kodu (https://github.com/ntolley/hnn_jove; patrz data/opt_baseline_config_correlation_best.json). Zachęcamy użytkowników do załadowania tego pliku konfiguracyjnego wraz z dostarczonymi plikami danych (data/pre-treatment.txt i data/post-treatment.txt) i zapoznania się z przykładowymi przepływami pracy w katalogu notebooks/, aby odtworzyć raportowane symulacje.
8. Zapisz zmodyfikowane ustawienia symulacji.
   1. Po zakończeniu modyfikacji w Krokach 3.5–3.7, kliknij kartę Simulation (Rysunek 4A) i wpisz „pre-treatment_handtuned” w polu tekstowym Name (Rysunek 4B).
9. Uruchom zmodyfikowaną symulację
   1. Kliknij przycisk Run, aby zasymulować zmodyfikowany zestaw parametrów.
   2. Sprawdź wygenerowany wykres w panelu rysunków (Rysunek 4F oraz Rysunek 7A–7D). Uzyskaj dostęp do poprzednich wykresów za pomocą odpowiednich kart rysunków (np. „Figure 1” i „Figure 2”).
10. Iteruj ręczne strojenie
    1. Kontynuuj iteracyjne ręczne strojenie, aby poprawić współczynnik korelacji.
    2. Powtórz Krok 3.4, aby ponownie wykreślić symulację z przebiegiem docelowym i ponownie obliczyć współczynnik korelacji.
11. Zapisz końcowe wyniki modelu.
    UWAGA: Protokół można przerwać po zapisaniu wyników symulacji. Wznów go, ładując zapisane pliki konfiguracyjne do oprogramowania.
    1. Kliknij przycisk Save Network, aby zapisać najlepiej dopasowany zestaw parametrów jako plik .json o nazwie „pre-treatment_handtuned.json”.
    2. Kliknij przycisk Save simulation, aby zapisać plik .txt o nazwie „pre-treatment_handtuned.txt”, który zawiera symulowany przebieg dipola (Rysunek 4D).
    3. Przenieś oba pliki do folderu projektu utworzonego w Kroku 2.5. Upewnij się, że nazwy plików odpowiadają nazwie symulacji w menu rozwijanym.
       UWAGA: Pliki są zapisywane w domyślnym katalogu pobierania przeglądarki internetowej użytej do uruchomienia GUI. Przenieś pliki ręcznie lub tymczasowo zmień katalog pobierania w przeglądarce.

Rycina 4. Porównanie symulowanego kształtu fali ERP kanonicznej sieci HNN z empirycznym ERP przed leczeniem. (A) Kategorie parametrów dostępnych za pomocą kart graficznego interfejsu użytkownika (GUI). (B) Parametry symulacji kontrolujące długość symulacji i liczbę prób. (C) Parametry wizualizacji kontrolujące wyświetlanie kształtu fali. (D) Panel sterowania symulacją służący do ładowania danych, uruchamiania symulacji i zapisywania wyników. (E) Histogramy impulsów pokazujące rozkłady egzogennych wejść sterujących w kanonicznej symulacji ERP. (F) Kształt fali dipola kanonicznej symulacji ERP (niebieski) nałożony na empiryczny słuchowy ERP (pomarańczowy) z pracy Kohl et al.⁴³. Początkowa symulacja nie pasuje do danych, z przesuniętym czasem wystąpienia szczytu i jego amplitudą (Corr < 0.95). Paradygmat eksperymentalny użyty do wygenerowania empirycznego ERP opisano w Kohl et al.⁴³: tony (1 kHz, czas trwania 50 ms, 10 ms wyciszenia/narastania) były prezentowane naprzemiennie do lewego i prawego ucha, z odstępami międzybodźcowymi 0.8–1.2 s przy poziomie 60 dB powyżej subiektywnego progu słyszalności. Kliknij tutaj, aby wyświetlić większą wersję tej ryciny.

Rycina 5. Modyfikacja czasu zewnętrznych sygnałów sterujących w celu dopasowania szczytów ERP. (A) Karta „External drives” w interfejsie GUI, służąca do konfiguracji wejść wywoływanych w modelu. (B–D) Dostosowanie parametrów średniego czasu dla poszczególnych zewnętrznych sygnałów sterujących w celu dopasowania symulowanych szczytów ERP do danych empirycznych. W szczególności: (B) proksymalny sygnał sterujący evprox1 dopasowany do ~60 ms, (C) dystalny sygnał sterujący evdist1 dopasowany do ~100 ms oraz (D) proksymalny sygnał sterujący evprox2 dopasowany do ~150 ms. Korekta parametru średniego czasu (wyróżniona) przesuwa moment wystąpienia symulowanych szczytów i poprawia zgodność z empirycznym przebiegiem fali. Dostosowania te przyczyniają się do lepszego dopasowania i zwiększenia korelacji z docelowym ERP (patrz Rycina 7B). Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rycina 6. Modyfikacja siły pobudzenia egzogennego w celu dostosowania amplitud szczytów ERP. (A i B) Wagi synaptyczne receptorów kwasu α-amino-3-hydroksy-5-metylu-4-izoksazolopropionowego (AMPA) i N-metylu-D-asparaginianu (NMDA) są modyfikowane za pomocą karty „External drives” w graficznym interfejsie użytkownika (GUI). (A) Dostosowanie wag synaptycznych dla pobudzenia dystalnego (evdist1), w tym przewodności AMPA i NMDA skierowanych do neuronów piramidalnych warstwy 2/3 (L2/3) i warstwy 5 (L5). (B) Dostosowanie wag synaptycznych dla pobudzenia proksymalnego (evprox2), wpływające głównie na przewodności AMPA w neuronach piramidalnych. W tym przykładzie wagi synaptyczne zostały zmniejszone dziesięciokrotnie w stosunku do wartości domyślnych, co doprowadziło do zmniejszenia amplitud szczytów ERP i lepszej zgodności z empirycznym przebiegiem fali (patrz Rycina 7C). Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

Rycina 7. Ręczne dostrajanie i optymalizacja parametrów modelu. Wszystkie symulacje przedstawiają 5 prób, z uśrednionym ERP (ciemnoniebieski) oraz poszczególnymi próbami (jasnoniebieski). (A) Kanoniczna symulacja ERP (niebieska) nałożona na ERP sprzed leczenia (pomarańczowa). (B) Korekta czasu wymuszenia egzogennego poprawia dopasowanie szczytów. (C) Zmniejszenie wag synaptycznych obniża amplitudę szczytów. (D) Automatyczna optymalizacja daje bliskie dopasowanie do empirycznej formy fali (Corr = 1.0), w tym zwiększoną zmienność czasu wymuszenia wywołanego. Kliknij tutaj, aby wyświetlić powiększoną wersję tej ryciny.

4. Ustalenie dopasowania modelu przed opracowaniem z optymalizacją parametrów

UWAGA: Kontrola losowego ziarna (seeding) w celu optymalizacji nie jest obecnie dostępna w interfejsie graficznym (GUI). Aby uzyskać powtarzalne przebiegi optymalizacji, należy użyć interfejsu Python API. Powiązane repozytorium kodu zawiera przykład implementacji (patrz code/baseline_optimization.py), w którym można ustawić stałe ziarno losowe, przekazując parametr seed do funkcji optymalizacji (np. optim.fit(..., seed=123)).

UWAGA: Ten przykład pokazuje, jak optymalizować docelowe parametry w celu oszacowania pojedynczych wartości, które zapewniają dokładne dopasowanie do przebiegu fali przy użyciu CMA-ES (nie należy mylić z SBI; oba są podejściami do dopasowania parametrów modelu, ale głównym wynikiem SBI jest rozkład). Przykład oszacowania rozkładów parametrów, które mogą odpowiadać przebiegom fal, przedstawiono w sekcji Wyniki. W przypadku ERP przed leczeniem należy rozpocząć od optymalizacji parametrów sterowania egzogennego przy założeniu, że parametry komórkowe i połączeń w lokalnej sieci w domyślnym modelu neokortykalnym HNN są stałe. Wieloskalowe przewidywania dostarczane przez HNN opisane w kroku 7 stanowią cele do walidacji tego założenia. W miarę pojawiania się nowych informacji ograniczających przewidywania modelu, struktura HNN umożliwia oszacowanie dowolnego zestawu parametrów.

1. Otwórz ustawienia optymalizacji
   1. Kliknij zakładkę Optimization w lewym górnym rogu interfejsu graficznego (Rycina 8A).
   2. Skonfiguruj ustawienia procesu optymalizacji, w tym liczbę iteracji, solver oraz funkcję celu.
      UWAGA: Domyślne ustawienia optymalizacji (Funkcja celu = „dipole_corr”; Solver = „cma”) są odpowiednie dla przebiegów ERP. Ta funkcja celu maksymalizuje współczynnik korelacji między symulowanymi a empirycznymi przebiegami. Zwiększ maksymalną liczbę iteracji w przypadku optymalizacji wielu parametrów. Współczynnik korelacji jest miarą bezwymiarową; dlatego przy użyciu „dipole_corr” należy dostosować czynnik skalowania po optymalizacji (Krok 4.7.1). Alternatywnie można użyć „dipole_rmse”, aby zminimalizować RMSE, w którym przypadku czynnik skalowania pozostaje stały.
   3. Kliknij pole tekstowe Max iterations i wpisz 100.
2. Wybierz parametry do optymalizacji
   1. Kliknij menu rozwijane dla wzbudzenia egzogennego, którego parametry mają zostać zoptymalizowane (Rycina 8A oraz Rycina 8B, czerwone koło).
   2. Wybierz parametry wzbudzenia do optymalizacji, zaznaczając pole wyboru w kolumnie „Optimized against?” (Rycina 8B).
3. Zdefiniuj ograniczenia parametrów
   1. Określ zakres wartości parametrów przeszukiwanych przez optymalizator, wpisując wartości w polach tekstowych Min i Max w sekcji Constraints (%) (Rycina 8B).
      UWAGA: Domyślne wartości 20% są odpowiednie dla symulacji, które wykazują już wysoki współczynnik korelacji (Corr > 0,9). Na przykład zastosowanie zakresu 20% do średniego czasu (Mean time) wynoszącego 65,53 ms daje granice 52,42–78,64 ms. W przypadku słabego dopasowania początkowego należy zwiększyć wartości procentowe Min i Max; może to jednak znacząco zwiększyć liczbę wymaganych symulacji.
4. Uruchom optymalizację
   1. Kliknij przycisk Run Optimization (Rycina 8A), aby wykonać procedurę optymalizacji.
5. Zapisz wyniki optymalizacji
   1. Kliknij przycisk Save Optimization History (Rycina 8A).
   2. Przenieś zapisany plik do folderu projektu utworzonego w Kroku 2.5.
      UWAGA: Wyniki optymalizacji można przechowywać i wykorzystywać ponownie. Protokół można w tym momencie przerwać i wznowić poprzez wczytanie zapisanej historii optymalizacji.
6. Oceń jakość optymalizacji
   1. Oceń jakość przeprowadzonej optymalizacji.
      UWAGA: Przy stosowaniu współczynnika korelacji jako miary dopasowania zaleca się kryterium zatrzymania Corr > 0,95, ponieważ zazwyczaj odzwierciedla to symulowany przebieg, który odtwarza główne piki i dołki docelowego ERP. Wczesne zatrzymanie nie jest obecnie wspierane, ale jest w fazie opracowywania. Zwiększ liczbę iteracji, jeśli kryterium zatrzymania nie zostało spełnione, ale strata nadal maleje co 10 iteracji.
7. Określ kolejne kroki na podstawie wyniku optymalizacji
   1. Jeśli uzyskano dobre dopasowanie do ERP przed leczeniem (tzn. Corr > 0,95), ponownie dostosuj czynnik skalowania, aby zminimalizować RMSE, i przejdź do Kroku 5.
      UWAGA: Jak opisano w Kroku 4.1, gdy jako funkcja celu stosowana jest opcja „dipole_corr”, należy ponownie dostosować czynnik skalowania po optymalizacji. W tym przykładzie czynnik skalowania został zmniejszony z domyślnej wartości 3000× (Rycina 7A–7C) do 1000× (Rycina 7D).
   2. Jeśli optymalizacja nie pozwoliła na uzyskanie dobrego dopasowania do ERP przed leczeniem, wróć do Kroku 4.2 i przeprowadź rozwiązywanie problemów poprzez zwiększenie maksymalnej liczby iteracji, poprawę ręcznie dostrojonego punktu startowego lub wybór alternatywnych parametrów do dostosowania.
      UWAGA: Szczegółowe wyjaśnienie kroków rozwiązywania problemów znajduje się w sekcji „Troubleshooting when fitting parameters to data features” w części Dyskusja.

Rysunek 8. Optymalizacja parametrów wymuszenia zewnętrznego w celu poprawy dopasowania do ERP przed leczeniem. (A) Zakładka optymalizacji w interfejsie graficznym (GUI) służąca do konfiguracji parametrów optymalizacji. (B) Wybór parametrów i zakresów ograniczeń dla optymalizacji. (C) Przykład wyniku optymalizacji wykazujący lepsze dopasowanie do empirycznych danych ERP z pracy Kohl i wsp.⁴³. (D) Krzywa straty optymalizacji wykazująca zbieżność po około 80 iteracjach. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

5. Ustalenie dopasowania modelu po przeprowadzeniu zabiegu

1. Zacznij od zoptymalizowanej symulacji ERP przed leczeniem. Ręcznie dostrój i zoptymalizuj interesujące parametry, aby pasowały do ERP po leczeniu.
2. Wczytaj empiryczną postać fali ERP po leczeniu
   1. Wczytaj empiryczną postać fali ERP po leczeniu z Kroku 1 do interfejsu GUI (ta sama procedura co w Kroku 3.2; Rysunek 9A).
3. Wczytaj zoptymalizowane parametry sprzed leczenia
   1. Wczytaj zoptymalizowane parametry ERP sprzed leczenia z Kroków 1–4 jako punkt wyjścia (Rysunek 9A).
4. Przeprowadź ręczne dostrajanie i optymalizację
   1. Przeprowadź ręczne dostrajanie i optymalizację parametrów (te same procedury co w Krokach 3.2–3.11 oraz Kroku 4) dla parametrów interesujących zidentyfikowanych w Kroku 1.7.
   2. Kontynuuj dostrajanie i optymalizację, aż do uzyskania wysokiej korelacji (Corr > 0.95) między ERP symulowanym a ERP po leczeniu.
      UWAGA: W celach ilustracyjnych na Rysunku 9B zastosowano ręczne dostrajanie parametru ukierunkowanego na sygnał (zmniejszono maksymalną przewodność GABA_B w sieci lokalnej), co pozwoliło na ściślejsze dopasowanie do danych po leczeniu. Nie przeprowadzono optymalizacji w celu oceny, w jakim stopniu ta zmiana parametru wyjaśnia dane. Sekcja „Representative Results” opisuje, jak za pomocą SBI szacować rozkłady wielu parametrów hipotezowanych jako parametry interesujące po leczeniu. SBI (szczegółowo omówione w Kroku 6) jest zalecane do rygorystycznych badań, ponieważ szacuje rozkłady parametrów wyjaśniających postać fali ERP, co umożliwia solidne porównania dopasowań parametrów.
5. Zapisz konfigurację modelu i porównaj parametry
   1. Zapisz konfigurację modelu i porównaj zoptymalizowane wartości parametrów interesujących pomiędzy warunkami przed leczeniem a po leczeniu (dane nie przedstawione).
   2. Powtórz Krok 3.11, aby wyeksportować plik .json z parametrami modelu. Przenieś plik do folderu projektu utworzonego w Kroku 2.5.
   3. Przejrzyj parametry wymuszenia egzogennego, klikając Load external drives (Rysunek 5A) i wybierając plik konfiguracji sieci sprzed leczenia lub po leczeniu.
   4. Przejrzyj parametry sieci lokalnej, klikając Load local network connectivity (Rysunek 9C) i wybierając plik konfiguracji sieci sprzed leczenia lub po leczeniu.
   5. Zidentyfikuj zmiany w wartościach parametrów pomiędzy konfiguracjami sieci sprzed leczenia a po leczeniu. Zinterpretuj te zmiany jako przewidywania modelu dotyczące mechanizmów biomarkerów po leczeniu.

Rysunek 9. Ocena siły synaptycznej kwasu gamma-aminomasłowego typu B (GABA_B) jako mechanizmu biomarkerów EEG po leczeniu. (A) Zoptymalizowana symulacja przed leczeniem (niebieska) nałożona na ERP po leczeniu (czerwona), wykazująca zmniejszone amplitudy szczytowe. (B) Redukcja siły synaptycznej GABA_B zmniejsza amplitudę N1, co sugeruje potencjalny mechanizm. (C) Panel GUI pokazujący miejsce modyfikacji lokalnej siły synaptycznej GABA_B. Kliknij tutaj, aby wyświetlić większą wersję tego rysunku.

6. Wykonanie kwantyfikacji niepewności za pomocą SBI i ocena separowalności przy użyciu interfejsu programistycznego aplikacji HNN-Python

UWAGA: SBI wymaga instalacji oddzielnego pakietu Python⁶². Należy zapoznać się z powiązanym repozytorium (https://github.com/ntolley/hnn_jove), w którym znajduje się przykład kodu szczegółowo opisujący sposób przeprowadzania wnioskowania parametrów w HNN przy użyciu pakietu oprogramowania SBI. Kod jest zorganizowany tak, aby odzwierciedlać kroki opisane w poniższym protokole. Pełna dyskusja na temat zastosowania SBI do modelu HNN została przedstawiona w⁵⁵.

1. Instalacja pakietu SBI
   1. Zainstaluj pakiet SBI, uruchamiając następującą komendę w terminalu z aktywnym środowiskiem Python: pip install sbi.
2. Zdefiniowanie przedziałów parametrów a priori
   1. Zidentyfikuj przedziały parametrów dla wybranego podzbioru parametrów ERP przed i po leczeniu, aby utworzyć ograniczony rozkład a priori do kwantyfikacji niepewności.
3. Generowanie zbioru danych treningowych.
   1. Zdefiniuj funkcję aktualizacji parametrów (stosując to samo podejście, co przy optymalizacji parametrów).
   2. Ustaw ziarno losowości (random seed) dla generowania próbek z rozkładu a priori, aby zapewnić powtarzalność. W przypadku używania NumPy do generowania losowych próbek, utwórz instancję generatora liczb pseudolosowych w skrypcie Python (np. rng = np.random.default_rng(123)) i używaj tego generatora do próbkowania.
      UWAGA: Powiązane repozytorium kodu (https://github.com/ntolley/hnn_jove) zawiera przykład użycia generatora losowego NumPy w pliku code/generate_simulations.py.
   3. Pobierz próbki parametrów z rozkładu a priori.
      UWAGA: Do wygenerowania reprezentatywnych wyników użyto 10 000 próbek.
   4. Wygeneruj zbiór danych symulowanych ERP, korzystając z pobranych wartości parametrów.
4. Wybór statystyk podsumowujących.
   1. Wybierz statystykę podsumowującą, która charakteryzuje przebieg fali EEG.
      UWAGA: Statystyka podsumowująca to dowolna wielkość, która oddaje kluczowe cechy przebiegu fali EEG. Częstym wyborem są czas wystąpienia oraz amplituda piku. W niniejszym opracowaniu do wyodrębnienia statystyk podsumowujących wykorzystano analizę głównych składowych (PCA) (tzn. ładunki czterech pierwszych składowych głównych). Pełną dyskusję na ten temat zaprezentowano w punkcie⁵⁵.
   2. Trenowanie sieci SBI
      UWAGA: W niniejszym samouczku dla obiektu estymatora a posteriori (neural posterior estimator) zastosowano domyślne parametry trenowania (np. density_estimator=”maf”, training_batch_size=200, learning_rate=0.0005) dostarczone z pakietem SBI. Parametry trenowania są opisane w dokumentacji SBI (https://sbi.readthedocs.io/en/stable/api_reference/_autosummary/sbi.inference.NPE_B.html).
   3. Ustaw globalne ziarno losowości PyTorch, aby zapewnić powtarzalność trenowania, dodając torch.manual_seed(0) w skrypcie Python po imporcie torch.
   4. Wytrenuj sieć SBI, aby mapowała kombinacje parametrów na symulowane przebiegi fal ERP.
      UWAGA: Wytrenowana sieć SBI jest obiektem Pythona, który przyjmuje jako wejście statystyki podsumowujące z danych EEG i zwraca rozkład parametrów (rozkład a posteriori). Jeśli trenowanie przebiegło pomyślnie, symulowanie parametrów z tego rozkładu w modelu HNN generuje przebiegi EEG podobne do danych empirycznych (posterior predictive check [PPC]).
   5. Generowanie próbek a posteriori i ocena dopasowania
   6. Podaj eksperymentalny przebieg fali EEG jako wejście warunkujące do wytrenowanej sieci.
   7. Pobierz próbki parametrów z rozkładu a posteriori uwarunkowanego eksperymentalnym przebiegiem fali EEG.
   8. Przeprowadź symulację dla próbek parametrów pobranych z rozkładu a posteriori.
   9. Oblicz podobieństwo między przebiegami symulowanymi a eksperymentalnym przebiegiem fali EEG podanym jako wejście.
      UWAGA: Procedura ta jest nazywana PPC. Dobrze wytrenowana sieć generuje symulacje, które ściśle odpowiadają przebiegowi empirycznemu (wysoka korelacja lub niskie RMSE). Jeśli PPC nie daje satysfakcjonujących symulacji, możliwe są dwie sytuacje: (1) hipotetyczne mechanizmy nie wyjaśniają biomarkera, co wymaga postawienia nowych hipotez i aktualizacji rozkładów a priori; lub (2) sieć SBI nie została pomyślnie wytrenowana. W takim przypadku należy zwiększyć budżet treningowy lub zmodyfikować statystyki podsumowujące.
   10. Jeśli symulowane ERP z pobranych rozkładów parametrów pasują do ERP przed i po leczeniu (PPC z Corr > 0,95), przejdź do kroku 6.8. W przeciwnym razie przejdź do kroku 6.7.
5. Rozwiązywanie problemów z trenowaniem sieci SBI
   UWAGA: Niepowodzenie PPC wskazuje na konieczność modyfikacji parametrów trenowania sieci SBI. Szczegółowe wyjaśnienia znajdują się w sekcji „Troubleshooting when fitting parameters to data features” w rozdziale Dyskusja.
   1. Zwiększ rozmiar zbioru danych treningowych.
   2. Zmodyfikuj cechy podsumowujące.
   3. Wybierz inną architekturę SBI do trenowania.
6. Wizualizacja rozkładów a posteriori i ocena rozdzielności
   1. Przekaż tablicę próbek parametrów z kroku 6.6.2 do funkcji pairplot i przypisz odrębne kolory do rozkładów odpowiadających każdemu warunkowi ERP.
      UWAGA: Powiązane repozytorium kodu zawiera funkcjonalność wykresowania pozwalającą odtworzyć Rysunek 10.
   2. Przeanalizuj panele diagonalne wygenerowanego wykresu pairplot pod kątem rozkładów niezachodzących na siebie. Oceń rozdzielność, obliczając OVL (Rysunek 10A). Parametry o silnie rozdzielnych rozkładach (OVL < 0,1) odpowiadają przewidywanym mechanizmom działania neuroterapii, które zmieniają się po leczeniu w stosunku do stanu przed leczeniem.
      UWAGA: OVL jest miarą kwantyfikującą rozdzielność rozkładów w zakresie (0,1), gdzie OVL = 0,0 oznacza brak nakładania się, a OVL = 1,0 oznacza całkowite nakładanie się^54,55. Kod do obliczania OVL znajduje się w powiązanym repozytorium kodu.

Rysunek 10. SBI do kwantyfikacji niepewności parametrów i identyfikacji mechanizmów neuroterapeutycznych. (A) Wizualizacja pairplot rozkładów parametrów oszacowanych przy użyciu SBI. Panele diagonalne (i–iv) przedstawiają rozkłady jednowymiarowe dla poszczególnych parametrów, w tym (i) synchronii wzgórzowo-korowej, (ii) przewodnictwa K_m w dendrytach, (iii) przewodnictwa GABA_B oraz (iv) siły sprzężenia zwrotnego kora-kora. Jednostki dla (i) wyrażone są jako mnożnik skalujący domyślną wartość parametru (przed leczeniem). Jednostki dla (ii-iv) wyrażone są jako mnożnik skalujący domyślną wartość parametru (przed leczeniem) w skali logarytmicznej. Rozkłady dla warunków przed leczeniem (niebieski) i po leczeniu (czerwony) wykazują różny stopień rozdzielności, przy czym synchronia wzgórzowo-korowa wykazuje najmniejsze nakładanie się (wartość nakładania, OVL = 0.07), co wskazuje na najsilniejszy efekt związany z leczeniem. Panele poza przekątną przedstawiają relacje dwuwymiarowe między parametrami. (B) Posterior predictive check (PPC) dla ERP przed leczeniem; symulowane przebiegi falowe (czarny) ściśle odpowiadają danym empirycznym (niebieski). (C) PPC dla ERP po leczeniu; symulowane przebiegi falowe (czarny) ściśle odpowiadają danym empirycznym (czerwony). Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

7. Przeprowadzenie badania, walidacji i dalszego ograniczania modelu

UWAGA: Ten krok przedstawia przykłady wizualizacji elementów symulowanej aktywności w interfejsie graficznym (GUI). Te wielkoskalowe szczegóły stanowią punkty odniesienia do walidacji i uzupełnienia przewidywań wynikających z modelu w dalszych eksperymentach^7,47. Niniejszy protokół nie zawiera wytycznych dotyczących wyboru przewidywań najlepiej nadających się do eksperymentów walidacyjnych ani sposobu przeprowadzenia takich eksperymentów (tj. Krok 7.3).

1. Wczytanie parametrów modelu i uruchomienie symulacji
   1. Wczytać parametry modelu zoptymalizowane dla warunków przed i po zabiegu, a następnie uruchomić symulacje.
      UWAGA: Można wczytać i przeanalizować parametry uzyskane z optymalizacji w krokach 4–5. Przykłady eksportu parametrów sieci wygenerowanych przez SBI w kroku 6 z interfejsu Python znajdują się w powiązanym repozytorium GitHub.
2. Analiza predykcji wieloskalowych
   1. Przeanalizować predykcje wieloskalowe na podstawie wyników symulacji.
   2. Wykres aktywności impulsowej na poziomie komórkowym
   3. Kliknąć kartę Visualization (Rysunek 4A).
   4. Kliknąć menu rozwijane Layout template i wybrać Dipole Layers-Spikes.
   5. W menu rozwijanym Dataset wybrać wyniki symulacji, które mają zostać przedstawione na wykresie.
   6. Kliknąć Make figure, aby zwizualizować aktywność impulsową przyczyniającą się do kształtu fali dipolowej.
      UWAGA: Niektóre cechy mikroukładu (np. LFP i CSD) są dostępne wyłącznie poprzez interfejs programistyczny (API) aplikacji HNN-Python. Samouczki oparte na kodzie dla tych funkcji są dostępne na stronie przykładów HNN (https://jonescompneurolab.github.io/hnn-core/stable/index.html).
3. Walidacja predykcji modelu z danymi empirycznymi
   1. Zidentyfikować istniejące zestawy danych i/lub zebrać nowe dane empiryczne (np. inwazyjną elektrofizjologię, laminarną MEG/EEG oraz spektroskopię rezonansu magnetycznego), aby przetestować predykcje wieloskalowe modelu.
   2. Porównać predykcje wieloskalowe modelu z empirycznymi zestawami danych.
   3. Jeśli predykcje wieloskalowe są zgodne z empirycznymi zestawami danych, model uznaje się za zwalidowany dla wybranej cechy mikroukładu.
   4. Jeśli predykcje wieloskalowe nie są zgodne z empirycznymi zestawami danych, zaktualizować domyślną sieć HNN poprzez ograniczenie jej nowymi danymi empirycznymi i powrócić do kroku 3.

Wyniki

W tej sekcji przedstawiono scenariusz, w którym neuroterapeutyk o nieznanym mechanizmie działania jest badany przy użyciu oprogramowania do modelowania HNN. Celem jest wykorzystanie sygnałów EEG sprzed i po zastosowaniu leczenia do wygenerowania prognoz dotyczących tego, w jaki sposób neuroterapeutyk modyfikuje obwody neuronalne. Wyniki są przedstawione w celach demonstracyjnych, aby zilustrować, w jaki sposób modelowanie HNN może zostać zastosowane do badania mechanizmów neuroterapeutycznych.

Opracowywanie hipotez mechanistycznych leżących u podstaw biomarkerów EEG ERP (Krok 1)

W niniejszym przykładzie wykorzystano hipotetyczny sensoryczny paradygmat ERP do zbadania, w jaki sposób neuroterapia zmienia sygnał (Krok 1). Rysunek 1A przedstawia słuchowy ERP przed leczeniem (kolor niebieski) obok hipotetycznego ERP po leczeniu (kolor czerwony; zobacz również Rysunek 9). Słuchowy ERP przed leczeniem to eksperymentalnie zarejestrowane dane z lokalizacją źródła pochodzące z pracy Kohl i wsp.⁴³, natomiast hipotetyczny ERP po leczeniu został wygenerowany poprzez przeskalowanie fali przed leczeniem za pomocą okna z taperowaniem gaussowskim. Jak pokazano, hipotetyczna neuroterapia powoduje znaczny spadek amplitudy komponentów P1, N1 i P2 w stosunku do ERP przed leczeniem.

Należy zauważyć, że w pracy Kohl i wsp.⁴³, z której pobrano dane ERP sprzed opracowania, symulacje HNN wykorzystywały model, w którym neurony piramidalne zostały wzbogacone o bardziej realistyczną dynamikę kanałów wapniowych niż w domyślnym modelu HNN. W rezultacie wyniki symulacji w pracy Kohl i wsp.⁴³ nieco różnią się od wyników przedstawionych tutaj. Model Kohla i wsp. z 2020 roku (oraz inne zaktualizowane modele HNN) są dostępne za pośrednictwem API Python (https://jonescompneurolab.github.io/hnn-core/stable/generated/hnn_core.calcium_model.html#hnn_core.calcium_model). Dostęp do takich rozszerzonych modeli poprzez interfejs GUI jest obecnie w fazie opracowywania.

Następnie należy zidentyfikować parametry modelu reprezentujące efekty związane z leczeniem (tj. parametry zainteresowania), które według hipotezy wyjaśniają, w jaki sposób neuroterapia redukuje amplitudy P1, N1 i P2 (Kroki 1.6–1.8). Szerokie kategorie potencjalnych mechanizmów neuronalnych (i odpowiadających im parametrów modelu) obejmują czas wejść synaptycznych egzogennych, lokalne przewodnictwa kanałów jonowych neuronów, lokalną łączność synaptyczną oraz egzogenną łączność synaptyczną (Rysunek 1B). W tym przykładzie potencjalne mechanizmy z każdej kategorii są oceniane za pomocą HNN, aby zbadać, jak zmiany tych parametrów wpływają na symulowany ERP.

Parametry analizowane

1. Odchylenie standardowe pierwszego (wzgórzowo-korowego) pobudzenia proksymalnego (tj. synchronii wzgórzowo-korowej), reprezentujące zmienność synchronizacji początkowych doprzewodnich sygnałów wejściowych z receptorów zmysłowych.
2. Przewodnictwo kanałów potasowych muskarynowych (K_m) w neuronach piramidowych warstwy 5 (L5), kontrolujące pobudliwość neuronów w taki sposób, że pobudliwość maleje wraz ze wzrostem przewodnictwa.
3. Siła lokalnych receptorów GABA_B, odpowiadająca powolnej synapsie hamującej przekazywanej przez interneurony do wszystkich komórek w lokalnej sieci.
4. Siła przewodnictwa sprzężenia zwrotnego (korowo-korowego) pobudzenia dystalnego, reprezentująca siłę wywołanego sensorycznie sygnału wejściowego sprzężenia zwrotnego o opóźnieniu ~100 ms do synaps AMPA i NMDA w warstwach supragranularnych.

Ustalenie dopasowania modelu ERP przed przygotowaniem (kroki 3–4)

Przeprowadź symulację wstępnego ERP, postępując zgodnie z krokami 3–4 (końcowa symulacja wstępnego przetwarzania przedstawiona na ryc. 8C). Pozytywny wynik sygnalizuje bliskie dopasowanie symulowanych i empirycznych przebiegów fal, co jest określane przez wysoki współczynnik korelacji i niski RMSE.

Ustalenie dopasowania modelu ERP po zakończeniu leczenia (Krok 5)

Użyj modelu ERP sprzed zabiegu jako punktu wyjścia, a następnie zastosuj ręczne dostrajanie i optymalizację parametrów, aby ustalić, czy interesujące parametry mogą odtworzyć empiryczny ERP po zabiegu. Prawidłowe dopasowanie wskazuje, że hipotetyczne parametry są wystarczające do wyjaśnienia zmian w przebiegu fali ERP związanych z zastosowanym zabiegiem.

Kwantyfikacja niepewności z użyciem SBI (Krok 6)

Ze względu na degenerację parametrów nieodłącznie występującą w modelach biofizycznych, kwantyfikacja niepewności przy użyciu SBI (Krok 6) jest niezbędna do tworzenia prognoz dotyczących zmian parametrów przed i po leczeniu. Krytycznym warunkiem wstępnym dla SBI jest uzyskanie dokładnych dopasowań do ERP przed i po leczeniu (Kroki 3–5). Jeśli nie uda się osiągnąć dokładnych dopasowań, próbki a posteriori wygenerowane przez SBI mogą nie odtwarzać empirycznych przebiegów fal, co prowadzi do niewiarygodnych prognoz.

Jeśli w krokach 3–5 nie uda się uzyskać poprawnego dopasowania, należy zweryfikować wybór interesujących parametrów oraz ich zakresy a priori przed zastosowaniem SBI.

W tym przykładzie SBI zostało zastosowane tylko do czterech interesujących nas parametrów po leczeniu, podczas gdy wszystkie pozostałe parametry pozostawiono bez zmian. Choć zastosowanie SBI do większego zestawu parametrów może poprawić odporność, znacząco zwiększa to koszty obliczeniowe (patrz Dyskusja).

SBI jest wykorzystywane do szacowania pełnych rozkładów parametrów, które generują symulowane ERP ściśle odpowiadające falom docelowym. W skrócie, SBI to podejście oparte na wnioskowaniu bayesowskim, które trenuje sieć neuronową w celu zmapowania wyników modelu na rozkłady parametrów modelu^52,53,55. Wytrenowana sieć jest następnie stosowana do empirycznych kształtów fal w celu wywnioskowania rozkładów parametrów zgodnych z danymi. Wymaga to postawienia hipotez a priori dotyczących zakresów parametrów.

W tym przykładzie zdefiniowano jednorodny rozkład a priori dla czterech interesujących parametrów: synchronii wzgórza-kory, przewodnictwa Km dendrytów neuronów piramidowych, lokalnego przewodnictwa GABAB oraz siły sprzężenia zwrotnego kora-kora. Granice a priori zdefiniowano jako wielokrotności wartości domyślnych: 0–5× dla synchronii wzgórza-kory oraz 10⁻¹–10¹× dla pozostałych parametrów.

Rysunek 10A przedstawia wynikowe rozkłady parametrów dla ERP przed i po zastosowaniu leczenia, zwizualizowane za pomocą wykresu par (pairplot). Panele na przekątnej pokazują rozkłady jednowymiarowe, natomiast panele poza przekątną przedstawiają zależności dwuwymiarowe. Przewidywania mechanistyczne odpowiadają parametrom o silnie rozdzielonych rozkładach pomiędzy badanymi warunkami.

Analiza rozkładów jednowymiarowych wykazuje, że synchronizacja wzgórzowo-korowa charakteryzuje się największą separowalnością przed i po leczeniu (najniższy wskaźnik OVL wynoszący 0,07) i wzrasta po zastosowaniu terapii (Rysunek 10A(iii), kolor czerwony). Wskazuje to, że model HNN przewiduje modulację synchronizacji wzgórzowo-korowej jako potencjalny mechanizm działania.

Posteriorna walidacja predykcyjna

Zweryfikuj wnioskowane rozkłady parametrów za pomocą PPC. Wygeneruj niezależne próbki parametrów z rozkładu a posteriori i zasymuluj odpowiadające im ERP. O pomyślnym PPC świadczy fakt, że zasymulowane przebiegi fal ściśle odpowiadają empirycznym ERP.

Jak pokazano na Ryc. 10B oraz Ryc. 10C, przebiegi fal zarówno przed leczeniem (Ryc. 10B, kolor niebieski), jak i po leczeniu (Ryc. 10C, kolor czerwony) ściśle odpowiadają symulacjom wygenerowanym z próbek a posteriori (kolor czarny), przy współczynnikach korelacji wynoszących odpowiednio 0,99 i 0,96 (średnia z 10 niezależnych próbek). Wyniki te potwierdzają, że wywnioskowane rozkłady parametrów pozwalają na dokładną rekonstrukcję przebiegów fal.

Przykład nieudanego PPC przedstawiono na Rysunku uzupełniającym 1. Przykład ten ma taką samą strukturę jak Rysunek 10 i wykorzystuje tę samą wytrenowaną sieć SBI; zastosowano jednak alternatywny przebieg fali po leczeniu, który nie jest dobrze reprezentowany w zbiorze treningowym (np. przebiegi ERP z dodatnim odchyleniem w latencji N1). Nieudany PPC jest wskazany na Rysunku uzupełniającym 1C, gdzie współczynnik korelacji jest niski (np. Corr < 0.95). Co istotne, rozkład a posteriori na Rysunku uzupełniającym 1A wykazuje silnie rozdzielone rozkłady parametrów. Bez przeprowadzenia PPC wyniki te mogłyby zostać błędnie zinterpretowane jako istotne różnice między warunkami przed i po leczeniu. Przykład ten podkreśla znaczenie przeprowadzania PPC równolegle z interpretacją rozkładów a posteriori, ponieważ wyniki z nieudanego PPC są niewiarygodne i nie powinny być dalej analizowane.

Badanie i walidacja modelu (Krok 7)

Dzięki zastosowaniu modelu HNN możliwe jest bezpośrednie badanie i wizualizacja aktywności na poziomie komórkowym i obwodowym, takiej jak aktywność impulsacyjna (spiking), która leży u podstaw każdej symulacji ERP (Krok 7.2.2). Rysunek 11A oraz Rysunek 11B przedstawiają symulowane ERP pobrane z rozkładów parametrów przed i po leczeniu, wraz z odpowiadającą im aktywnością impulsacyjną specyficzną dla komórek (Rysunek 11C oraz Rysunek 11D).

Rysunek 11. Aktywność impulsacyjna na poziomie komórkowym leżąca u podstaw generowania biomarkerów EEG. (A) ERP przed leczeniem (niebieski) z pojedynczą symulacją a posteriori (czarny). (B) ERP po leczeniu (czerwony) z odpowiadającą mu symulacją a posteriori (czarny). (C) Symulowana aktywność impulsacyjna leżąca u podstaw ERP przed leczeniem. (D) Symulowana aktywność impulsacyjna leżąca u podstaw ERP po leczeniu. Kliknij tutaj, aby wyświetlić powiększoną wersję tego rysunku.

Kształty fal są wizualizowane bez wygładzania, aby podkreślić wkład czasu wystąpienia kolców w dipol prądowy. W eksperymentalnych sygnałach EEG duże populacje neuronalne generują sygnały uśrednione przestrzennie, które wydają się gładsze. Ponieważ HNN symuluje mniejszą populację (200 neuronów piramidowych), stosuje się wygładzanie, aby przybliżyć aktywność w większej skali (>100 000 neuronów).

Wyraźną różnicą między warunkami jest zmniejszona aktywność impulsowa w neuronach piramidowych warstwy L5 po zastosowaniu leczenia (Rycina 11C i Rycina 11D, czerwona kropka). Należy zauważyć, że Rycina 11 przedstawia pojedynczą próbkę z rozkładu a posteriori; w celu uzyskania wiarygodnych prognoz należy przeanalizować wiele próbek. Wyniki te dowodzą, że hipotetyczny lek neuroterapeutyczny zmienia wieloskalową aktywność obwodową, co skutkuje zmniejszeniem amplitud P1–N1–P2.

Takie przewidywania można bezpośrednio przetestować za pomocą inwazyjnej elektrofizjologii (np. rejestracji za pomocą sond laminarnych o wysokiej gęstości) lub innych metod obrazowania (Krok 7.3). Nowo pozyskane dane mogą być następnie wykorzystane do dalszego ograniczania przewidywań modelu. Chociaż niniejszy protokół koncentruje się na dopasowywaniu danych EEG w skali makro w celu wnioskowania o aktywności mikroobwodów, strukturę tę można również zastosować w odwrotny sposób, dopasowując dane z mikroobwodów (np. potencjały czynnościowe, LFP/CSD) w celu wnioskowania o sygnałach EEG w skali makro.

Rycina uzupełniająca 1. Przykład nieudanego posterior predictive check w procesie SBI. Wykresy są zorganizowane analogicznie do Ryciny 10. Dane przed leczeniem (niebieskie) są identyczne jak na Rycini 10. Hipotetyczne dane po leczeniu zostały wygenerowane w ten sam sposób co wcześniej (forma fali pomnożona przez okno z zanikiem gaussowskim), ale przekształcone tak, aby uzyskać dodatni szczyt, który nie jest dobrze reprezentowany w zbiorze treningowym symulacji HNN. (A) Wizualizacja pairplot rozkładów parametrów oszacowanych przy użyciu SBI. Panele diagonalne (i–iv) pokazują rozkłady jednowymiarowe dla poszczególnych parametrów, w tym (i) synchronii wzgórzowo-korowej, (ii) przewodnictwa K_m w dendrytach, (iii) przewodnictwa GABA_B oraz (iv) siły sprzężenia zwrotnego kora-kora. Rozkłady dla warunków przed leczeniem (niebieskie) i po leczeniu (czerwone) wykazują wysoką rozdzielność dla wszystkich parametrów (OVL < 0.1). Panele pozadiagonalne pokazują zależności dwuwymiarowe między parametrami. (B) Posterior predictive check (PPC) dla ERP przed leczeniem; symulowane formy fal (czarne) ściśle odpowiadają danym empirycznym (niebieskie). (C) PPC dla ERP po leczeniu; symulowane formy fal (czarne) są bardzo odmienne od danych empirycznych (czerwone), a wartość Corr < 0.95 wskazuje na nieudany PPC.Kliknij tutaj, aby pobrać ten plik.

Dyskusja

Modelowanie neuronowe obliczeniowe biomarkerów EEG może umożliwić głębsze zrozumienie, w jaki sposób terapeutyki OUN rekonfigurują obwody neuronalne i zapewnić przewidywania dotyczące procesów biologicznych leżących u podstaw efektów terapeutycznych. Prezentowany tutaj przepływ pracy pokazuje, jak powszechnie mierzone biomarkery EEG, czyli wyniki ewolucjonowanych potencjałów słuchowych (ERP), wraz z modelowaniem biofizycznym przy użyciu HNN, mogą być wykorzystane jako okno na mechanizmy, przez które lek wpływa na aktywność neuronalną. Łącząc makroskalowe pomiary EEG z leżącymi u podstaw procesami komórkowymi i obwodowymi, ten protokół zapewnia strukturalnie ukierunkowany i oparty na hipotezach rachunek dla interpretacji mechanizus. Co ważne, podejście to nie ogranicza się do ERP i może być rozszerzone na badanie innych lokalnych sygnałów EEG, w tym niskoczęstotliwościowych oscylacji neuronalnych^40,63 i przejściowych zdarzeń spektralnych^7,47,64, rozszerzając tym samym jego zastosowanie na różne biomarkery elektrofizjologiczne i paradygmaty eksperymentalne.

W porównaniu z innymi ramami modelowania neuronowego EEG, HNN oferuje równowagę między złożonością modelu a wydajnością obliczeniową, co jest szczególnie korzystne dla iteracyjnego testowania hipotez. Na przykład The Virtual Brain umożliwia symulację dużych sieci mózgowych, które generują przestrzenno-czasowe sygnały EEG^34,65. Jednakże, aby osiągnąć modelowanie całego mózgu, aktywność neuronalna jest reprezentowana za pomocą uproszczonych formuł matematycznych, które eliminują szczegółowe cechy komórkowe, takie jak morfologię neuronów piramidalnych i ograniczają zdolność do bezpośredniego łączenia parametrów modelu z komórkowymi mechanizmami działania leków. Z drugiej strony, duże skalowe morfologicznie i fizjologicznie szczegółowe modele mogą symulować sygnały EEG z dużą rzeczywistością biologiczną^66,67,68,69, ale przy dużych kosztach obliczeniowych, często wymagających kilku godzin obliczeń na symulację zaledwie kilku sekund aktywności neuronalnej. To obciążenie obliczeniowe może ograniczać dostępność i spowolnić iteracyjny proces wymagany do generowania i testowania hipotez. HNN zajmuje położenie pośrednie (Rysunek 2), umożliwiając symulację zlokalizowanego obwodu kory nowej z wystarczającą ilością biologicznych szczegółów, aby generować prognozy na poziomie komórkowym i obwodowym, jednocześnie zachowując wydajność obliczeniową (czyli symulacje na poziomie sekund), co sprawia, że jest dobrze przystosowany do integracji z warunkami eksperymentalnymi.

Mimo tych zalet, kilka ograniczeń należy wziąć pod uwagę podczas stosowania EEG i modelowania neuronowego biofizycznego do badania chorób mózgu i mechanizmów działania leków. Właściwości komórkowe i obwodowe, które generują sygnały EEG, nie uwzględniają pełnego spektrum procesów biologicznych, na które wpływają interwencje farmakologiczne. Na przykład, reakcje układowe lub immunologiczne mogą nie wpływać bezpośrednio na sygnały EEG i dlatego mogą nie odzwierciedlać się w wynikach modelu. Ponadto, hipotezy mechanizus są często pochodne z badań na zwierzętach, które mogą nie być w pełni przekładane na funkcjonowanie ludzkiego mózgu, szczególnie w zaburzeniach neuropsychiatrycznych, gdzie wyniki kliniczne są oparte na ocenach behawioralnych i poznawczych^70,71. Innym ważnym wyzwaniem jest rozróżnienie między ostrymi a przewlekłymi efektami farmakologicznymi. Podczas gdy ostre interakcje leku-receptora są stosunkowo dobrze scharakteryzowane, długoterminowe adaptacje wywołane przedłużoną ekspozycją na leki są mniej dobrze zrozumiane i mogą nie być w pełni uchwycone w obecnych ramach modelowania. Co więcej, model HNN reprezentuje pojedynczą zlokalizowaną kanoniczną sieć kory nowej, podczas gdy neuroterapeutyki i choroby OUN często wywierają rozproszone efekty na wielu obszarach mózgu. Chociaż wpływy z innych regionów mogą być przybliżone poprzez zmiany w czasie i sile exogenych wejść, bezpośrednie empiryczne charakteryzowanie tych obwodów układowych jest często ograniczone, co ograniczając interpretację modelu.

Degeneracja parametrów stanowi fundamentalne wyzwanie we wszystkich biofizycznych modelach neuronowych, ponieważ wiele konfiguracji parametrów może wytwarzać podobne wyniki modelu. W tym protokole, SBI jest używane do rozwiązania tego problemu poprzez oszacowanie dystrybucji parametrów, które generują fale ERP zgodne z danymi empirycznymi (Rysunek 10). To podejście umożliwia sprecyzowanie niepewności w parametrach modelu, zapewniając bardziej solidny rachunek dla interpretacji mechanizus niż estymacje jednoprzedziałowe. Jednak ze względu na trudności obliczeniowe, SBI jest stosowany do ograniczonego podzbioru parametrów odpowiadających hipotezom mechanizmów działania leków, a założenia dotyczące nieoszacowanych parametrów mogą wpływać na wynikającą dynamikę sieci. Rozszerzenie wywnioskowania do większych przestrzeni parametr

Materiały

Lista materiałów użytych w tym artykule

Nazwa	Firma	Numer katalogowy	Komentarze
Anaconda Python	Anaconda, Inc.	N.D.	Dystrybucja Pythona; Wersja Pythona ≥3.9 i <3.14
Komputer stacjonarny	N.D.	N.D.	System operacyjny: Windows ≥10, Linux lub macOS. Minimalne zalecane sprzęt: ≥16 GB RAM, ≥8 rdzeni CPU
EEGLAB	EEGLAB Developers	N.D.	Opcjonalne narzędzie MATLAB do przetwarzania EEG i analizy ERP
FieldTrip	Donders Institute for Brain, Cognition and Behaviour, Radboud University	N.D.	Opcjonalne narzędzie MATLAB do analizy EEG/MEG
Human Neocortical Neurosolver (HNN-core)	HNN Developers	N.D.	Biofizyczne oprogramowanie do modelowania neuronów; Wersja ≥0.6.0 używana w tym badaniu
MATLAB	MathWorks	N.D.	Wymagane do uruchomienia EEGLAB i FieldTrip (jeśli jest używane)
MNE-Python	MNE Developers	N.D.	Używane do przetwarzania EEG i lokalizacji źródeł
NumPy	NumPy Developers	N.D.	Używane do obliczeń numerycznych i generowania liczb losowych
Pixi (menedżer pakietów/środowiska)	Prefix.dev	N.D.	Używane do zarządzania zależnościami w powiązanym repozytorium kodu
PyTorch	PyTorch Developers	N.D.	Używane do trenowania sieci neuronowych SBI i ustawiania losowych ziaren
Pakiet SBI (Simulation-Based Inference)	SBI Developers	N.D.	Pakiet Pythona do estymacji parametrów i ilościowania niepewności
Windows Subsystem for Linux (WSL2)	Microsoft Corporation	N.D.	Wymagany tylko dla instalacji opartych na Windows