Моделирование мембранного потенциала клетки сино-атриального узла Текст научной статьи по специальности «Математика»

Вестник Сыктывкарского университета. Сер Л. Вып. 9. 2009

УДК 539.3

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

В.Н. Тарасов

Работа сердца связана с регулярным сокращенны мышц образующих сердечные стенки. Ритмичное сокращение этих мышц вызвано распространением электрического потенциала действия (ПД) по сердечным волокнам. Большинство этих волокон способны самовозбуждаться. Этим он отличается от ПД в нервных клетках: там для возбужения ПД необходим некоторый сигнал (ток инициализации). Участок сердца, дающий возбуждение максимальной величины- это сино-атриальный узел (СА). В настоящее работе предлагается математическая модель электрической активности клеток С А.

Построение математической модели

Известно, что клетки СА генерируют ПД, который определяется различными концентрациями ионов на внутренней и внешней сторонах мембраны. Таким образом мембрану можно рассматривать как конденсатор, имежду зарядом и напряжением выполнено равенство

Я = си, (1)

где заряд, II— напряжение, с-электроемкость мембраны. Дифференцируя равенство (1) получим

¿и

1 = (2)

где /- ток, проходящий через мембрану. Измеряемой величиной обычно является разность потенциалов на сторонах мембраны (мембранный потенциал ), график типичноного потенциала в зависимости от времени приведен на рис 1. Функцию, изображенную на рис 1 обозначим через Ф(£) (экспериментальная кривая). Данный график заимствован из [2].

© Тарасов В.Н., 2009.

рис. 1

рис. 2

Из результатов экспериментов следует, что формирование мембранного потенциала представляет собой некоторый автоколебательный пр-цесс. Общий ток через мембрану складывается из ионных токов (калиевый, натриевый, кальциевый и другие) ионы проходят через каналы в мемране (для каждого типа ионов свои каналы). Проводимость этих каналов зависит от мембранного потенциала, кроме того величина ионных токов зависит, разумеется, от концентрации ионов на внутренней и внешней стороне мембраны. Также имеются ионные насосы, которые переносят ионы против градиента электрохимического потенциала, при этом затрпачивается энергия.

В клетках С А происходит генерация незатухающих колебаний А. А. Андронов назвал такие системы автоколебательными [1], и также показал что автоколебания невозможны в системах, описываемых диф-ферециальными уравнениями первого порядка и в линейных системах (т.е. в системах, описываемых линейными дифференциальными уравнениями).

Из сказанного выше можно предположить, что клетки С А представляют собой некоторую нелинейную электрическую схпему (с точки зрения генерации мембранного потенциала), описываемой системой дифференциальных уравнений:

где Д(?7, /), /2 (С/, I) достатачное число раз дифференцируемые функции. Будем предполагать, что они предсталяют собой многочлены двух переменных некоторой степени. Предполагается также, что система является автономной, т.е. правые части системы не содержат явно время.

Теперь задача ставится следующим образом: подобрать коэффициенты полиномов в правых частях системы (3) таким образом, чтобы

(3)

функция U(t) аироксимировала наблюдаемую кривую, изображенную на рис.1 с достаточной степенью точности. При этом желательно, чтобы степень полиномов была как можно меньше.

На рис.1 время измеряется в милисекундах t G [0,1000], потенциал в миливольтах (mv). Удобнее работать с безразмерными величинами, поэтому сделаем замену переменных

1000т Ггло п , ч Ф(т)

* = —, г € [0,2т], =

В данном случае Ф(0) = —66.8mv и является минимальным значением экспериментального потенциала Ф(£).

Одним из простейших и самых изветных уравнений, описывающих нелинейные автоколебания является уравнение Ван Дер Поля

^-„(1-„»)£ + „„ = „. (4)

Это уравнение используется, например, в [6], однако применение его там ничем не обосновано. Решение уравнеия (4) симметрично относительно оси т и форма автоколебаний не совпадает с наблюдаемой.

Для моделирования мембранного потенциала будем использовать несколько более сложное уравнение

П 1 / (17/

_ = &(„_ /32)2 + (03 + /34(и - /32) +fo(u- /32)2) — - Аз (« - fo), (5) или эквивалентную ему систему

% = рг(и- /32)2 + (/З3 + - /32) + /55(и - /32)2)у - /36(и - /32)

Тг=У-

(6)

Уравнение (5) имеет единственное неустойчивое положение равновесия и — ¡52« Решение системы (6) зависит от коэффициетнтов /3 = (/?!, /5в), т.е. и(т) — и(т;/3). Для того, чтобы определить коэффи-

циенты /3 возьмем точки Т{ — 27гг/т, г £ [0 : т], где т натуральное число, и введем в рассмотрение функцию

т -

/(/3) = тш(у/(и(ъ /3) - Фз))2 + (п - Т,)2+

г=0 3

т -

+ гшп(у/(и(тг, Р) - + (П - г,)2. (7)

э=о

Вектор коэффициентов находим, решая задачу оптимизации:

/(/3) = ппп/(/3).

(8)

Р

Функция /(/3), вобще говоря не является непрерывно дииферецируе-мой, поэтому для решения задачи (8) применялся метод Монте-Карло. Для вычисления значений функции /3) в точках сетки использовался метод Рунге-Кутта четвертого порядка точности, причем на интервале интегрирования [0 : 2тг] число шагов ТУ = 4000, в формуле (7) т = 400. Программа была составлена на языке проограммирования ФОРТРАН-ЭО.

В результате решения задачи оптимизации было получено следующее решение.

/3 = (-25.98179, -0.284866, 2.4379215, -5.9363316,35.8858104,19.5746535)

На рис.2 а- экспериментальная кривая (в безразмерном виде), Ь-функция п(т), являющеяся решением системы уравнений (6). На рис.3 а производная экспериментальной кривой, Ь -у(г), производная п(т), полученная в результате решения системы дифференциальных уравнений (6). На рис.4 Ь- фазовая кривая (и(т), у(т) т Е [0, 2тг]) -предельный цикл системы (6), а- экспериментальная кривая (^(г), ^Т(т) т Е [0, 2тг]) . Используя методы [1] можно показать, что предельный цикл системы (6) является ассимтотически устойчивым.

Для дифференцирования функций, определяемых по экспериментальным данным использовался следующий метод. Пусть

Ь

рис. 3

рис. 4

Ясно, что ф{х) (р(х) при с —У оо. Тогда

Замечание 1. Коэффициенты системы уравнений (6) (или, что тоже самое (7)) получаются большими по модулю. Не смотря на величину коэффициентов при интегрировании системы (6) неустойчивости счета не наблюдается. Результаты работы Фортран-программы проверены в системе МАРЬЕ 8, при помощи которой нарисованы графики.

Замечание 2. Функция /(/3) отличается от обычной функции метода наименьших квадратов, потому что решение и(т; /3) на некоторых участках имеет большие производные, и использование более простой целевой функции вида /(/3) = (Р(тз))2 не пРивел0 ни како-

му разумному результату. К тому же расстояние между двумя кривыми разумно определить следующим образом. Пусть

71 = {^1 = = Ыт)}, 72 = {Ф\ = ФЛт),Ф2 = ^2(т)}, т е [а, Ь]

-две кривые, заданные параметрическим образом. Тогда

Иными словами для любой точки, лежащей на кривой 71 отыскивается ближайшая ей точка кривой 72 и наоборот и все эти расстояния суммируются вдоль всей кривой. Такое определение расстояния не зависит от параметризации кривой.

Замечание 3. Математических моделей электрической активности мембран существует достатачно много. Чаще всего для их построения модифицируется (усложняется с целью более детального описания ) модель Ходжкина-Хаксли. Такие модели предлагаются например в [2,3].Однако , как и в модели Ходжкина-Хаксли, вероятности открывания и закрывания ионных каналов (в модели Ходжкина-Хаксли переменные т, 77, К) не не удается отождествить с каими-либо реальными компонентами мебраны ( [4]). Таким образом все эти модели также являются феноменологическими. Также автору данной работы не уда-лоссь получить авоколебания используя уравнения модели как [2], так и [3]. Повидимому способ описания мембранных потенциалов в клетках , способных самовожбуждаться, и в нейронах должен быть различным.

пни

ге[а,ь]

тт - ФЛт))2 + (М*) - Мт))2)Лт.

Литература

1. Андронов А.А., Витт А.А., Хайкин С.Э. Теория колебаний. М.: Физматгиз. 1959. 999 с.

2. Zhang Н, Holden A.V., Kodama I, Honjo H, Lei M, Varghese T, Boyett M.R. (2000) Mathematical models of action potentials in the periphery and center of the rabbit sinoatrial node. Am J Physiol Heart Circ Physiol 279:H397-H421

3. R. Wilder, H.J. Jonqsma and A.C.G. van Ginneken. Pacemaker activity of the rabbit sinoatrial node// Biophys. J. Biophysical Society Volume 60, 1991. P. 1202-1216.

4. Буравцев B.H. Модель авторегуляции в системе аксона кальмара// Биофизика, 1999. Т 44, выи 5. С. 892-897.

5. Карнаухов А.В., Карнаухова Е.В. Новый метод идентификаци-ии систем на основе критерия минимальной квадратичной невязки для задач биофизики// Биофизика, 2004. Т 49, выи 1. С. 88-97.

6. Богачева И.Н., Кучер В.И., Щербакова Н.А., Мусиенко П.Е., Герасименко Ю.П. Математическое моделирование процессов формирования локомоторных патертонов при эпидуральной стимуляции спинного мозга с учетом периферической обратной связи// Биофизика, 2005. Т 50, выи 6. С. 1125-1130.

Summary

Tarasov V.N. Modeling the cell membrane potential of the sinoatrial node

The work of the heart associated with a regular contractions of muscle forming the heart wall. The heart rate caused by the spread of the electrical action potential (AP) on heart fibers. Most of these fibers can self-excite. This is in contrast to AP in the nerve cells. For self-excitation AP in neuronal cells requires a signal (current initialization). Sector of the heart, giving the excitation maximum value - it is sinoatrial node. In this paper a mathematical model of electrical activity of cells of the sinotrial node is proposed.

Отдел математики КНЦ УрО РАН

Поступила 28.03.2009