CC BY
12
Вычислительные технологии
Том 20, № 1, 2015
Оптимизация вычисления концентрации вещества в уравнении диффузии с источником, зависящим от времени
Я. А. Сулейманов1'*, Ф.М. Гафаров1, Н. Р. Хуснутдинов1, Н.А. Емельянова2 1 Казанский федеральный университет, Россия
2Казанский филиал Московского государственного университета путей сообщения, Россия *Контактный e-mail: [email protected]
Разработан метод сокращения числа итераций вычисления концентрации вещества. Метод применён для численного решения уравнения диффузии с источником, зависящим от времени. Предложенный подход представляет собой инструмент, с помощью которого можно сократить время численного анализа на 13-26 %.
Ключевые слова: уравнение диффузии, оптимизация вычислений, нейронные сети.
Введение
Численное моделирование сложных систем, состояния которых в каждый момент времени зависят от множества параметров, представляет собой достаточно сложный вычислительный процесс. Для оптимизации таких задач существуют различные методы: оптимизация алгоритма программного кода, перевод алгоритма из последовательной в параллельную парадигму вычислений, упрощение самой задачи.
В представленной работе приводится метод, позволяющий сократить время численных экспериментов при моделировании эволюции нейронных сетей. В качестве эволюции нейронных сетей рассматривается структурная пластичность в нейронных сетях. Структурная пластичность связана с анатомической структурой нейронов и определяется связями между нейронами [1]. Рост и ветвление аксонов и дендритов приводят к анатомическому изменению нейронных сетей. Моделирование роста аксона каждого нейрона представляет собой сложный и длительный процесс, требующий использования громоздких вычислений. Ключевым явлением для роста аксона является диффузия вещества AGM (axon growth molecules) [2-4], которое испускается нейронами и воспринимается кончиком аксона как управляющий сигнал. Для того чтобы в заданной точке определить значение концентрации AGM, используется стандартное уравнение диффузии. Поскольку численное решение уравнения диффузии с источником, зависящим от времени, представляет собой сложный процесс, то необходимость оптимизации работы программы появляется ещё на стадии компьютерной реализации математической модели. Для проведения необходимой оптимизации разработан метод, который позволяет анализировать и сокращать число итераций при вычислении концентрации вещества AGM.
© ивт со ран, 2015
1. Объект исследования
Объект исследования — модель, описывающая динамическое образование межнейронных связей за счёт вещества, выделяемого самими нейронами с учётом ветвления аксонов. Для описания распространения вещества АОМ в пространстве в модели используется стандартное уравнение диффузии. Концентрация вещества АОМ сг, выделившегося из ¿-го нейрона, удовлетворяет уравнению диффузии
дс ■
- П2Аасг + ксг = Зг{гз,1). (1)
Здесь Б2 — коэффициент диффузии, к — коэффициент деградации, , ¿) — функция источника, зависящая от времени t. В качестве источника в модели рассматривается нейрон. Решением уравнения (1) является выражение для концентрации
сг(г,- - тг,г) = а <ИнСй(Г] - - tk),
где jг(tk) — активность нейрона, а — количество вещества АОМ, выделяющегося из нейрона за единицу времени, - гг^ - tk) — функция Грина, которая выглядит
следующим образом:
1 2
^' = («1^е-и-^ (3)
В модель включён механизм ветвления аксона. Ветвление аксона происходит путём образования новой ветви из оси аксона [5-8]. Для определения точки ветвления на оси аксона и его начального направления роста необходимо иметь данные о значении концентрации АОМ в одной из множества точек на оси аксона и в точках на сфере вокруг точки ветвления.
Уравнение движения аксона нейрона, описываемого радиус-вектором gn, имеет вид
= А^'г - jШ) £ Vcfc(gn - гк, (4)
к= 1
где А — коэффициент, описывающий чувствительность и подвижность аксона, в^г -jth) — ступенчатая функция, зависящая от активности нейрона. Пороговый параметр j1Н определяет движение аксона. Аксон движется, если активность jг нейрона больше порогового значения активности jИндекс п означает номер конуса роста (номер ветви). В начальной момент времени t нейрон обладает одним конусом роста п =1 на аксоне нейрона.
При проведении численного анализа моделируются нейронные сети, состоящие из множества нейронов, расположенных в трёхмерном пространстве. Нейрон обладает активностью, которая влияет на рост аксона [9,10]. При высоких значениях активности нейрона аксон не растёт и, наоборот, при низких значениях — происходит его рост. Таким образом, система из множества нейронов со временем образует между собой связи. Аксон нейрона растёт в направлении наибольшего градиента концентрации АОМ. Результаты численных экспериментов были опубликованы в работах [11-13].
t
2. Теоретическое описание метода сокращения числа итераций
В выражении (2) с каждым увеличением времени t на значение At число итераций при вычислении концентрации AGM увеличивается на значение At/tstep, где tstep — шаг интегрирования. Значение концентрации AGM вычисляется отдельно для каждой ветви аксона. При проведении численного анализа из-за процесса ветвления аксона число ветвей с течением времени возрастает, что вызывает увеличение числа итераций при вычислении концентрации AGM. Таким образом, при вычислении концентрации AGM на конце ветви нейронной сети, состоящей из N нейронов, каждый из которых на аксоне имеет N,r ветвей, число итераций равно
N NNbr t
Niter = —-, (5)
tstep
где Nbr — число ветвей, на которых происходит ветвление. Данное выражение показывает число итераций, необходимое для вычисления концентрации AGM в момент времени t. Необходимо учитывать, что рост аксона может длиться днями. Например, для нейронной сети, состоящей из восьми нейронов с одним кончиком аксона, которая развивается в течение двух дней, т.е. t = 172 800с, tstep = 1 с, а число итераций для вычисления концентрации AGM в момент времени t равно 1 382 400. Здесь время t измеряется в секундах, Nbr = 1. Эволюция нейронной сети может длиться более двух дней, что соответственно потребует большего времени вычислений.
Кроме движения аксона в модели реализован процесс ветвления, который связан с вычислением концентрации AGM. В отличие от вычисления концентрации AGM во время роста аксона, где расчёт производится в одной точке, в процессе ветвления для определения направления роста новой ветви вычисления ведутся на множестве точек, расположенных на сфере
м NNiatNiotNbrt
Niter = -1-, (6)
tstep
где Niat — число точек по широте, Niot — число точек по долготе.
3. Сокращение числа итераций при вычислении концентрации ЛОМ
При вычислении интеграла (2) с каждым изменением времени £ на значение ДЬ число итераций увеличивается в арифметической прогрессии. Для уменьшения числа итераций необходимо осуществить ряд процедур, которые приведут в конечном итоге к сокращению некоторых членов в вычисляемом интеграле. В программе для численных расчётов разработаны и применяются две процедуры: SetAnalysisIntegral и AnalysisIntegral. Процедура SetAnalysisIntegral позволяет включать или выключать анализ интеграла (2). Если время £ кратно параметру кь, то анализ интеграла включается: в памяти выделяется пространство для динамического массива С, размер которого определяется как Ь/Ьзьер + 1; при вычислении концентрации AGM значение концентрации вещества на каждой итерации сохраняется в массив С. Параметр кь определяет промежуток времени через который необходимо проводить анализ интеграла. Например, если этот интервал составляет 500 с, то параметр кь должен равняться 500. Применение процедуры AnalysisIntegral позволяет проводить анализ массива С и устанавливать новое
значение переменной Ъ0 • В обычном случае при вычислении интеграла (2) предел интегрирования находится на отрезке от 0 до Ъ. При сокращении числа итераций нижний предел интегрирования принимает значение Ъ0 • Значение параметра Ъ0 определяется следующим образом. Сначала вычисляется сумма 5 всех элементов массива С. Далее следует вычисление суммы 5а элементов массива С с наложением условия: элементы суммируются до тех пор, пока не будет выполнено условие
Рь < - ms • 5а| < Ре. (7)
Параметры рь, тв и ре вычисляются подбором, описанным ниже в разделе 4, 5 — сумма всех элементов массива С, тогда как 5а является суммой первых элементов массива С, значения которых не удовлетворяют условию (7). После выполнения указанного условия запоминается номер итерации в переменной пцет, далее получаем значение Ъ0:
Ъ0 = пИет • ЪвЬер. (8)
Чтобы во время численных экспериментов было возможным применение метода сокращения числа итераций, необходимо найти значения параметров рь, тв и ре.
4. Подбор параметров для сокращения числа итераций при вычислении концентрации ЛСМ
Метод сокращения числа итераций позволяет удалять значения первых итераций, суммой которых можно пренебречь. Для определения значения параметров рь, тв и ре проведено численное моделирование системы из двух нейронов (рис. 1). У одного из них происходит рост аксона. Второй нейрон является активным, поэтому притягивает аксон первого нейрона. Во время роста аксона при вычислении градиента концентрации АСЫ значение концентрации АСЫ на каждой итерации сохранялось в файл. Сохранения значений числа итераций происходило через каждый 500 с. Рост аксона до цели длился 61 500 с. Всего было создано 124 файла со значениями концентраций, вычисленными на каждой итерации.
Рис. 1. Три этапа эволюции нейронной сети из двух нейронов: а — 27400 с, б — 53 900 с, в 61500 с
а
в
I
б
II
б
Рис. 2. Графики со значениями концентрации ЛОМ на каждой итерации: а — 1000 с, б 12 000 с, в — 61 500 с (II — в увеличенном масштабе)
а
в
а
в
Значения концентрации AGM на каждой итерации для различных моментов времени £ представлены в виде графиков (рис. 2), из которых хорошо видно, что концентрации на большей части итераций либо равны нулю, либо столь близко приближены к нулевому, что ими можно пренебречь. Для момента времени £ = 1000 с к числу итераций, которые можно не учитывать, мы отнесли те, которые лежат в диапазоне от 0 до 300 (см. рис. 2, I, а) ив увеличенном масштабе — от 0 до 200 (рис. 2, II, а) При £ = 12 000 с число итераций, которые можно не учитывать, лежат в диапазоне от 0 до 5600 (рис. 2Д, б) или в увеличенном масштабе — от 0 до 3500 (рис. 2, II, б). При £ = 61 500 с можно не учитывать число итераций от 0 до 32 000 (рис. 2, I, в) или в увеличенном масштабе — от 0 до 24 000 (рис. 2, II, в).
Значения параметров ш8 и рь должны быть такими, чтобы выражение
|£ ёП - £ ёПм)| (9)
было наименьшим, а число итераций Мцег наибольшим. Здесь ^ ёп — сумма радиус-векторов сегментов аксона, ^ ёП(м) — сумма радиус-векторов сегментов аксона с учётом удалённых итераций. Значение ёп определяется дифференциальным уравнением (4). Индекс п означает номер конуса роста, г — номер нейрона.
Нахождение значений параметров с помощью разработанной компьютерной программы представляет собой следующий алгоритм. Считываются файлы со значениями концентраций AGM для каждой итерации. При помощи уравнения (4) вычисляются ёП и ёП(М) при различных значениях параметров шз и ре. Параметр ре принимает все значения концентраций итераций от 1000 до 50 000, которые кратны 1000. Параметр шз принимает четыре значения: 10, 100, 1000 и 10 000. Остальные используемые параметры приводятся в табл. 1.
Таблица 1. Значения параметров численных экспериментов
Параметр Значение Единица измерения
Количество ЛОМ, испускаемое нейроном 10-5 нМ/с
за единицу времени а [14-16]
Коэффициент диффузии Б2 [14] 6•10-5 см2 /с
Коэффициент, описывающий чувствитель- 4•10-6 см2 /нМс
ность и подвижность аксонов, А [15]
Коэффициент деградации к 10-3 с-1
Д£ 100 с
Шаг интегрирования £8^ер 1 с
Коэффициент к 500
Далее для всех наборов параметров ш8 и ре вычисляется выражение (9), результат сохраняется в массив С1. Затем производится поиск значений параметров ш8 и ре, для которых число итераций максимально и разность (9) минимальна. Параметр рь принимает значение 1 для всех значений ш8 и ре.
Результаты первого численного эксперимента для ш8 = 10 000 и ре = 33 000 ^ 50 000 представлены на рис. 3. На горизонтальной оси от 0 до 200 отложены номера элементов массива С1, каждый из которых определяется результатами численного анализа при различных значениях параметров ш8 и ре. Например, число 50 означает номер элемента, в котором содержатся результаты численного анализа с параметрами ш8 = 10 и ре = 1000. На вертикальной оси от 0 до 0.46 отложены значения разности (9). На горизонтальной оси от 106 до 3.2 • 106 представлены значения числа итераций для различных значений параметров ш8 и ре.
Во втором эксперименте ш8 присваиваем значение 10 000, поскольку из первого эксперимента при таком значении данного параметра был получен необходимый результат. Параметр ре принимает значения от 1000 до 100 000, кратные 1000. Полученный результат представлен на рис. 4.
Рис. 3. Трёхмерный график значений числа итераций разности (9) при различных значениях параметров ре и т8. Численный эксперимент выполнен при т8 = 10 000, ре принимает значения от 10 000 до 50000. Горизонтальная ось от 0 до 200 означает номер элемента массива С1
Номер элемента массива CI О
Рис. 4. Трёхмерный график значений числа итераций разности (9) с такими же, что и на рис. 3, значениями параметров ре и ш8, находящимися при этом в других диапазонах и имеющими меньшее число наборов. ш8 = 10000, ре принимает значения от 1000 до 100000, кратные 1000
Из первого и второго численных экспериментов максимальное число итераций было получено при pe = 33 000, при этом разность (9) равна 0.000 045 590 087, число сокращённых итераций Niter = 3 591 247. В первом эксперименте минимальное значение разности (9) равно 0.000 044 846 805, при этом число сокращённых итераций составляет Niter = 3 591247.
5. Экспериментальная часть
Для подтверждения эффективности используемого метода оптимизации вычислений смоделированы нейронные сети с двумя нейронами. Моделирование производилось четыре раза с повтором каждого эксперимента 10 раз при различных значениях параметров ms и pe. Первый эксперимент проводился без применения метода сокращения числа итераций, в остальных трёх параметры ms и pe имели следующие значения:
1. pe = 5000, ms = 10 000;
2. pe = 33 000, ms = 10 000;
3. pe = 100 000, ms = 10 000.
В экспериментах измерялось время от начала вычислений при запуске программы до момента, когда аксон нейрона 0 достигал своей цели — нейрона 1 (см. рис. 1). Измерение времени осуществлялось с помощью функции omp_get_wtime(), которая относится к библиотеке OpenMP. Фиксировались моменты времени t\ при запуске программы и t2, когда аксон достигал своей цели. Время, необходимое для роста аксона до нейрона 1, t = t2 — t\. Эксперименты проводились на компьютере с процессором Intel Core i7-2600 (4 ядра) 3.4 ГГц, 8 Гб ОЗУ, 64-bit Windows 7. Полученные результаты приведены в табл. 2.
В эксперименте без применения метода сокращения числа итераций аксон достиг своей цели за 61 600 с. Среднее значение реального времени роста аксона равно 32.055 с. В остальных трёх случаях аксон достиг своей цели, так же как и в случае без применения метода сокращения числа итераций, за 61 500 с. Это значит, что разница между временем роста аксона с применением и без применения метода сокращения числа итераций незначительна. При измерении же реального времени роста аксона получены
Таблица 2. Результаты четырёх экспериментов с различными значениями ш3 и ре
Номер Без сокра- Ре = 5000, Ре = 33 000, Ре = 100 000,
п / п щения числа итераций ш3 = 10 000 ш3 = 10 000 ш3 = 10 000
1 32.52 28.18 23.64 25.91
2 32.12 29.08 23.53 28.11
3 32.45 27.84 23.59 28.50
4 31.51 26.68 23.66 27.91
5 32.07 28.02 23.48 28.15
6 32.43 26.93 23.33 28.60
7 31.39 27.22 23.73 27.47
8 31.72 27.49 23.47 27.45
9 32.19 28.13 23.70 28.59
10 32.15 29.05 23.71 27.35
Среднее
значение 32.055 27.862 23.584 27.804
различные значения. Так, наименьшее из них было при ре = 33 000 и ш3 = 10 000, и среднее значение равнялось 23.584 с. Это на 26.42 % меньше, чем без применения метода сокращения числа итераций. При ре = 5000 и 100 000 среднее значение времени равно 27.862 и 27.804 с, что соответственно на 13.08 и 13.26 % меньше, чем без применения метода сокращения числа итераций. Поскольку при ре = 33 000 и ш3 = 10 000 число сокращённых итераций максимально, а значение реального времени роста аксона минимально, то эти значения ре и ш3 были выбраны для применения в методе сокращения числа итераций при вычислении концентрации АОМ.
Таким образом, представленный метод сокращения числа итераций при вычислении градиента концентрации АОМ позволяет сократить время численного анализа на четверть. Данный метод можно использовать в задачах, где необходимо вычислить концентрацию вещества или температуру при помощи уравнения диффузии (теплопроводности) с источником, зависящим от времени.
Заключение
Описанный метод был разработан для модели эволюции нейронных сетей. Численное моделирование занимает достаточно длительное время даже в простом случае моделирования системы из восьми нейронов. Реальные биологические нейронные сети состоят из сотен и тысяч нейронов. Если попытаться численно реализовать модель, состоящую из нескольких сотен нейронов, на обычном компьютере, то время численного эксперимента может растянуться до бесконечности. Предложенная модель реализована с применением метода параллельных вычислений, что эффективно, когда численный анализ реализуется на компьютерах с большим количеством вычислительных ядер. В работе [13] моделировались нейронные сети, состоящие из 12 и 18 нейронов. В численных экспериментах использовались методы параллельных вычислений. Это хотя и позволило сократить время проведения численного эксперимента, но незначительно, поскольку численные анализы проводились на компьютере с четырьмя вычислительными ядрами. Метод сокращения числа итераций при вычислении концентрации АОМ даёт возмож-
ность сократить время численного экспериментах без увеличения вычислительных ресурсов. Применение представленного метода совместно с параллельными вычислениями позволит сократить время численных экспериментов на многоядерных компьютерах.
Список литературы / References
[1] Chklovskii, D.B., Mel, B.W., Svoboda, K. Cortical rewiring and information storage // Nature. 2004. Vol. 14, No. 431. P. 782-788.
[2] Keynes, R., Cook, G.M.W. Axon guidance molecules // Cell. 1995. Vol.83, No.2. P.161-169.
[3] Dickson, B.J. Molecular mechanisms of axon guidance // Science. 2002. Vol. 298, No. 5600. P. 1959-1964.
[4] Tessier-Lavigne, M., Goodman, C.S. The molecular biology of axon guidance // Ibid. 1996. Vol. 274, No. 5290. P. 1123-1133.
[5] Szebenyi, G., Callaway, J.L., Dent, E.W., Kalil, K. Interstitial branches develop from active regions of the axon demarcated by the primary growth cone during pausing behaviors // J. Neurosci. 1998. Vol. 18, No. 19. P. 7930-7940.
[6] Dent, E.W., Callaway, J.L., Szebenyi, G. Reorganization and movement of microtubules in axonal growth cones and developing interstitial branches // Ibid. 1999. Vol. 19, No. 20. P. 8894-8908.
[7] Kalil, K., Szebenyi, G., Dent, E.W. Common mechanisms underlying growth cone guidance and axon branching // Ibid. 2000. Vol. 44, No. 2. P. 145-158.
[8] Dent, E.W., Tang, F., Kalil, K. Axon guidance by growth cones and branches: common cytoskeletal and signaling mechanisms // Neuroscientist. 2003. Vol. 9, No. 5. P. 343-353.
[9] Balkowiec, A., Katz, D.M. Activity-dependent release of endogenous brain-derived neurotrophic factor from primary sensory neurons detected by ELISA in situ //J. Neurosci. 2000. Vol. 20, No. 19. P. 7417-7423.
[10] Henley, J., Poo, M.M. Guiding neuronal growth cones using Ca2+ signals // Trends Cell Biol. 2004. Vol. 14, No. 6. P. 320-330.
[11] Gafarov, F., Khusnutdinov, N., Galimyanov, F. Morpholess neurons compromise the development of cortical connectivity // J. Integr. Neurosci. 2009. Vol. 8, No. 1. P. 35-48.
[12] Сулейманов Я.А. Математическое моделирование структурной пластичности в нейронных сетях // Вестник ТГГПУ. 2011. Т. 25, № 3. С. 37-41.
Suleymanov, Ya.A. Mathematical modeling of structural plasticity in neural networks // Bulletin of the TSHPU. 2011. Vol. 25, No. 3. P. 37-41 (in Russian).
[13] Suleymanov, Y., Gafarov, F., Khusnutdinov, N. Modeling of interstitial branching of axonal networks //J. Integr. Neurosci. 2013. Vol. 12, No. 1. P. 103-116.
[14] Goodhill, G.J. Diffusion in axon guidance // Europ. J. Neurosci. 1997. Vol. 9, No. 7. P. 1414-1421.
[15] Rosoff, W.J., Urbach, J.S., Esrick, M.A. A new chemotaxis assay shows the extreme sensitivity of axons to molecular gradients // Nat. Neurosci. 2004. Vol. 7, No. 6. P. 678-682.
[16] Gundersen, R.W., Barrett, J.N. Characterization of the turning response of dorsal root neurites toward nerve growth factor //J. Cell Biol. 1980. Vol. 87, No. 3. P. 546-554.
Поступила в 'редакцию 2 апреля 2014 г. с доработки — 9 декабря 2014 г.
62 H. A. CynemuaHOB, 0. M. ra^apoB, H. P. XycHyTguHOB, H. A. EuenbSHOBa
Optimization of calculation of matter concentration in the diffusion equation with a time-dependent source
suleymanov, yanis a.1'*, gafarov, fail m.1, khusnutdinov, nail r.1, Yemelyanova, Natalya a.2
1 Kazan Federal University, Kazan, 420008, Russia
2Kazan branch of Moscow State University of Railway Engineering, Kazan, 420008, Russia * Corresponding author: Suleymanov, Yanis A., e-mail: [email protected]
Purpose. This paper addresses an acceleration of numerical experiments in the study of the evolution of neural networks based on structural plasticity. Structural plasticity is related to the anatomical structure of neurons and connections between them. Modeling of the growth for neuron axons is a complex and lengthy process, requiring the use of cumbersome calculations. Axon growth is controlled by diffusion substance AGM, which is released by neurons and the axon tip is perceived it as a guidance signal. Since the numerical solution of the diffusion equation with time-dependent source is a complex process and resource-intensive process, it is necessary to optimize the program.
Methodology. In this paper we have developed a method which reduces the time of numerical experiments for the modeling of the evolution of neural networks. This method allows us to reduce the number of iterations in the calculation of the concentration of the AGM substance. The developed model is implemented using the method of parallel computing. The method is applied for the numerical solution of the diffusion equation with a time-dependent source.
Findings. During the numerical analysis we simulated a neural network, which consists of a set of neurons arranged in the three-dimensional space. The model also contains the mechanism of axon branching. Axon branching occurs through the formation of a new branch in the axis of the axon. The growth and branching of axons and dendrites leads to anatomical changes of the neural networks (structural plasticity).
Originality\value Application of the proposed method for reduction of the number of iterations in the process for calculation of the concentration of AGM significantly reduces the time of numerical experiments without increasing computational resources. This method can be used in applications which use the diffusion equation (conduction) with a time-dependent source.
Keywords: diffusion equation, numerical optimization, neural networks.
Received 2 April 2014 Received in revised form 9 December 2014
© ICT SB RAS, 2015
