CC BY
14
УДК 543.4:544.2
Я. А. Сулейманов, Ф. М. Гафаров, Н. А. Емельянова, В. Н. Александров
ОПТИМИЗАЦИЯ ВЫЧИСЛЕНИЯ ГРАДИЕНТА КОНЦЕНТРАЦИИ ВЕЩЕСТВА В УРАВНЕНИИ ДИФФУЗИИ С ИСТОЧНИКОМ, ЗАВИСЯЩИМ ОТ ВРЕМЕНИ
Ключевые слова: уравнение диффузии, оптимизация вычислений, нейронные сети.
Разработан метод сокращения числа итераций вычисления градиента концентрации вещества. Метод применен для численного решения уравнения диффузии с источником, зависящим от времени. Предложенный подход представляет собой инструмент, с помощью которого можно сократить время численного анализа на 13-26%.
Keywords: diffusion equation, numerical optimization, neural networks.
The method of iteration reduction for matter gradient concentration calculation is developed. The method is applied for numerical solution of diffusion equation with depending on time source. The time of numerical analysis is decreased on 13-26 % by taking into account the suggested approach.
Введение
Численное моделирование сложных систем состояния, которых в каждый момент времени зависят от множества параметров, как и квантово-химических расчетов [1,2], представляет собой достаточно сложный вычислительный процесс. Для оптимизации таких задач существуют различные методы: оптимизация алгоритма программного кода, перевод алгоритма из последовательной в параллельную парадигму вычислений, упрощение самой задачи.
В представленной работе приводится метод, позволяющий сократить время численных экспериментов при моделировании эволюции нейронных сетей. В качестве эволюции нейронных сетей рассматривается структурная пластичность в нейронных сетях. Структурная пластичность связана с анатомической структурой нейронов и определяется связями между нейронами [3]. Рост и ветвление аксонов и дендритов приводят к анатомическому изменению нейронных сетей. Моделирование роста аксона каждого нейрона представляет собой сложный и длительный процесс, требующий использования громоздких вычислений. Ключевым явлением для роста аксона является диффузия вещества AGM (axon growth molecules) [4-6], которое испускается нейронами и воспринимается кончиком аксона как управляющий сигнал. Для того чтобы в заданной точке определить значение концентрации AGM, используем стандартное уравнение диффузии. Поскольку численное решение уравнения диффузии с источником, зависящим от времени, представляет собой сложный процесс, то необходимость оптимизации работы программы появляется еще на стадии компьютерной реализации математической модели. Для проведения необходимой оптимизации разработан метод, который позволяет анализировать и сокращать число итераций при вычислении концентрации вещества AGM.
Объект исследования
Объектом исследования в работе является модель, описывающая динамическое образование межнейронных связей за счет вещества, выделяемого самими нейронами с учётом ветвления аксонов. В
модели используется стандартное уравнение диффузии
dr
-j- - D2r, - kr = Jl(rj,t),
(1)
где й - коэффициент диффузии, Ji(rj,t) - функция источника, зависящего от времени t. Источником в модели является нейрон. Решением уравнения (1) является выражение для концентрации
ri(ri - r
,t) = ajdtkGd(rj -r„t - tk)j(tk) ,
(2)
где j i(tk) - активность нейрона, а - количество вещества AGM, выделяющегося из нейрона за единицу времени, - г - tk) - функция Грина, которая выглядит следующим образом:
- 1 Г2 вс1(г^) =-^-ттуexp(-М---).
(4tD2 п )а/2 4tD2
В модель включен механизм ветвления аксона. Ветвление аксона происходит путем образования новой ветви из оси аксона [7-10]. Для определения точки ветвления на оси аксона и его начального направления роста необходимо иметь данные о значении концентрации AGM в одной из множества точек на оси аксона и в точках на сфере вокруг точки ветвления.
Движение аксона описывается дифференциальным уравнением
dt
•=¡Fa,)^™^, - rk,t). (3)
При проведении численного анализа моделируются нейронные сети, которые состоят из множества нейронов, расположенных в трехмерном пространстве. Нейрон обладает активностью, которая влияет на рост аксона [11,12]. При высоких значениях активности у нейрона аксон не растет и, наоборот, при низких значениях активности нейрона происходит рост аксона. Таким образом, система из множества нейронов со временем образует связи между собой. Аксон нейрона растет в направлении наибольшего градиента концентрации AGM. Ре-
k=1
зультаты численных экспериментов были опубликованы в работах [13-15].
Теоретическое описание метода сокращения итераций
В выражении (2) с каждым увеличением времени ( на значение Д, число итераций при вычислении концентрации ЛвМ увеличивается на значение Дt / tstep , где tstep - шаг интегрирования. Значение концентрации ЛвМ вычисляется отдельно для каждой ветви аксона. При проведении численного анализа число ветвей увеличивается с течением времени из-за процесса ветвления аксона. Увеличение числа ветвей приводит к увеличению числа итераций вычисления концентрации ЛвМ. Таким образом, при вычислении концентрации ЛвМ на конце ветви нейронной сети, состоящей из N нейронов, каждый из которых на аксоне имеет Nbr ветвей, число итераций равно
NNJ
N e =■
step
Данное выражение показывает число итераций, необходимое для вычисления либо значения концентрации ЛвМ, либо градиента концентрации ЛвМ в момент времени t. Необходимо учитывать, что рост аксона может длиться днями. Например, для нейронной сети, состоящей из 8 нейронов с одним кончиком аксона, которая развивается 2 дня, то есть, t =172800 с, значение параметра tstep =1 с, число итераций для вычисления концентрации ЛвМ равно 1382400 итераций в момент времени t. Здесь время t измеряется в секундах, а число Nbr =1. Эволюция нейронной сети может длиться больше двух дней, что, соответственно, потребует большего количества времени.
Кроме движения аксона в модели реализован процесс ветвления, который связан с вычислением концентрации ЛвМ. В отличие от вычисления концентрации ЛвМ во время роста аксона, где расчет производится в одной точке, в процессе ветвления для определения направления роста новой ветви вычисления ведутся на множестве точек, расположенных на сфере
Niter = ■
1 step
где N^1 - число точек по широте, N|0t - число точек по долготе, Nbr - число ветвей, на которых происходит ветвление.
Сокращение количества итераций при вычислении концентрации ДОМ При вычислении интеграла (2) с каждым изменением времени t на значение Дt число итераций увеличивается в арифметической прогрессии. Для уменьшения числа итераций необходимо провести ряд процедур, которые приведут в конечном итоге к сокращению некоторых членов в вычисляемом интеграле. В программе для численных расчетов разработаны и применяются две процедуры: SetAnalysisIntegral и Лпа1у8181Ш^га1. Процедура SetAna1ysisIntegra1 позволяет включать или выклю-
чать анализ интеграла (2). Если время t кратно параметру к,, то анализ интеграла включается: в памяти выделяется пространство для динамического массива С, размер которого определяется как t / +1; при вычислении концентрации ЛвМ,
значение каждой итерации сохраняется в массив С . Применение процедуры Ana1ysisIntegra1 позволяет проводить анализ массива С и устанавливать новое значение переменной ,0 . В обычном случае, при вычислении интеграла (2), предел интегрирования находится на отрезке от 0 до t. При сокращении итераций нижний предел интегрирования изменяется и принимает значение ,0 . Значение параметра ,0 определяется следующим образом. Сначала вычисляется сумма 5 всех элементов массива С. Далее следует вычисление суммы 5а элементов массива С с наложением условия: элементы суммируются до тех пор, пока не будет выполнено условие:
Pb S - ms ■ Sa\ ^ Pe ■
(4)
После выполнения условия запоминается номер итерации в переменной п е , далее получаем значе-
ние t0 :
'о = niter ' ^step ■
(5)
Для того чтобы во время численных экспериментов было возможным применение метода сокращения итераций, необходимо найти значения параметров
Рь , т И Ре .
Подбор параметров для сокращения итераций при вычислении концентрации ДОМ
Метод сокращения итераций позволяет удалять значения первых итераций, суммой которых можно пренебречь.
1. . 0;52 0. 1. 0.
0.
0,15
0.
0. 0. 0.
0.
, 0,52 1. 0.
0,15
0. 0. 0.
Рис. 1 - Два этапа эволюции нейронной сети из двух нейронов: справа - 53900 с, слева - 61600 с
Для определения значения параметров рь, т8 и ре, проведено численное моделирование системы из двух нейронов (рис. 1). У одного из двух нейронов происходит процесс роста аксона. Второй нейрон является активным, поэтому притягивает аксон первого нейрона. Во время роста аксона при вычислении градиента концентрации ЛвМ, значение каждой итерации сохранялось в файл. Сохранения значений итераций происходило в течение 500 с, затем создавался новый файл. Рост аксона до цели
длился 61600 с. Всего было создано 124 файла со значениями каждой итерации.
Рис. 2 - На рисунке изображен график значений каждой итерации при вычислении концентрации ЛОМ в момент времени t =1000 с
Значения каждой итерации для различных моментов времени t представлены в виде графиков на рисунках 2 и 3. На графиках хорошо видно, что значения большей части итераций равно нулю, либо достаточно близко к нулевому так, что ими можно пренебречь. Для момента времени t =1000 с к числу итераций, которые можно не учитывать, мы отнесли те, которые лежат в диапазоне от 0 до 300 (рис. 2). Для значения t =61500 с число итераций, которые можно не учитывать, лежит в диапазоне от 0 до 24000 (рис. 3).
Рис. 3 - График значений каждой итерации при вычислении концентрации ЛОМ в момент времени t=61500 с
Значения параметров должны быть такими, чтобы выражение дП дП(М)\ было наименьшим, а число итераций N¡ter наибольшим. Значение ^ дП является суммой радиус-векторов сегментов аксона, ^ 9П(М) - сумма радиус-векторов сегментов аксона с учётом удаленных итераций. Значение дП определяется дифференциальным уравнением (3).
Нахождение значений параметров с помощью разработанной компьютерной программы представляет собой следующий алгоритм. Считы-ваются файлы со значениями каждой итерации. При
помощи уравнения (3) вычисляется значение дП и значение дП(М). Значения этих переменных вычисляются при различных значениях параметров т8 и РС). Параметр рС) принимает все значения итераций от 1000 до 50000, которые кратны 1000. Параметр принимает четыре значения: 10, 100, 1000 и 10000. Остальные используемые параметры приводятся в таблице 1.
Таблица 1 - Значения параметров численных экспериментов
Параметры Значение Единица измерения
Количество ЛвМ, испускаемое нейроном за единицу времени а [16, 17, 18] 10-5 нМ/с
Коэффициент диффузии й2 [16] 6-10-5 см2/с
Коэффициент, описывающий чувствительность и подвижность аксонов А [17] 4-10-6 см5/нМс
Коэффициент деградации к 10-3 с-1
Д, 100 с
Шаг интегрирования tstep 1 с
Далее для всех наборов параметров и вычисляем выражение дП дП(М)\, результат
сохраняем в массив С1. Затем производится поиск значений параметров т8 и рс,, для которых число
итераций максимально и разность дП - ^ д"(М^ минимальна. Параметр рь принимает значение 1 для всех значений параметров т8 и ре.
Результаты первого численного эксперимента для параметров тз =10000 и = 33000^50000 представлены на рисунке 4. На горизонтальной оси со значениями от 0 до 200 отложены номера элементов массива С1, каждый из которых определяется результатами численного анализа при различных значениях параметров т8 и рС). Например, число 50 обозначает номер элемента, в котором содержатся результаты численного анализа с параметрами тз =10 и рС! =1000. На вертикальной оси со значениями от 0 до 0,46 отложены значения разности дП дП(М)\. На горизонтальной оси со
значениями от 106 до 3,2-106 отложены значения числа итераций для различных значений параметров
ms и Рe.
Во втором численном эксперименте параметру т5 присваиваем значение 10000, поскольку из первого численного эксперимента при таком значении был получен необходимый результат. Параметр ре принимает значения от 1000 до 100000, кратные 1000.
1 х 10Б
i7-2600 (4 ядра) 3,4 GHz, 8 Gb ОЗУ, 64-bit Windows 7. Результаты экспериментов приведены в табл. 2.
Таблица 2 - Результаты 4-х экспериментов с различными параметрами ms и pe
150
100
50 о
Рис. 4 - Трехмерный график значения числа итераций разности д? 9П(М)\ при различных значениях параметров ms и pe
Из первого и второго численного экспериментов максимальное число итераций было получено при pe =33000. При таком pe разность
9? 9П(м\ равна 0,000045590087, число сокращенных итераций Niter =3591247. В первом численном эксперименте минимальное значение разности
д? дп(и)\ равно 0,000044846805, при этом число сокращенных итераций Niter =3591247.
Экспериментальная часть
Для подтверждения эффективности используемого метода оптимизации вычислений смоделированы нейронные сети с двумя нейронами. Моделирование производилось четыре раза при различных значениях параметров ms и pe . Первый эксперимент проводился без применения метода сокращения итераций. Остальные три эксперимента проводились со следующими значениями параметров ms и pe :
2. pe=5000, ms=10000;
3. pe=33000, ms=10000;
4. pe=100000, ms=10000.
Каждый эксперимент проводился 10 раз. В экспериментах измерялось время от начала вычислений при запуске программы до момента, когда аксон нейрона 0 достигал своей цели - нейрона 1 (рис. 1). Измерение времени осуществлялось с помощью функции omp_get_wtime(), которая относится к библиотеке OpenMP. Фиксировалось значение времени t1 при запуске программы и значение времени t2 , когда аксон достигал своей цели. Время, необходимое для роста аксона до нейрона 1, t = t2 -11. Эксперименты проводились на компьютере с процессором Intel Core
№ экспе ри-мента Без сокращения итераций pe = 5000, ms = 10000 pe = 33000, ms = 10000 pe = 100000, ms = 10000
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
нее
зна-
чение
В эксперименте без применения метода сокращения итераций аксон достиг своей цели за 61600 с. Среднее значение реального времени роста аксона до своей цели равно 32,055 с. В остальных трех экспериментах аксон достиг своей цели, так же как и в случае без применения метода сокращения итераций за 61600 с. Это значит, что разница значения времени роста аксона с применением и без применения метода сокращения итераций незначительна. При измерении же реального времени роста аксона получились различные значения времени. Так, наименьшее значение получилось, когда ре = 33000 и тз = 10000, и среднее значение равно 23,584 с. Это на 26,42 % меньше, чем без применения метода сокращения итераций. Когда параметр ре принимает значение 5000 и 100000, среднее значение времени равно 27,862 с и 27,804 с, что на 13,08 % и 13,26 % соответственно меньше, чем без применения метода сокращения итераций. Поскольку при ре = 33000 и тз = 10000 число сокращенных итераций
максимально, а значение реального времени роста аксона минимально, то эти значения были выбраны для применения в методе сокращения итераций при вычислении концентрации ЛвМ.
Вывод
Разработан метод сокращения итераций при вычислении градиента концентрации ЛвМ. Метод сокращения итераций позволяет сократить время численного анализа на четверть. Данный метод можно использовать в задачах, где необходимо вычислить концентрацию вещества или температуру при помощи уравнения диффузии (теплопроводности) с источником, зависящим от времени.
Заключение
Метод, описанный в работе, был разработан для модели эволюции нейронных сетей. Численное моделирование занимает достаточно долгое время даже в простом случае моделирования системы из 8 нейронов. Реальные биологические нейронные сети состоят из сотен и тысяч нейронов. Если попытаться численно реализовать модель, состоящую из нескольких сотен нейронов, на обычном компьютере, то время численного эксперимента может растянуться до бесконечности. Наша модель реализована с применением метода параллельных вычислений. Применение параллельных вычислений эффективно, когда численный анализ реализуется на компьютерах с большим количеством вычислительных ядер. В нашей работе [15] моделировались нейронные сети, состоящие из 12 и 18 нейронов. В численных экспериментах использовались методы параллельных вычислений. Это позволило сократить время проведения численного эксперимента, но не значительно, поскольку численные анализы проводились на компьютере с 4 вычислительными ядрами. Применение метода сокращения итераций при вычислении концентрации ЛвМ позволило сократить время численного экспериментах без увеличения вычислительных ресурсов. Применение предложенного в работе метода совместно с параллельными вычислениями позволит сократить время численных экспериментов на многоядерных компьютерах.
Литература
1. Т.Н. Гришаева, А.Н Маслий, Вестник Казан. технол. ун-та, 15, 12, 7-11 (2012).
2. А.Н. Маслий, Э.И. Мадиров, Вестник Казан. технол. ун-та, 16, 23, 12-18 (2013).
3. D.B. Chklovskii, B.W. Mel, K. Svoboda, Nature, 14, 431, 782-788 (2004).
4. R. Keynes, G.M.W. Cook, Cell, 83, 2, 161-169 (1995).
5. B.J. Dickson, Science, 298, 5600, 1959-1964 (2002).
6. M. Tessier-Lavigne, C.S. Goodman, Science, 274, 5290, 1123-1133 (1996).
7. G. Szebenyi, J.L. Callaway, E.W. Dent, K. Kalil, J Neurosci., 18, 19, 7930-7940 (1998).
8. E.W. Dent, J.L. Callaway, G. Szebenyi, J Neurosci., 19, 20, 8894-8908 (1999).
9. K. Kalil, G. Szebenyi, E.W. Dent, J Neurobiol, 44, 2, 145-158 (2000).
10. E.W. Dent, F. Tang, K. Kalil, Neuroscientist, 9, 5, 343353, (2003).
11. A. Balkowiec, D.M. Katz, J Neurosci., 20, 19, 74177423 (2000).
12. J. Henley, M.M. Poo, Trends Cell Biol., 14, 6, 320-330 (2004).
13. F. Gafarov, N. Khusnutdinov, F. Galimyanov, J Integr Neurosci., 8, 1, 35-48 (2009).
14. Я. А. Сулейманов, Вестник ТГГПУ, 25, 3, 37-41 (2011).
15. Y. Suleymanov, F. Gafarov, N. Khusnutdinov, J Integr Neurosci., 12, 1, 103-116 (2013).
16. G. J. Goodhill, Eur J Neurosci., 9, 7, 1414-1421 (1997).
17. W.J. Rosoff, J.S. Urbach, M.A. Esrick, Nat Neurosci., 7, 6, 678-682 (2004).
18. R. W. Gundersen, J.N. Barrett, J Cell Biol., 87, 3, 546554 (1980).
© Я. А. Сулейманов - асп. каф. информационных систем КФУ, yaniscom@mai1.ru; Ф. М. Гафаров - к.ф.м.н, доцент той же кафедры, fgafarov@yandex.ru; Н. А. Емельянова - к.ф.м.н, доцент Казанского филиал Московского госуд. ун-та путей сообщения (МИИТ), pr170@mai1.ru; В. Н. Александров - к.т.н, доц. каф. химии и технологии высокомолекулярных соединений КНИТУ, 1abgor@kstu.ru.
