ИЗТЕГЛЕТЕ PDF
ABSTRACT. Advanced neural network simulation packages such as NEST and NEURON have been adapted for high performance computer systems. However, these tools rely on differential equations and interpolation methods to describe the dynamics and changes in the state of neurons and synapses, which are computationally and communication intensive. Simulations performed with the NEST of the supercomputer K and JUQUEEN take about 40 minutes of computational time to simulate 1 second activity of rat visual cortex neurons. To significantly reduce simulation time, new mathematical tools and a different organization of the computational process are needed to work effectively across a variety of multi-core hardware platforms. This article introduces an original new technique that effectively utilizes the capabilities of Exascale Supercomputer and aims to reduce computing and communication times, as well as the need for high memory that comes along with most basic simulation tools. In principle, the new procedures included in the proposed algorithm are described.
First, we used coarse and fine-grained procedures, which only allows us to accurately calculate the Hodgkin-Huxley model once. It describes the changing electrical characteristics of the neural membrane. The resulting solution is repeatedly used to calculate the local membrane potential generated by each synapse. It has been experimentally found that the membrane potential is the sum of local (individual) potentials generated by the axons associated with each cell.
Second, the cell in our simulation is represented as a structure consisting of all signals coming from all associated with its axons generated at a fixed time interval. This allows us to represent the input of each cell as a vector consisting of all spikes sent about all the related neurons. This representation allows each cell to be separated as a single unit and all the neural network cells to be modeled simultaneously and independently of one another, without exchanging information between them during the simulation interval. Generally, within one interval, the modeling process consists of three consecutive phases – calculating the potential of the cell from the direct signals to it; to it is added the potential resulting from the effects of the feedback signals and the adjacent group of cells. The neurons of the network are grouped together in a similar way, grouped in the brain (Neocortical columns).
The effect of describing the network with the proposed topology is that the modeling process of the equivalent silicone “model” consists of multiple independent paralyzed (threads) with data exchange only at the end of the simulation interval.
OPENMP and one-sided MPI communications were used in parallelization to optimize the overlap between calculation and data communication between cores in one processor and between supercomputer processors.
Because the simulation is not driven by time or event, but is instead driven by the cell, the required simulation time for the network drops sharply and therefore the actual performance of the processor cores is close to their theoretical performance. This in turn provides the ability to perform multi-threaded and simultaneously large-scale neural simulation.
Keywords: Hodgkin-Huxley model, large-scale simulations, massive parallel computing, spiking neural network, multi-core architecture, high-performance machines.
1. Въведение
Съвременната невробиология комбинира теоретични и експериментални подходи с изчислително моделиране, за да изучи свойствата на нервната система в различни пространствени скали на човешката мозъчна организация, за да се предвиди как молекулите, клетките и невронните мрежи генерират човешкото поведение. За целта са разработени мощни изчислителни симулационни среди като NEST (Gewaltig и Diesmann, 2007), NEURON (Hines и Carnevale, 2006), GENESIS (Bower и Beeman, 2012) и BRIAN (Stimberg et al., 2014). Всеки софтуерен пакет набляга на различни пространствени нива на биофизично описание на невроните и синапсите, за да се създадат еквивалентни модели на мозъка. Прецизното моделиране донякъде позволява да се проследят процесите в мозъка и да се изучи неговата биофизика, но то е пряко свързано със симулации на мрежи с десетки до стотици милиони неврони и синапси, което поставя уникални изисквани към изчислителния процес – разпределението му (мащабираност) между десетки хиляди процесори, разпадането му на няколко стотин хиляди паралелни процеса (мащабен фин изчислен паралелизъм) и слаба свързаност между тях, (Frings et al., 2009). Например, високите разходи за изчисление и памет произтичат от високата свързаност и сложното свързване на невронните мрежи (Jordan et al., 2018; Ippen et al., 2017). От комуникационна гледна точка, както посочват Morrison et al., времето за отчитане на събитията става значително с увеличаването както на броя на изчислителните възли, така и на симулираните неврони (Morrison et al., 2005).
Стремежът е да се проектират много големите невронни мрежи на мозъка върху функционално еквивалентна софтуерна среда, която да се имплантира в специално проектирана за нея суперкомпютърна платформа и да се създаде високоскоростен паралелен алгоритъм за тяхното моделиране с производителност, по-голяма от тази на мозъка. Това е съвместната цел на два от най-големите Европейски проекти – Partnership for Advanced Computing in Europe (PRACE)[1], чиято цел е възможности за научно и инженерни изследвания и разработки с високо въздействие във всички дисциплини за повишаване на европейската конкурентоспособност в полза на обществото чрез осигуряването на високопроизводителни ресурси за изчисляване и управление на данни, както и Human Brain Project (HBP)[2], чиято цел е изграждане на изследователска инфраструктура за подпомагане на невронауката, медицината и компютърните технологии. PRACE представлява сдружение на 47 университета и изследователски центрове от 27 Европейски държави (включително България). HBP пък от своя страна обединява 116 партньора от 21 европейски страни и още 2 две страни извън ЕС и с бюджет от над 1.2 милиарда евро за периода 2015-2020 година.
Като продължение на тази работа, в тази статия представяме нов алгоритъм за ефективна симулация на изключително големи невронни мрежи, който да намали изчислителните и комуникационни времена, както и съществените изисквания за памет, необходими на такива симулации.
Настоящата статия е организирана по следния начин. Започваме с кратък преглед на физиологичната структура на неврона за осигуряване на невробиологичната основа за модела. След това е представено въвеждането на модел на описание на неврона на Ходжкин-Хъксли, използван в настоящата работа. След това описваме нашия метод и организация на данните. Завършваме със заключения и бъдещи направления.
2. Кратък преглед на структурата и функцията на неврона
В следващото кратко описание се фокусираме върху структурните характеристики на невроните. Въпреки че съществуват различни видове неврони, които се различават по своята морфология, молекулярна структура, генна експресия и други характеристики, всички те се състоят от клетъчно тяло, аксони и крайни разклонения на аксона (дендрити) (фигура 1).
По своята форма невроните се делят на сферични, зърнести, пирамидални, звездовидни и вретенообразни. Тяхната големина варира от 5 мкм (10-6 метра) за малките зърнести клетки до 120 -150 мкм за големите пирамидални неврони. Последните са над 80% от клетките на човешкия мозък. По морфологичен признак, невроните се разделят на: мулти-полярни (с един аксон и няколко дендрита) -като преобладаващата част от клетките в неокортекса са от този вид; еднополярни (с един израстък) – този тип неврони се намират предимно в сензорната част на троичния нерв (trigeminus) в средния мозък; псевдо еднополярни, които са групирани в гръбначния мозък и ганглиите между прешлените; и накрая са биполярни (с един аксон и един дендрит), които се намират в сетивните органи.
 |
| Фигура 1: Структура на неврона. 1 – тяло; 2 и 5 – крайни разклонения на аксона (дендрити ); 3 – аксон; 4 – миелинова обвивка. |
Невронът се състои от тяло с клетъчно ядро (soma) и два вида израстъци: дендрити и аксон. Между влизащите аксони и дендритите има междина, която се нарича синапс. На повърхността си пресинаптичната мембрана на аксона има много йонни канали. Аксоните могат да варират в диаметъра си и миелиновата обвивка. Те пренасят определен тип електрически импулс, наречен потенциал на действие. Последният предизвиква локална промяна в потенциала на клетъчната мембрана. По време на състояние на покой на неврон, електрическият потенциал през мембраната е около -70 mV. Като сума от локалните потенциали мембранният потенциал може да достигне стойност над 110 mV. След това клетката излъчва импулс към другите свързани с нея неврони (фиг.2), потенциалът на мембраната бързо намалява и достига стойността му на покой. По време на покоя невронът не реагира на сигналите , пренесени към него от синапсите.
Един аксон може да се свърже с много други неврони и да предизвика отговор във всички тях едновременно. Обикновено синапсите предават сигнал по посока – от пресинаптичната мембрана към постсинаптичната мембрана. Преобразуването на сигнала се осъществява чрез освобождаване на невротрансмитери в синаптичната цепнатина. Под негово въздействие някои везикули се отварят и да отделят невротрансмитери, които се свързват се с рецепторите на мембраната на неврона. Това променя йонната пропускливост клетъчна мембрана, като се предизвиква каскада от обмяна на йони от вътрешността на клетката към средата около нея и обратно.
Фигура 2. Потенциалът за действие индуцира освобождаването на невротрансмитери през синаптичната цепнатина, генерира електрически сигнал в постсинаптичния неврон. (Снимка от: (Thomas Splettstoesser, 2015) CC BY-SA 4.0).
3. Модел на неврона
Един от най-известните модели, описващи биофизичните процеси, които протичат във възбудения неврон е модела на Ходжкин-Хъксли (Hodgkin and Huxley, 1952). Той се счита за един от най-биологично реалистичните модели в изчислителната невронаука. Представлява нелинейна система от обикновени диференциални уравнения (ODE), състояща се от няколко уравнения, които отчитат динамиката на йонните канали в клетъчната мембрана и изменението на тяхната проходимост в зависимост от разликата на потенциала между нейната вътрешна и външна стена (трансмембранен потенциал). Например, натриевите йони са повече извън клетката, калиевите йони са повече вътре в нея. Потенциалът на мембраната се определя от разликата на концентрацията на натриеви и калиеви йони вътре в клетката и в околната среда, като:
където 
представлява ток, докато С е мембранният капацитет. Структури 
представляват максималната проводимост на натриевия, калиев канал и хлорния канал, който формира тока на утечка. 
and 
са потенциалите, предизвикани от промяната на концетрация за натриеви и калиеви йони, докато 
е обратен потенциал за проводимостта. Пропускащите променливи 
служат като регулатори на проходимостта на йонните канали и са нелинейна функция на потенциала на мембраната :
Зависимите от напрежението скорости на отваряне и затваряне са дадени като:
Като експериментално измерените константи на мембраната са:
Фигура 3. Потенциал за действие на Ходжкин-Хъксли. Фигура от: (Bower and Beeman, 2012)
4. Моделиране биологични невронни мрежи
Симулация на цифров аналог на анатомична невронни мрежи на суперкомпютър Blue Gene/P с 8000 процесора, е осъществена през 2009 година. Моделът включва 31 000 неврона от 17 вида, които помежду си има 37 милиона синапси. Това представлява една колона от зрителния кортекс на плъх. Потенциалът на мембраната на всеки неврон е изчислен посредством уравненията на Ходжкин-Хъксли. Целият модел, създаден от екип от над 80 изследователя под ръководството на проф. Хенри Маркрам, е базиран на резултати от 30 годишни изследвания в областта на физиологията, морфологията, анатомията и клетъчната биология на невроните и на микроструктурите на кортекса (Markram, 2015).
 |
| Фигура 4. Мрежа с 31 хил. неврона и 37 мил. Синапса. Фигурата е взета от Markram, Henry, et al. “Reconstruction and simulation of neocortical microcircuitry.” Cell 163.2 (2015): 456-492. |
Нарастването на производителността на суперкомпютрите позволява през 2012 година да се моделират 100 колони от кортекса, съставени общо от около 106 неврона и 109 синапса. Като пример може да се посочи, че този брой неврони и синапси съставят размера на мозъка на пчелата. През 2013 пък е успешно моделирана невронна мрежа състояща се от 1,86 милиарда неврона и 11,1 трилиона синапса на суперкомпютъра в Япония – К-computer, който има 82,944 процесора ,663,552 ядра, както и на Blue Gene/Q, намиращ се в Суперкомпютърния изследователски център в Julich, който пък се състои от 24,576 процесора (393,216 ядра) и оперативна памет от 3.93 Терабайта (Diesmann, 2013). Симулираната мрежа на двата суперкомпютъра представлява мрежата на зрителния център (visual cortex) на плъха, която е с аналогична изчислителна сложност. За това са използвани два вида софтуер – японския Spiking network simulation code for petascale computers, както и немския пакет Neural Simulation Tools (NEST). За моделирането на една секунда от активност на анатомична невронна мрежа на плъха са необходими 40 минути машинно време, което прави невъзможна нейната симулацията в реално време (Kunkel, 2012).
Причината е, че всеки импулс, пренесен от аксона през синапса до мембраната на неврона, локално променя потенциала й. Тъй като в рамките на 300 милисекунди преминават между 100-150 сигнала, за да се изчисли стойността на възбудения от всеки от сигналите мембранен потенциал, е необходимо да се реши системата от четири диференциални уравнения. Симулацията на 1 секунда биологично време на мрежа с 11.1 трилиона синапса налага решението на 1.11 и 1.65 x 1015 системи. Съответно това представлява между 26 и 40×1018 операции, което е лимитиращ фактор, тъй като дори на най-мощните суперкомпютри те отнемат 40 минути изчислително време.
В следващия раздел е описана нашия модел , новата му структура и нашия оригинален алгоритъм за неговото решаване. Той се базира на по-внимателно проследяване на физиологичните реакции във веригата аксон->синапс->мембрана на биологичния неврон.
5. Нов високоскоростен алгоритъм за моделирането на големи биологични невронни мрежи.
Преди да се инициализира симулацията, е необходимо дефиниране на мрежовата топология, което може да се мисли като графика, така че, G = (V,E) където V и E съответно представляват невронните и връзките между тях. Допълнителното описание на мрежовата топология включва няколко характеристики като принадлежност на неврони към определена популация (или колона от кортекса, ако такава е необходима), нейният глобален идентификатор и дали е възбуждащ или инхибиращ. Описани са и специфичните връзки между невроните – с глобален идентификатор, теглото и времето на забавяне на предаването, което може да бъде различно за всеки неврон.
По този начин всеки неврон 
се представя като независим обект. Той се описва чрез структура, в който са записани векторите на спайкове 
които е получил в предишния интервал от време с дължина 300 ms от свързани с него неврони 
. Като 
, където M представлява броя на неврони в мрежата, 
е броя на синапсите, през които сигналите от аксоните достигат i клетка. В следващия интервал с дължина от 300 милисекунди се изчисляват едновременно и независимо потенциалите на мембраните на всички клетки. Стъпката на изчисление е 0.1 ms. По този начин алгоритмът за моделиране на мрежата се разделя на хиляди паралелни и независими изчислителни процеси.
Синапа представлява структура съдържаща информация за синоптичното тегло (
), времето на закъснение на сигнала пренасян от аксон 
, свързващ клетка 
с клетка 
, която изпраща сигнала. Структурата съдържа и моментите на възникване на спайковете 
,които се изпращат през аксона от клетката i до j клетка през синапса в интервала 0 – 300 милисекунди и е записан като вектор във вид:
където 
е броят на спайковете, пренесени от i клетката до j клетката.
От модела на Ходжкин-Хъксли (HH), с предварително зададена точност, се изчисляват 60 дискретни стойности на потенциала на мембраната и се запомнят в паметта.
За всички спайкове, изпратен до клетката от клетки, свързани с нея, се изчислява предизвиканото от тях локално изменение на потенциала на мембраната на клетка, където p=1 за сигнал от стимулираща клетка и p=-1 за сигнал от задържаща (inhibitory) клетка.
Следва сумирането на всички локални потенциали и се изчислява глобалния потенциал на мембраната на клетката в 3000 точки в интервал 0-300 милисекунди (стъпка 0,1 милисекунда):
След изчислението на глобалния потенциал се проверява дали в някоя от точките е преминат прага (threshold) при който се генерира сигнал, заложен за всяка клетка от -50V. Ако този праг е преминат се записва времето на възникване на импулса във вектор.
При симулацията трябва да се изчислят и сигналите, които се предават от рекурентна връзка, както и тези, които са от синапси, които са външни за колоната. Това се извършва при стъпката за изчисление на глобалния потенциал. След записването на новите сигнали следва и преди обновяването на масивите на синапсите към изчислението на глобалния потенциал се прибавят сигналите идващи от рекурентните връзки, както и тези от външните. Следва преизчисляване на всеки от глобалните потенциали и отново се проверява отново дали сигнала е преминал предварително определения праг и вече се записват генерираните импулси. С генерираните спайкове се обновяват масивите на синапсите с новите и се преминава към следващите 300 милисекунди.
Цялостното описание на алгоритъма може да се резюмира по следния начин:
- Само веднъж се решава модела на Ходжкин-Хъксли. Решението е еталонен потенциал и се записва в паметта.
- Всеки 300 ms се обновяват масивите на синапсите. Старите масиви се заместват с новите, в които са записани всички сигнали (spikes), генерирани от всички клетки на мрежата в интервала от 300 ms.
- Локалният потенциал, породен от всеки импулс, се изчислява като еталонният потенциал и се умножава по интензитета на импулса . Трябва да се отчете и момента, в който сигналът е достигнал мембраната. Не се решава системата диференциални уравнения.
4. Сумират се локалните потенциали и ако сумата надвиши определен праг, клетката генерира импулс (spike). Моментът на генерация се записва във вектора S(i,j). Интервалът е 300 ms. На фигура 5 е показана връзка между два неврона и вектора на времената на възникване на импулсите. - Връщане в точка 2.
Нашият метод драстично намалява обем на изчисленията по сравнение с този на класическия алгоритъм , описан в секция 3. За да се изчисли един локален потенциал са необходими 40 събирания и умножения. За една секунда, един аксон пренася между 100 и 150 сигнала. За да се симулира една секунда от работата на мрежа с 11.1 трилиона синапса трябва да се извършат от 89 до 133 x 10-15 операции. Това са между 220-260 пъти по-малко операции от тези на досега използваните алгоритми. При тях на всяка стъпка се решава уравнение на Ходжкин-Хъксли по метода на Рунге-Кута от четвърти ред (над 260 операции).
 |
|
Фигура 5. На фигурата е изобразена връзка между два неврона и вектора на времената на възникване на импулсите |
Също така глобалните потенциали на всички неврони се изчисляват паралелно и независимо един от друг. Изчислителните процеси не обменят данни по между си, тоест няма комуникации. Практически се достига максималната производителност на машината.
Посредством представеният модел невроните се обособяват като независими обекти, представянето на връзките между тях като матрица и записа на импулсите (spikes) като вектори, ни позволява бързо да обновим информацията в мрежата. Също така, на всеки 300 ms старите вектори се заместват с новите, като се използва локалния идентификатор, глобалния идентификатор и новия вектор (показани на фигура 6).
|
| Фигура 6. Представяне на неврона, връзките и спайковете |
Заключение
Настоящата статия представи нов метод за симулирането на големи невронни мрежи. Специфичната структура на данни отделя невроните, връзките и шиповете, необходими за изчисляване на локалния и глобалния потенциал, образувайки независими структури, които да бъдат реализирани и изчислявани паралелно и независимо. Предимство на предлагания метод е, че елементите от изчислителния модел могат да бъдат реализирани в различни хардуерни архитектури – например нови процесори като AMD EPYC2 или Fujitsu A64FX, както и графични процесорни единици (GPUS).
ФИНАНСИРАНЕ
Тази работа е финансово подпомогната от НЦСП-България и PRACE-5IP.
БИБЛИОГРАФИЯ
Bower, J. M. and Beeman, D. (2012). The book of GENESIS: exploring realistic neural models with the GEneral NEural SImulation System (Springer Science & Business Media)
Frings, W., Wolf, F., and Petkov, V. (2009). Scalable massively parallel i/o to task-local files. In High Performance Computing Networking, Storage and Analysis, Proceedings of the Conference on (IEEE), 1–11
Fujita, K., Okuno, S., and Kashimori, Y. (2018). Evaluation of the computational efficacy in gpu-accelerated simulations of spiking neurons. Computing , 1–20
Gewaltig, M.-O. and Diesmann, M. (2007). Nest (neural simulation tool). Scholarpedia 2, 1430
Hines, M. L. and Carnevale, N. T. (2006). The neuron simulation environment. NEURON 9
Hodgkin, A. L. and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. The Journal of physiology 117, 500–544
Ippen, T., Eppler, J. M., Plesser, H. E., and Diesmann, M. (2017). Constructing neuronal network models in massively parallel environments. Frontiers in neuroinformatics 11, 30
Jarrell, T. A., Wang, Y., Bloniarz, A. E., Brittin, C. A., Xu, M., Thomson, J. N., et al. (2012). The connectome of a decision-making neural network. Science 337, 437–444. doi:10.1126/science.1221762
Jordan, J., Ippen, T., Helias, M., Kitayama, I., Sato, M., Igarashi, J., et al. (2018). Extremely scalable spiking neuronal network simulation code: from laptops to exascale computers. Frontiers in neuroinformatics 12, 2
Lechner, M., Grosu, R., and Hasani, R. M. (2017). Worm-level control through search-based reinforcement learning. arXiv preprint arXiv:1711.03467
Morrison, A., Mehring, C., Geisel, T., Aertsen, A., and Diesmann, M. (2005). Advancing the boundaries of high-connectivity network simulation with distributed computing. Neural computation 17, 1776–1801
Paugam-Moisy, H. and Bohte, S. (2012). Computing with spiking neuron networks. In Handbook of natural computing (Springer). 335–376
Preissl, R., Wong, T. M., Datta, P., Flickner, M., Singh, R., Esser, S. K., et al. (2012). Compass: A scalable simulator for an architecture for cognitive computing. In Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis (IEEE Computer Society Press), 54
Stimberg, M., Goodman, D. F., Benichoux, V., and Brette, R. (2014). Equation-oriented specification of neural models for simulations. Frontiers in Neuroinformatics 8
[Dataset] Thomas Splettstoesser (2015). Synapse schematic. [Online; accessed September 30, 2018]
Llort, Germán, et al. “On the usefulness of object tracking techniques in performance analysis.” Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis. ACM, 2013.
Markram, Henry, et al. “Reconstruction and simulation of neocortical microcircuitry.” Cell 163.2 (2015): 456-492.
Diesmann, Markus. “The road to brain-scale simulations on K.” Biosupercomput. Newslett 8.8 (2013).
Kunkel, Susanne, et al. “Meeting the memory challenges of brain-scale network simulation.” Frontiers in neuroinformatics 5 (2012): 35.
[1] http://www.prace-ri.eu/
[2] https://www.humanbrainproject.eu/en/
Leave a Reply