К имитационному моделированию нервного импульса

Сер. 10. 2011. Вып. 4

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

УДК 531:518:577 В. С. Новоселов

К ИМИТАЦИОННОМУ МОДЕЛИРОВАНИЮ НЕРВНОГО ИМПУЛЬСА

1. Информационные нервные импульсы зарождаются в шейках сенсорных нейронов, передаются с помощью химического синапса на сеть ассоциативных нейронов и поступают в центральную нервную систему [1, 2]. При этом сохраняются их профиль и величина амплитуды. Команды центральной нервной системы на мышцы и различные внутренние органы передаются также с помощью подобных импульсов. Для управления целенаправленными движениями человека [3] важное значение имеет создание практически приемлемой теории передачи информации в реальных нейронных сетях. В ряде исследований [4-6] для построения макроскопической модели распространения нервных импульсов по нейронному волокну предполагается использовать длинноволновое уравнение Кортевега-де Фриса (KdV). В работах [7-14] создана глубокая теория уединенных волн (солитонов) этого уравнения, получены точные аналитические решения как для изолированных волн, так и для серии волн. Однако все еще отсутствует чисто физическое обоснование применимости уравнения KdV к описанию динамики нервного импульса в том виде, как это сделано, например, в книге [14] для уединенной волны на мелкой воде, акустической волны в коллоидной плазме и зарождающейся магнитной волны в намагниченной плазме. В настоящей статье предлагается физическое доказательство применимости уравнения KdV для описания распространения нервного импульса в немиелинированном аксоне.

При построении динамики нервного импульса находят широкое использование различные варианты уравнений в частных производных типа реакция-диффузия [15-19] или обобщенного уравнения Фоккера-Планка [5]. Рассматриваются модели с одной пространственной координатой x и осуществляется переход к волновой фазовой переменной £ = x — vt, где v = const>0 - скорость распространения волны, t - время. Уравнения в частных производных относительно x и t переходят в обыкновенные уравнения для переменной £. Есть существенное отличие между профилем информационной волны, который можно считать практически симметричным, и профилем нервного импульса в мышечной ткани сердца, который имеет короткий и крутой передний фронт, а также длительный и пологий задний фронт. Точных решений таких уравнений нет. Авторы приводят профили волн, полученные численно, но не указывают, при каких значениях параметров они получены. Однако автоволны в рассматриваемых уравнениях возникают при определенных трудно уловимых значениях параметров, и повторить вычисления авторов очень трудно. Для выявления структуры волны недостаточно

Новоселов Виктор Сергеевич — доктор физико-математических наук, профессор кафедры механики управляемого движения факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 190. Научные направления: аналитическая механика, динамика управляемых систем, биомеханика, вероятностные процессы. E-mail: novoselov@apmath.spbu.ru.

© В. С. Новоселов, 2011

используются кусочно-линейные модели, решения для которых можно получить в конечном виде. В настоящей статье описываются три кусочно-линейные модели. Первая модель является сильно упрощенной. Она сводится к исследованию линейного дифференциального уравнения второго порядка. Вторая модель имитирует уравнения модели ФитцХью-Нагумо [17, 19], третья - уравнения модели Ходжкина-Хаксли [16-19].

2. На основе анализа баланса энергии элемента нервного волокна длиной Дх можно предложить следующее физическое обоснование применения солитонного длинноволнового KdV-уравнения в динамике нервного импульса. Укажем особенности функционирования нервных клеток [3, 5, 16-19]. Постоянно работающие клеточные насосы засасывают ионы К+ и выталкивают ионы в пропорции 2:3, тем самым поддерживают на мембране нейрона равновесный электрический потенциал ир. При этом концентрация ионов внутри клетки оказывается существенно заниженной, а концентрация ионов К+ завышенной. Если внешнее возбуждение достигает значения порогового потенциала ип, то включается механизм самовозбуждения клетки (открываются натриевые каналы). Входящие на участке нейрона длиной Дх ионы повышают за время ДЬ электрическую энергию клетки на величину идг, которая расходуется на изменение и - энергии диффузионного потока частиц вдоль оси ж аксоплазмы и на и2 - энергию перезарядки конденсатора мембраны на участке Дх, т. е.

идг = и1 + и2.

Энергию входящих за ДЬ на участке Дх положительных зарядов Д^ будем вычислять по формуле из книги [20]

идг =

здесь иср - некоторая постоянная величина, которая будет рассматриваться как неопределенный параметр. Обозначаем через с удельную емкость мембраны нейрона, приходящуюся на 1 см2 ее поверхности. Эта величина порядка 1 мкф • см~2 (см. [16, 18]). Пусть с - емкость мембраны на единицу длины, т. е. с = 2пгс, где г - радиус нервного волокна, равный по работе [19] для немиелинированного аксона 3-5 мкм, а для волокна, покрытого миелином, - от 4 до 20 мкм. Емкость участка Дх мембраны сДх. Поскольку емкость конденсатора равна отношению его заряда к разности потенциалов на обкладках [21], то для силы тока входящих ионов можно принять

ди

I = г—Ах, дд = ш.

На основании формул (2) и (3) имеем

ди

~т‘

Оценим изменение и за время ДЬ. Основную роль при этом играют положительные ионы и К+, несущие единичный положительный заряд е, поэтому удельное на единицу длины аксоплазмы количество электричества будет

Qуд = еМЛ цБ:

где Мл - число Авогадро; ц - количество вещества в молях на единицу объема; Б -площадь поперечного сечения нервного волокна. При открытии натриевых каналов и быстром входе в клетку ионов диффузионный ток вдоль аксоплазмы изменяется

£/д

не на единичном расстоянии, а на какой-то длине I. Следовательно, надо рассматривать количество электричества

Q = Qуд I = еМлцБ1.

Отсюда получаем выражение силы тока диффузионного перемещения частиц

т 9 1\Г Я1дд

1' = Ж = ‘Мл'181-Ш

и мощности этого тока

дц

дЬ

Здесь Ш - разность потенциалов в аксоплазме. Считаем эту величину, как и I, положительной постоянной, хотя и неопределенной.

Пусть уединенная волна перемещается в положительную сторону оси х. Запишем изменение электрической энергии диффузионного потока частиц в сечении ж + Дх за промежуток времени ДЬ следующим образом:

( д д \

II1 = еШНаБ

Ввиду предполагаемого сохранения формы волны при ее перемещении вдоль оси х,

следует считать

ц

Согласно уравнению диффузии, например из книги [5],

д1 = В^

дЬ дх2’ 1 >

где Б - постоянный для конкретного нервного волокна коэффициент диффузии. Выражение (5) на основании (6) и (7) представим в виде

и, = -,от^л8

Весь заряд нервного волокна располагается на внутренней поверхности мембраны. Он будет создавать потенциал и(х, Ь), который можно определить по формуле для конденсатора единичной длины

си = ецМлБ.

Формулу (8) запишем так:

д2 д2

111 = -сВШ

Определим энергию перезарядки конденсатора. Для этого запишем выражение энергии конденсатора (см. [20]) на участке от сечения х до сечения х + Ах, имеющего емкость сАх, для двух моментов времени і и і + Аі

UK(t) = — cAx(u(x,t) — Up)2,

£/к(£ + Д£) = -сАх(и(х, t + Д£) — мр)2.

За время Дг передний фронт волны переместится на величину Дж и, ввиду сохранения профиля волны,

u(x + Ax, t + At) = u(x, t).

Поэтому имеем

1.

Энергия перезарядки выделенного участка конденсатора равна

UK(t) = -cAx(u(x + Ax, t + At) — up)2.

U2 = UK(t + At) — UK(t) = ——cAx [_(u(x + Ax, t + At) — up)2 — (u(x, t + At) — up)2] .

Полагая Дж = vAt, где v = const > 0 есть скорость распространения уединенной волны, и учитывая непрерывность мембранного потенциала, приходим к выражению

U2 = ~2^v [_

Подставим в формулу (1) представления (4), (9) и (10), поделим на с(иср — ир) и перейдем к пределу при Дж ^ 0 и At ^ 0. Получим уравнение KdV

д(и — ир) д , ч9 „ д3 . .

dt ^adx^U _Мр') + dx3^U~Up^ = ^

где приняты обозначения для положительных коэффициентов

1

2

Параметр х = (lW)-1 = const > 0 является характеристикой кривизны фронта волны. При уменьшении l или W понижается мощность оттока положительных ионов от выбранного сечения. Чем меньше отток ионов вдоль аксона, тем выше в данном сечении мембранный потенциал. Рассматриваем х как второй неопределенный параметр. Уравнение (11) имеет точное односолитонное решение [5-14]

-v

u — up = За 1vch 2

ill. . —p2 V 2

Амплитуда А = 3а-1у нервного импульса оказывается одинаковой как для аксонов нервных волокон, так и для мышечных тканей сердца [2, 5, 16, 17]. Она равна 110 мВ. С помощью формулы (12) получаем иср — ир = = 18.33 мВ. Заметим, что для сер-

дечной мышцы [17]: ир = —90 мВ, ип = —75 мВ, ип — ир = 15 мВ, для аксона [16]: ир = —70 мВ, ип ~ —50 мВ, ип — ир « 20 мВ. Поэтому рассматриваемый первый неопределенный параметр имеет величины ип — ир, и его можно приближенно принимать равным

этому значению. Средняя ширина солитона (13) или средняя ширина одного фронта волны по работам [5, 6] равна Ь = . Учитывая обозначение (12), запишем

4Б

X

Ъ2ь(иСр — ир)'

Поэтому второй неопределенный параметр или параметр крутизны фронта волны обратно пропорционален квадрату средней ширины фронта волны. Укажем известные числовые значения [17, 19]. Коэффициент диффузии Б для химически возбудимой среды - примерно 10-4см2(с)-1, для мышечной ткани сердца - 0.1 см2(с)-1. Скорость V распространения нервного импульса (в мс-1): для миелинированной нервной ткани -30-100, для немиелинированной - 15-20, для мышечной ткани сердца - 0.3-4.

Формула (13) дает представление уединенной волны симметричного профиля с одинаковыми коэффициентами кривизны как для переднего, так и для заднего фронта. Такое решение может применяться для моделирования распространения информационных импульсов в нервных волокнах. Сердечный импульс является несимметричным. Он имеет крутой передний фронт и размазанный задний фронт с отношением характерных времен для указанных фронтов порядка 1:300 [17]. Такой нервный импульс можно моделировать формулой (13), но с различными параметрами крутизны для переднего и заднего фронтов.

3. К общепринятому [15-19] уравнению в частных производных типа реакция-диффузия можно прийти на основе следующего рассмотрения. Исходя из основной формулы конденсатора [20], запишем

ди

ст=3-

Здесь введено ранее не используемое обозначение ] для плотности тока через 1 см2 поверхности мембраны. Кроме мембранного тока ] будет еще диффузионный ток ионов вдоль аксоплазмы. Если закрыть все поры и каналы мембраны, то продольный ток определяется уравнением диффузии

ди г)^и пд2и

дЬ дх2 ’ дЬ дх2 ’

где обозначено, что Б = Бс. Суммируя указанные токи, приходим к следующему уравнению для секундного изменения потенциала:

ди г,д2и

с~т=3+ д^2' ^ ^

По терминологии книги [5] его можно рассматривать как обобщенное уравнение Фоккера-Планка

ди д2и .

~т= д^+С

Оно выражает продольную диффузию с учетом перемещения частиц через мембрану.

Рассмотрим упрощенную математическую модель распространения возмущения в активной среде. Пусть V = и — ир, Уп = ип — ир. Запишем граничное условие V(±го) = 0. Допускаем, что процесс автогенерации нервного импульса начинается при

V ^ Уп и заканчивается при достижении наибольшего значения потенциала У^ах = А. Упрощенным аналогом уравнения (14) будет

dV - д2 V

с—— = D-7—— /гУ + F, к = const > 0.

dt дх2

Здесь F = h = const > 0 при автогенерации и F = 0 - в ее отсутствие. Второй член правой части выражает медленный ток поперечного дрейфа ионов через поверхность мембраны, третий - быстрый натриевый ток, вызванный открытием каналов для входа ионов Na+. Выполним переход к волновой фазовой переменной £ = х — vt, v = const > 0:

D V" + cvV' — kV + F = 0.

Здесь и далее штрихами обозначаем производные по £. Характеристическое уравнение однородного уравнения имеет корни

Aij2 = ±^cwZ)_1 (V1 + 4Dkc-2v~2 =F lj , Ai > 0, A2 <0, Ai < |A2|.

Пусть за точку отсчета принято положение, в котором V достигает наибольшего значения A. Здесь х = 0, t = 0, £ = 0. Обозначим через £п > 0 пороговое значение волновой фазовой координаты, для которой V(£п) = Vn. Можно задать последнее слагаемое правой части уравнения (15) в виде

0 при £ < 0,

h при 0 £ < £п

0 при £ > £п.

На каждом из трех участков, отвечающих формуле (16), общее решение уравнения (15) запишется так:

V = С1 еЛ1« + С2еЛ^ + F.

Для первого участка переднего фронта (£ > £п), в силу граничного условия V(ж) = 0, имеем

С11 =0, С21 = Ке-Л2«-, V = VпeЛ2

Отсюда V/І£=le+ = X2Vп < 0. В нулевой точке второго участка (0 ^ ^ £п) переднего

фронта V(0) = А, V/|^=0+ = 0. Поэтому получим

V = С12еЛ1« + С22еЛ^ + к,

С12 = —А2(А — к)(А1 — А2) 1, С22 = А1 (А — к)(А1 — А2) 1.

С помощью формул (18) и (19) запишем условие непрерывности потенциала в точке £п

Можно потребовать также выполнения непрерывности V/ в этой точке

А1

Уравнения (20) и (21) дают условия на выбор величин A, D, C, Vn, k, v, h, £п. При этом должно быть h > A, что отвечает отрицательности V' для переднего фронта волны. Пусть £п << 1, что характерно для некоторых мышечных клеток сердца [16]. Опуская в (20) члены порядка £^, получим

£n = y/2Dk-l(A-Vn)(h-A)-K Для заднего фронта волны, представленного одним участком £ < 0, имеем условия

F = 0, V (—ж) = 0, V (0) = A.

По формуле (17) находим

V = C23 еЛ1«, C23 = A.

При этом V'|?=0- = AiA > 0. Чтобы добиться обращения в нуль указанной производной в состоянии наибольшего возбуждения, следует повысить порядок дифференциального уравнения для волновой фазовой координаты.

4. Рассмотрим кусочно-линейную систему, отвечающую модели ФитцХью-Нагумо [17, 19] вида

ж-мм-

Здесь V - мембранный потенциал, f (V, Зи) - плотность всего мембранного тока, Зи -плотность медленного дрейфа частиц (ионов) через поры мембраны. Примем представления

f = F — k1V + ]и, g = —k2V, k1 = const > 0, k2 = const > 0,

где F - быстрый натриевый ток, определяемый формулой (16); kiV - ток утечки по закону Ома.

Перейдем к волновой фазовой переменной

D V'' + cvV' — kiV + Зи + F = 0, vf'^ — k2 V = 0.

Для получения эквивалентного дифференциального уравнения третьего порядка надо дифференцировать обобщенную функцию F = h0(£) — h0(£ — £п), где в(£) - функция Хевисайда. Используя обозначение работы [5], приходим к уравнению

D V''' + cvV'' — k1V' + k2v-1 V = —h6

Характеристическое уравнение однородной системы

D A3 + cvA2 — k1 A + k2v-1 = 0

по теореме Декарта имеет один отрицательный вещественный корень A3. Два других корня могут быть вещественными положительными или сопряженными комплексными. Будем исходить из первого случая и примем A1 > A2 > 0.

Как следует из (18), отрицательный вещественный корень обеспечит возрастание

V от 0 до Vn на этапе предпорогового возбуждения при уменьшении £ от ж до £п по формуле V = Vnexp[A3(£ — £п)]. Для исследования последующего второго участка

переднего фронта найдем по формуле свертки Бореля из работы [5] частное решение уравнения

5

V = 1 Ш(£ - г)/(г)3,т, /(т) = В-\-Н5(0 + Н5(£ — )),

о

где импульсная переходная функция равна Ш(£ — т) = в(£ — т)г(£ — т). Здесь г = С1вЛ15 + С2еХ25 + СэвЛз5 есть решение однородного уравнения, отвечающего уравнению (23) при условии г5=о = 0, г'5=о = 0, г" ^=о = 1. Для указанных постоянных С*, г = 1, 2, 3, получим значения

С\ = —(Ах — А2) 1 (Аэ — А2) 1, 1 ^ 2 ^ 3.

С помощью оператора тождественного отображения [5] находим

ээ

V = — НВ-1в

*=1 *=1

Запишем общее решение уравнения (23) в виде

V = С12вЛ15 + С22вЛ2 5 + С32вХз 5 + V.

Учтем граничные условия второго участка

з з

V5=о+ =$3 С*2 = А, V'|5=о+ = £ А*С*2 = 0,

*=1 *=1

3

V 15=5- = ]Т

*=1

3

У%-5- = Х)(С<2 + Н В-1С*)А*вЛ5 = АэУп.

*=1

Три из перечисленных соотношений позволяют определить постоянные С12, С22 и Сэ2. Подстановка полученных выражений в оставшееся соотношение дает условие на выбор величин А, В, с, Уп, &1, ^2, V, Н, £п.

Изменение потенциала на заднем фронте волны определяется формулой (25) при

V = 0. Принимаем граничные условия: V(—те) = 0, V(0) = А, V'(0) = 0. В результате имеем

С13 = —А1АА — А2) 1, С2Э = А2А(А1 — А2) 1, С33 = 0.

Если А1 и А2 будут сопряженными комплексными числами, то задний фронт заканчивается колеблющимся затухающим хвостом. Заметим, что в этом случае изложение решения для переднего фронта также должно быть иным, так как экстремумы не будут совпадать с амплитудами колебаний [21].

5. Для представления бегущего импульса в немиелинированном нервном волокне следует обратиться к рассмотрению уединенной волны модели Ходжкина-Хаксли. Основное уравнение этой модели [15-19] имеет вид (14), в котором плотность мембранного тока представлена суммой натриевого и калиевого токов, а ток утечки как

j = jNa + jK + je- На основании экспериментов с одиночными нервными импульсами немиелинированного аксона гигантского кальмара авторам модели удалось формализовать зависимость проводимостей быстрых натриевых и калиевых токов от мембранного потенциала. Прямая физическая интерпретация такой зависимости затруднительна [16]. При этом используемые в ней характеристики специфичны для каждой конкретной нервной системы.

Если заменить плотности натриевого и калиевого токов их средними значениями, а для плотности тока утечки принять представления п. 4, то можно прийти к кусочнолинейной системе уравнений (22), в которой

F = М(£) - hd(£ - £п) - h*d(£ - £*) + Ъ*в(£), h = const > 0, h* = const > 0,

£п > 0, £* < 0, V(£п) = V(£*) = К.

Эта модель содержит четыре линейных участка. Из них два для переднего фронта волны — допорогового возбуждения при £ > £п и быстрого входящего тока ионов Na+ с характеристикой h при 0+ ^ ^ £п, а также два участка для заднего фронта волны - быстрого входа через мембрану ионов K+ с характеристикой h* при £* ^ ^ 0-

и затухания импульса при закрытых натриевом и калиевом каналах при £ < £*. Передний фронт волны определяется выбором тех же параметров, как и в п. 4, с помощью формул (18), (24)-(26). При этом последнее соотношение формулы (26) дает условие на выбор величин A, D, с, Vп, к\, к2, v, h, £п.

Решение для первого участка заднего фронта можно получить, исходя из решения второго участка переднего фронта. Действительно, для отрицательных £ имеем уравнение

DV"' + cvV" - k1V' + k2v-1V = -h*6

Введем в рассмотрение вспомогательную волновую фазовую переменную п = £ - £* и вспомогательную величину £п = -£* > 0. Левая часть уравнения (27) сохранит тот же вид, если штрихами обозначить производные по п. Правая часть этого уравнения запишется так: -h*(n - £*) + h*6(n). Уравнение (27) примет вид уравнения (23) при замене £ на п, h на -h*, £п на £*. Получим решение (24), (25), в котором надо выполнить ту же замену. Запишем граничные условия

3

V |,=о+ = V |?=е+ = £ Ci3 = К,

i=1

3

V In=?r = V |5=о- = ]T(Ci3 + h* D-1Ci )в-^' = A,

i=1

3

V'Iri=&- = V'|^=o- ^(Ci3 - h*D-1Ci)\ie-x^ = 0, i=1

которые при выборе £* однозначно определяют постоянные Ci3, i =1, 2, 3. Для производной в точке контакта со вторым участком заднего фронта имеем

V%=„+ = V 15=е+ = Е АС*3-

1=1

Решение на втором участке заднего фронта V = 5^3=1 С^4вЛ^ должно удовлетворять граничным условиям

V (-«>)= 0, V (£*-)= к, Пе_ = V |е+ .

Отсюда получаем

СивЛ^ + С24 вЛ2«* = К, С34 = 0,

3

С14А1вЛ1«* + С24А2вЛ2** = ^\С\з.

г=1

Постоянные С^4, г = 1, 2, 3, определяются однозначно. При постановке задачи и в процессе ее решения высказывались различные ограничения на выбор параметров. Отметим основные условия: V, < А, последнее соотношение формулы

д = ^ + \д2 К °’

где р = —\>с?у‘2Ь~‘2 — д = ^с3-у3_0~3 + ^к^суВ^2 + к2У~1Ё>~1 . Условие (28),

например, выполняется при достаточно малом к2. Если параметры выбраны так, что это неравенство заменяется на обратное, то два корня будут сопряженными комплексными. Задний фронт становится осциллирующим. Полученные в п. 2-5 явные решения кусочно-линейных уравнений могут служить начальными приближениями для более строгих постановок.

6. В качестве примера рассмотрим задачу определения характеристических чисел, отвечающих распространению потенциала действия в проводящей системе сердца в области предсердий и желудочков. Как выше указывалось, можно принять с = 10-6 фсм-2, = 10-7фс-1. Скорость V распространения возбуждения в указанной об-

ласти - 2-4 мс-1 [17]. Примем V = 300 смс-1. Подвижность к1 равна обратной величине удельного электрического сопротивления. В книге [16] приведено значение удельного сопротивления в покое для мембраны нейрона порядка 1000 Ом • см2. Для среднего значения неравновесного удельного сопротивления мембраны клетки проводящей системы сердца примем несколько большее значение р = 1111 Ом • см2. Тогда к1 = р-1 = 9 • 1(Г4 Ом-1 • см^2. Получим р = -3009 • 103 см~2, д = 2009 • 106 + ± • 105к2 см~3. Неравенство (28) равносильно соотношению ^д\р\^ > из которого следует к2 < 0.192 фсм-2 • с-2.

Аналитическое определение корней характеристического уравнения А3 + ^А2 + а2А + аз = 0, где для рассматриваемой задачи = 3 • 103, а2 = 9 • 103, аз = ^ • 105к2, выполним с помощью перехода к неполному кубическому уравнению у3 + ру + д = 0, используя замену А = у — ^а\. Для любого значения к2, удовлетворяющего полученному условию, можно найти тригонометрическое решение для корней у^ и А^, г = 1, 2, 3. В предельном случае при к2 = 0.192 получаем у\^2 = ~щ(—р)5 соэ = 1001.499,

уз = (-p)z = —2002.998. При этом = Уг — 1000. Непосредственное вычисле-

ние по тригонометрическому методу корней при 0 < к2 < 0.192 приводит к большим ошибкам, так как связано с определением угла а по cos а в области минимального значения последнего. Примем для эталонной задачи А1 = 2, А2 = 1, A3 = -3003. По теореме Виета получаем

А1 + А2 + A3 = —3000 = &1, А1А2А3 = —6006 = —аз,

А1А2 + А1А3 + А2Аз = -9007 ж а2.

Из второго соотношения находим неизвестный коэффициент поперечного дрейфа частиц к2 = 0.1818 фсм-2 • с-2.

Литература

1. Розенблатт Ф. Принципы нейродинамики. Перцептроны и теория механизмов мозга / пер. с англ. В. Я. Алтаева; под ред. С. М. Осовца. М.: Мир, 1965. 480 с.

2. Веденов А. А. Моделирование элементов мышления. М.: Наука, 1988. 160 с. (Современные проблемы физики).

3. Коренев Г. В. Введение в механику человека. М.: Наука, 1977. 264 с.

4. Филиппов А. Т. Многоликий солитон. М.: Наука, 1990. 298 с.

5. Новоселов В. С. Статистическая динамика: учеб. пособие. СПб.: Изд-во С.-Петерб. гос. ун-та, 2009. 393 с.

6. Новоселов В. С. Смешанное длинноволновое уравнение // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2010. Вып. 1. С. 65—72.

7. Hirota R. Exact N-soliton solutions of the wave equation of long waves in shallow-water and nonlinear lattices// J. Math. Phys. 1973. Vol. 14, N 7. P. 810—814.

8. Кунин И. Л. Теория упругих сред с микроструктурой. М.: Наука, 1975. 416 с.

9. Захаров В. Е. Метод обратной задачи рассеяния // Кунин И. Л. Теория упругих сред с микроструктурой. М.: Наука, 1975. С. 226-273.

10. Маханьков П. Г. Солитоны и численный эксперимент // Физика элементарных частиц атомного ядра. 1983. Т. 4, вып. 1. С. 123-180.

11. Фаддеев Л. Д., Тахтаджян Л. А. Гамильтонов подход в теории солитонов. М.: Наука, 1986.

528 с.

12. Додд Р., Эйлбек Д., Гиббон Д., Моррис Х. Солитоны и нелинейные волновые уравнения / пер. с англ. В. П. Гурария; под ред. А. Б. Шабата. М.: Мир, 1988. 694 с.

13. Takahashi R., Ohkawa T. Numerical experiment on interaction of solitons describing recurrence of initial state// Computational Mechanics (Springer-Verlag). 1989. P. 273-281.

14. Balashov V., Bladimirov V. Solitary waves in dispersive complex media. Berlin: Springer, 2005. 292 p.

15. Ризниченко Г. Ю. Математические модели в биофизике и экологии. М.; Ижевск: Ин-т ком-пьют. технологий, 2003. 184 с.

16. Жаботинский А. М. Концентрационные автоколебания. М.: Наука, 1974. 179 с.

17. Кринский В. И., Медвинский А. Б., Панфилов А. В. Эволюция автоволновых вихрей. М.: Знание, 1986. 48 с. (Сер. Математика, кибернетика. № 8).

18. Покровский А. Н. Процессы управления в нервных клетках. Л.: Изд-во Ленингр. ун-та, 1987.

85 с.

19. Keener J., Sneyd J. Mathematical physiology // Interdisciplinary Applied Mathematic. 1998. Vol. 8. P. 268-296.

20. Тамм И. Е. Основы теории электричества: учеб. пособие для вузов. М.: Наука, 1989. 504 с.

21. Новоселов В. С., Королев В. С. Аналитическая механика управляемой системы: учеб. пособие. СПб.: Изд-во С.-Петерб. гос. ун-та, 2005. 298 с.

Статья рекомендована к печати проф. Ю. М. Далем.

Статья принята к печати 19 мая 2011 г.