CC BY
52УДК 531/534:57
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ КОРРЕКЦИИ ВЫХОДНОГО СИГНАЛА С ГРАВИТОИНЕРЦИАЛЬНОГО МЕХАНОРЕЦЕПТОРА
ВЕСТИБУЛЯРНОГО АППАРАТА
В. А. Садовничий1, В. В. Александров2, Т. Б. Александрова3, A.A. Коник4, В. Б. Пахомов5, Г. Ю. Сидоренко6, Э. Сото7, К. В. Тихонова8, Н.Э. Шуленина9
Представлена модифицированная упрощенная математическая модель Ходжкина-Хаксли активности первичного афферентного нейрона вестибулярной системы, и проведен анализ этой модели.
Ключевые слова: биомехатронная система, вестибулярный механорецептор, вестибу-лосенсорный конфликт.
A modified simplified Hodgkin-Huxley mathematical model of activity of the primary-afferent neuron of the vestibular system and an analysis of this model are presented.
Key words: biomechatronic system, vestibular mechanoreceptor, vestibular sensory conflict.
Введение. Развитие систем полуавтоматического управления движением за последние тридцать лет привело к формированию нового направления в механике управляемых систем — биомехатроники [1]. Под биомехатронной системой (ВМС) будем понимать систему состоящую из трех частей — подвижного объекта, электронной части и человека или части его организма — и предназначенную для оказания помощи человеку с целью измерения характеристик объекта и (или) управления движением (реальным или виртуальным).
Ярким примером прогресса в разработке биомехатронных систем является создание и эксплуатация авиационных и космических тренажеров. При этом математическое обеспечение такого рода ВМС занимает одну из главных позиций [2].
В данной статье рассматривается возможность создания конкретной ВМС — корректора выходной информации гравитоинерциального механорецептора [3] вестибулярного аппарата.
В 1992-2001 гг. И.Б. Козловская и ее коллеги в экспериментах на орбитальной станции "Мир" и МКС обнаружили, что развитие вестибулосенсорного конфликта в условиях микрогравитации приводит к увеличению времени установки взора в 3 и более раз по сравнению с Землей [4]. К.В. Тихоновой было показано [5], что использование микроакселерометра, установленного на шлеме космонавта, может уменьшить эту задержку при моделировании гальванической стимуляции. В 2007 г. японские специалисты (см. [6]) представили результаты экспериментов по возможности изменения персональной ориентации в пространстве с помощью гальванической стимуляции вестибулярного аппарата. И.В. Орлов с сотрудниками в 2003 г. опубликовали статью [7] о возможности коррекции вертикальной позы человека с помощью биологической обратной связи.
1 Садовничий Виктор Антонович — доктор физ.-мат. наук, проф., академик РАН, ректор МГУ, e-mail: vladimiralexandrov366Qhotmail.com.
2 Александров Владимир Васильевич — доктор физ.-мат. наук, проф., зав. каф. прикладной механики и управления мех.-мат. ф-та МГУ, e-mail: vladimiralexandrov366Qhotmail.com.
3 Александрова Тамара Борисовна — канд. биол. наук, ст. науч. сотр. Центра новых информационных технологий МГУ, e-mail: vladimiralexandrov366Qhotmail.com.
4 Коник Анастасия Алексеевна — студ. каф. прикладной механики и управления мех.-мат. ф-та МГУ, e-mail: vladimiralexandrov366Qhotmail.com.
5 Пахомов Виктор Борисович — науч. сотр. НИИ механики МГУ, e-mail: vladimiralexandrov366Qhotmail.com.
6 Сидоренко Галина Юрьевна — канд. физ.-мат. наук, ассист. каф. прикладной механики и управления мех.-мат. ф-та МГУ, e-mail: gingrulQyandex.ru.
7 Сото Энрике — PhD, зав. лаб. Автономного университета г. Пуэбла, Мексика, e-mail: esotoQsiu.buap.mx.
8 Тихонова Катерина Владимировна — науч. сотр. Института математических исследований сложных систем МГУ, e-mail: vladimiralexandrov366Qhotmail.com.
9 Шуленина Нейля Энверовна — канд. физ.-мат. наук, ст. науч. сотр. лаб. математического обеспечения имитационных динамических систем мех.-мат. ф-та МГУ, e-mail: neshulQyandex.ru.
1. Постановка задачи. Рассмотрим функциональную схему гравитоинерцнального механорецеп-тора [3] с горизонтальной осью чувствительности, являющегося составной частью отолитова органа вестибулярного аппарата (рис. 1). Выходной сигнал этого инерцнального биосснсора нрсдставляст собой чередующиеся пачки импульсов (снайков) с неременной частотой следования.
а Динамика отолитовой массы X Механо-элсктричсская транедукция 'тг Динамика волосковой клетки V\ Синаптическая /Syn Активация первичного нейрона
трансмиссия
Рис. 1. Функциональная схема гравитоинерциального механорецептора: а — линейное ускорение, х — отклонение отолитовой массы, /тг — ток транедукции, V — мембранный потенциал волосков ой клетки, У2 — мембранный
потенциал первичного нейрона
Таким образом, наличие снайков и изменение временных интервалов между ними является описанием выходной информации о механическом стимуле.
Самая простая схема ВМС коррекции выходного сигнала состоит из трех
Горизонтальный
микро-аксслсромстр
ка
Аппарат гальванической стимуляции
'galv
Афферентный первичный нейрон
Г,( О
Рис. 2. Функциональная схема биомехатроииой системы коррекции
выходного сигнала
блоков (рис. 2). В электронной части этой биомехатронной системы (блок 2) должен быть реализован алгоритм преобразования выходного идеального сигнала микроакселерометра ускорения а в гальванический стимул /gaiv, корректирующий выходной сигнал V2(t) первичного нейрона. В настоящее время есть две гипотезы о воздействии гальванического тока на вестибулярный аппарат [8|: непосредственное влияние на первичные нейроны, когда имеет место совместное воздействие еинаптичеекого и гальванического токов; опосредованное влияние, оказываемое вначале на волосковые клетки, а затем на первичные нейроны. Схема ВМС, представленная на рис. 2, соответствует первой гипотезе. В рамках этой схемы необходимо ответить на вопрос: что делать с инструментальными погрешностям«, реально присутствующими в выходном сигнале микроакселерометра?
Рассмотрим эту задачу в следующей постановке. В результате гальванической стимуляции на входе первичного нейрона присутствует аналог еинаптичеекого тока Iga\v, состоящий из суммы двух компонент, а именно детерминированной составляющей (кусочно-постоянного сигнала с акселерометра, соответствующего механическому стимулу) и белого шума постоянной интенсивности, являющегося традиционной моделью аддитивных инструментальных погрешностей. Требуется составить математическую модель активности первичного нейрона и получить статистическую оценку влияния /gaiv на выходной сигнал с этой модели.
2. Модифицированная модель Ходжкина—Хаксли. В качестве математической модели активности первичного нейрона рассмотрим детерминированное описание в среднем стохастического процесса образования импульсов (епайков), представляющих собой релаксационные автоколебания постоянной амплитуды и переменной во времени частоты. Распространение этих автоколебаний в виде автоволн не рассматривается. За основу взята классическая модель Ходжкина Хаксли [9], которая была модифицирована и упрощена. Модификации были сделаны по результатам проведенных экспериментов в лаборатории нейрофизиологии Автономного университета штата Пуэбла (Мексика). Перечислим эти модификации и упрощения.
A. Данные, используемые в модели, были получены в опытах на млекопитающих.
B. При описании калиевого тока /к возникла необходимость во введении параметр а инактивации Нк и в связи с этим использования еще одного уравнения Колмогорова для описания динамики этого вероятностного параметра. В данной статье рассматривается зависимость только от потенциала действия Нк = Нк (V2).
C. Математически было обосновано уменьшение порядка математической модели в случае наличия малого параметра тт (постоянной времени натриевого тока /Na) при производной подсистемы, описывающей динамику параметра активации m этого тока.
D. Интеграл n + hNa = C = const, используемый физиологами при упрощении модели Ходжкина-Хаксли, модифицирован, согласно данным экспериментов, в виде n + hNa = C(V2).
Таким образом, модифицированная модель Ходжкина Хаксли представляется системой в форме Ко-
ши третьего порядка с малым параметром в правой части уравнения, описывающего динамику вероятностного параметра Лк [10]. Явная зависимость этого параметра от времени в данной статье не рассматривается, и поэтому упрощенная и модифицированная модель Ходжкина-Хаксли является системой второго порядка (малый параметр равен нулю):
С2 • = Ьуп — — 1к — 1ь2,
Тп(У2) ■ = {Поо(У2) - п)(5ю,
где
/ка = С(У2) - и)(У2 - Ека), /к = дки4Лк<х(У2 - Ек), 1ь2 = дь2(У2 - ЕЬ2),
С (У2) = и^(У2) + Лк
а те (У2).
С2
/Ма — натриевый ток, дка — максимальное значение проводимости и Ека — потенциал реверса натриевого тока;
/к — калиевый ток, дк — максимальная проводимость и Ек — потенциал реверса калиевого тока; /¿2 — ток утечки в месте формирования афферентного импульса, дь2 — его максимальная проводимость, а ЕЬ2 — потенциал реверса тока утечки; и
ствия частиц активации калиевого тока в калиевых каналах;
Лк — параметр, описывающий процесс инактивации калиевого тока, являющийся вероятностью отсутствия частиц инактивации калиевого тока;
Лка — параметр, описывающий процесс инактивации натриевого тока, являющийся вероятностью отсутствия частиц инактивации натриевого тока;
тп — постоянная времени процесса активации калиевого тока;
и<х, — стационарные значения процессов активации калиевого и натриевого тока соответственно; го) Лк— стационарные значения процессов инактивации натриевого и калиевого тока соответственно;
Ql0 — коэффициент температурной зависимости, необходимость введения которого связана с разностью физиологической температуры и комнатной (см. таблицу), при которой проводились эксперименты по определению параметров модели.
Функциональные параметры модифицированной и упрощенной модели активности первичного нейрона имеют вид
т°°(1У = ТТ-/-№+33.8) А' ^.оо(^) = —--,
1 + ехр ) 1 + ехр ^)
ПхШ =-1--= о.«мсв-о.таа9
1+е*р(^Ш2)' 1+ехр(й±а5|Е)
Тп(У2) =-т—-т-——С(У2)=п00(У2) + }1Шоо(У2), <210(а,Т)=аУ 10 Л
еХр(^)+ехр(М)
Параметр Численное значение Доверительный интервал Размерность
с2 1 _ мкФ см
■ЕфГа 52 _ мВ
Ек -84 _ мВ
ЕЬ2 -63 _ мВ
да а 2,3 2-8 мСм/см2
ок 2,4 1-2,6 мСм/см2
дь2 0,03 0,02-0,16 мСм/см2
а 3 3-4 _
Т 37 _ град
То 20 _ град
Численные параметры представлены в таблице.
3. Анализ необходимости присутствия белого шума в ВМС коррекции выходного сигнала с гравитоинерциального механо-рецептора. Будем считать, что измерения микроакселерометра идеальны и представлены кусочно-постоянной функцией ка(^), состоящей из двух ступенек. Значение первой ступеньки берется из левой окрестности первой точки бифур-
2
-40 -20 0
Рис. 3. Особью точки н продольный цикл математической модели: кружочек числеиио найденная неустойчивая особая точка при /а^т = 1,8 мкА/см2 (неустойчивый фокус); линия — предельный цикл, соответствующий этому же значению входного сигнала: звездочка единственный аттрактор при /Бцт = 0,5 мкА/см2 (устойчивый узел)
ступеньки из правой окрестности. Таким образом, точка бифуркации рассматривается как порог чувствительности. Выходной сигнал второго блока в схеме ВМС имеет вид /stim = ka(t) + rG(t), где второе слагаемое есть белый шум (или его аппроксимация) с постоянной спектральной плотностью, равной r2. В соответствии с постановкой задачи необходимо ответить на вопрос: может ли белый шум оказывать положительное влияние на реакцию первичного афферентного нейрона в ответ на механический стимул?
1. Положительное влияние в большинстве случаев связано с возможностью существования стохастического резонанса [11]. Например, в нелинейных бнетабнльных еи-
r
приводящую к увеличению отношения сигнала на выходе системы к шуму или к необходимости отыскания максимума соответствующего функционала (статистической оценки спектральной плотности выходного сигнала и др.). Рассматриваемая модель при /stim = ka(t) является нелинейной и бнетабильной в обобщенном смысле, так как на первой ступеньке существует одна асимптотически устойчивая точка покоя, а на второй ступеньке существует более сложный аттрактор, являющийся предельным циклом (рис. 3).
Следует отметить очень удобное расположение этих аттракторов на фазовой плоскости (рис. 3), не препятствующее, а, наоборот, способствующее переходу из окрестности точечного аттрактора в окрестность предельного цикла.
Таким образом, на первый взгляд нет острой необходимости привлекать стохастический вспомогательный процесс для преодоления препятствий при переходе из окрестности одного аттрактора в окрестность другого, как это бывает при наличии сепаратрис [11].
2. Для моделирования стохастической составляющей rG(t), опираясь на [11], используем процесс
Орнштейна Уленбека U(t) с корреляционной функцией Ru(s) = ¿exp^—дисперсией ^ и односторонней спектральной плотностью ФitS(s) = i+^V (которая связана со спектральной плотностью следующим образом: Ф®5'^) = Фu(s) + Ф[г(—s)). Тогда процесс G(t) = (2c)^U(t) будет характеризоваться односторонней спектральной плотностью Ф^5^) = .
Результаты, представленные ниже, соответствуют значению параметра c = 0,05, при котором можем рассматривать rG(t) как процесс со спектральной плотностью, близкой к постоянному значению y2 = 2cr2. Заметим, что размерность коэффициента r равна мкА/см2, а процесс G(t) безразмерен.
Для G(t) применима аппроксимация рядом Каца Шинозукн Gn(í) = E¿Iicos(w¿W + ¥>o¿)>
где фазы равномерно распределены в интервале [0; 2п], распределение частот (w) характеризуется плотностью вероятности /w(s) = ^ Ф^5^) = ^ 1 _^¿2S2 ■ Расчеты в среде Matlab были выполнены при значении N = 100.
3. На рис. 4 показаны два графика в отсутствие и при наличии белого шума. Появление одиночных спайков на первой ступеньке, соответствующей отсутствию необходимого стимула, свидетельствует о хорошей чувствительности рассматриваемой модели первичного нейрона. На второй ступеньке можно заметить увеличение числа спайков на рассматриваемом интервале.
В связи с этим возьмем в качестве функционала статистическую оценку математического ожидания числа спайков на фиксированном временном интервале второй ступеньки при заданной интенсивности.
Проводилось интегрирование системы уравнений второго порядка с использованием функции Matlab ode45() при постоянном значении ka(t) = 1,8 мкА/см2, различных реализациях rG(t), соответствующих различным значениям r, и при одном и том же начальном условии (V2,n) = (-45,66; 0,11). В каждом
r
r
(the number of spikes)). Наличие максимума позволяет утверждать о возможности подбора оптимальной интенсивности белого шума во втором блоке корректора. На рис. 5, а представлена зависимость еред-
r
тервалов, соответствующая значению r = 636,4 мкА/см2. Гистограммы строились путем объединения подсчетов по 10 решениям, аналогичным рассмотренным выше, но найденным на интервале 150 мс при ka(t) = 1,5 мкА/см2 и c = 0,0001. Таким образом, с учетом шума на входе на выходе получили решения, характеризующиеся наличием спайков, но с иным распределением межспайковых интервалов, нежели при rG(t) = 0. Случаю отсутствия шума соответствует решение с постоянными межспайковыми интервалами в 22,3 мс и гистограмма, состоящая из одного столбца.
Рис. 4. Реакция первичного нейрона на входной сигнал при отсутствии и при наличии белого шума: а решение системы без шума, б решение системы с шумом интенсивности г = 22,14 мкА/см2 для входного сигнала — ступеньки со значениями, описанными
для рис. 3
Рис. 5. Результаты статистического анализа: а зависимость статистической оценки математического ожидания числа спайков на временном интервале второй ступеньки от инг
г = 636,4 мкА/см2
Наличие нескольких максимумов в гистограмме на рис. 5, б говорит о том, что полученная модель активности афферентного первичного нейрона лучше модели Фитцхью Нагумо, так как качественно совпадает с экспериментами на первичном нейроне (гистограммы мультимодальны), в отличие от гистограммы, построенной при воздействии белого шума на модель Фитцхью Нагумо, которая является унимодальной [11].
4. Заключение. Таким образом, при идеальных показаниях акселерометра алгоритм во втором блоке корректора формируется по формуле = коа(Ь) + ГоС(£), где го есть оптимальная интенсивность белого шума, соответствующая постоянному ускорению на временном интервале \Ък, ^+1]- Полученный результат является первым шагом в проектировании различного рода биомехатронных систем, позволяющих корректировать персональную ориентацию в пространстве: стабилизацию взора в экстремальных условиях (микрогравитация и др.); вертикальную позу для лиц с вестибулярной дисфункцией и для лиц пожилого возраста. Этот же результат можно использовать для разработки алгоритмов управления динамическими имитаторами с целью тестирования прототипов вестибулярных протезов. Кроме того, представленная модель активности первичного афферентного нейрона вестибулярного аппарата, полученная в результате клеточных экспериментов на млекопитающих, дает возможность сравнить данный результат ее анализа с другими результатами по численному анализу классической модели Ходжкина Хаксли (см. библиографию в [12]).
Работа выполнена при поддержке государственного контракта № 11.519.11.2045 и гранта РФФИ № 13 01 00515.
СПИСОК ЛИТЕРАТУРЫ
1. Садовничий В.А., Александров В.В., Лемак С.С., Шкель A.M. О биомехатронике: Мат-лы научной школы-конференции "Мобильные роботы и мехатронные системы". М.: Изд-во МГУ, 2005. Ч. 1. 53-60.
2. Александров В.В., Воронин Л.И., Глазков Ю.Н., Ишлинский А.Ю., Садовничий В.А. Математические задачи динамической имитации аэрокосмических полетов. М.: Изд-во МГУ, 1995.
3. Александров В.В., Александрова Т.В., Мигунов С. С. Математическая модель гравитоинерциального механоре-цептора // Вестн. Моск. ун-та. Матем. Механ. 2006. № 2. 59-64.
4. Tomilovskaya E.S., Berger М., Gerstenbrand F., Kozlovskaya I. В. Effects of long-duration space flight on target acquisition // Acta Astronautica. 2011. 68. 1454-1464.
5. Тихонова К.В. Коррекция стабилизации взора в экстремальных условиях визуального управления движением: Магистер. дис. М., 2011.
6. Садовничий В.А., Александров В.В., Александрова Т.Е., Сото Э., Сидоренко Г.Ю., Тихонова К.В. Об автоматической коррекции вестибулосепсорпого конфликта в условиях невесомости, основанной на принципе гальванической стимуляции и на компьютерном моделировании // Интеграл. 2012. № 2. 70-74.
7. Орлов И.В., Гусев В.М., Долгобродов С.Г., ГЦупляков B.C. О возможности коррекции вертикальной позиции человека с помощью биологической обратной связи // Сенсорные системы. 2003. 17, № 1. 58-67.
8. Kim K.S., Minor L.B., Delia Santina С.С., Lasker D.M. Variation in response dynamics of regular and irregular vestibular-nerve afferents during sinusoidal head rotations and currents in the chinchilla // Exp. Brain Res. 2011. 210, N 3-4. 643-649.
9. Hodgkin A.L., Huxley A.F. A quantitative description of membrane current and excitation in nerve //J. Physiol. 1952. 117. 500-544.
10. Александров В.В., Михалева Е.Ю., Сото Э., Гарсиа-Тамайо F. Модификация математической модели Ходжки-на-Хаксли первичного нейрона вестибулярного аппарата // Вестн. Моск. ун-та. Матем. Механ. 2006. № 5. 65-68.
11. Simiu Е. Chaotic Transitions in Deterministic and Stochastic Dynamical Systems. Princeton: Princeton Univ. Press, 2002.
12. Tuckwell H.C., Jost J. Analysis of Inverse Stochastic Resonance and the Long-Term Firing of Hodgkin-Huxley Neurons with Gaussian White Noise. Leipzig: Max Planck Inst., 2012.
Поступила в редакцию 22.03.2013
УДК 531.8
О МЕТОДАХ ИНЕРЦИАЛЬНОЙ ГРАВИМЕТРИИ Ю.В. Болотин1, А. А. Голован2
В обзорной статье приводится краткое описание предыстории, современного состояния и перспектив развития инерциальной гравиметрии в России, как их видят авторы статьи — непосредственные свидетели становления этого нового направления геофизики. Основное внимание уделяется информационной, алгоритмической, составляющей авиационной гравиметрии.
Ключевые слова: инерциальная гравиметрия, аэрогравиметрия, алгоритм оценивания.
The paper discusses the history and modern state of inertial gravimetry in Russia, as seen by the authors who have witnessed developments in this field of geophysics. More attention is given to information and algorithmic aspects of airborne gravimetry.
Key words: inertial gravimetry, airborne gravimetry, estimation algorithm.
Введение. Инерциальная гравиметрия — это прикладная наука об определении силы тяжести Земли по движению чувствительной массы. Основоположником инерциальной гравиметрии можно считать
1 Болотин Юрий Владимирович — доктор физ.-мат. наук, проф. каф. прикладной механики и управления мех.-мат. ф-та МГУ, e-mail: ybolotinQyandex.ru.
2 Голован Андрей Андреевич — доктор физ.-мат. наук, зав. лаб. управления и навигации МГУ, e-mail: aagolovanQyandex.ru.
