CC BY
70
УДК 531:518:577
Вестник СПбГУ. Сер. 10, 2013, вып. 4
В. С. Новоселов
О МАТЕМАТИЧЕСКОЙ МОДЕЛИ ВОЗБУЖДЕНИЯ КЛЕТОК СЕРДЦА
1. Мышечные клетки сердца выполняют три функции. Одни клетки находятся в колебательном режиме, другие возбуждаются и распространяют электрический потенциал. Наконец, остальные возбуждаются и под действием указанного потенциала сокращаются, способствуя этим выталкиванию крови. Каждая клетка специализирована на контакт с близлежащими клетками.
Активный потенциал генерируется в клетках первого типа, называемых водителями ритма, пейсмекерами или Р-клетками. Они входят в состав всех отделов проводящей системы сердца [1—6]. Передача потенциала от Р-клеток к передающим клеткам, а также на заключительном этапе к клеткам миокарда осуществляется, как правило, с помощью электрического синапса, при котором пре- и послесинаптическая мембраны фактически соприкасаются. Но даже в этой узкой синаптической щели имеются белковые мостики, т. е. каналы с беспрепятственным перемещением ионов и небольших молекул. В статье [7] дано физическое обоснование применения длинноволнового уравнения Кортевега-де Фриса для математического описания распространения бегущего импульса в нервных волокнах и рассмотрены кусочно-линейные имитирующие модели. В работе [8] построена автоколебательная математическая модель пейсмекера. В настоящей статье, являющейся продолжением отмеченных двух исследований, разработана математическая модель возбуждения клеток миокарда с учетом структуры бегущего импульса и кинетических уравнений мышечного сокращения.
2. Как следует из графиков, приводимых в литературных источниках, например [9], профиль бегущей волны передающих клеток проводящей системы сердца (второй тип клеток) примерно симметричный, имеет амплитуду около 80 мВ, длительность приблизительно 0.2 с. Клетки второго типа выполняют роль аксона нейрона [5, 6]. Однако у аксона характер сигнала определяется видом серии одинаковых узких импульсов. В клетках второго типа проводящей системы сердца каждую секунду возбуждается существенно менее крутая одиночная волна. За исключением клеток, входящих в волокна Пуркинье и непосредственно примыкающих к сократительным мышечным клеткам желудочков, для математического ожидания можно использовать уравнения Ходжкина-Хаксли, Кортевега-де Фриса и ФитцХью-Нагумо [7, 9]. Особое положение занимают волокна Пуркинье. Они содержат ансамбли как пейсмекерных клеток, так и передающих клеток. Причем последние через электрический синапс должны передавать на контактную сократительную клетку достаточное количество положительных ионов, чтобы привести в возбуждение все связанные с ней сократительные клетки. Как следует из [8], в этом большую помощь могут оказать пейсмекерные клетки.
3. Клетки третьего типа, составляющие миокард предсердий и желудочков, называемые кардиомиоцитами, имеют сократительный механизм, подобный поперечнополосатым (скелетным) мышцам. Механизм сокращения этих мышц сравнительно хорошо изучен [5, 6, 10]. Сокращение единицы мышечной структуры (саркомера) происходит за счет скольжения тонких протофибрилл (актиновых) относительно толстых
Новоселов Виктор Сергеевич — доктор физико-математических наук, профессор, 199034, Санкт-Петербургский государственный университет; e-mail: novoselov@apmath.spbu.ru.
© В. С. Новоселов, 2013
(миозиновых), которое осуществляется в результате соединения головок (мостиков) миозина, предварительно активизированных молекулами АТФ (аденозинтрифосфат), с активным центром тонкой протофибриллы. Чтобы указанный центр активизировать, требуется некоторое количество ионов Са2+, которые, присоединяясь к закрывающему центр тропонину, вызывают смещение последнего. Отсоединение миозиновой головки от актина требует присоединения еще одной молекулы АТФ.
В книгах [5, 6] на основе рассмотрения соответствующей ферментативной реакции построена статистическая математическая модель сокращения поперечно-полосатых мышц. Пусть в перекрывающей области саркомера между актином и миозином число миозиновых головок совпадает с числом способных к активизации центров актина и равно N(£). Среди них х\ - число миозиновых головок, связанных с актином, Х2 - число миозиновых головок, не связанных с актином, но присоединивших АТФ, хз - число свободных головок, Х4 - число активизированных центров, не связанных с миозино-выми головками, Х5 - число свободных центров актина, Хб - число активизирующих групп ионов Са2+, Х7 - число свободных молекул АТФ, Х8 - число молекул АДФ (аде-нозиндифосфат). Приведем дифференциальные уравнения модели
^ = ах2хА - /Зх1хг + ^жхЛ^1,
_ I ЛГ-1
dt = 1Х3Х7 - a.X2Xi + —X2N , ^=I3x1X7-1X3X7 + ^X3N-\ (1)
^ = ях5х6 - ах2хА + yixiN-1, _ „,__„„._„.„ лг-1
/Зх\х7 — жх 5Ж6 + N~
Хх8 ~ (Зх-ixr - 7*3*7 — х2)N~
^ = 2/3xixr -ХХ8-
Здесь а, в, Y, к - неотрицательные кусочно-постоянные скорости бимолекулярных реакций, х - вероятность ресинтеза АДФ в АТФ в митохондриях клетки в единицу времени, приходящаяся на одну молекулу АДФ. В рассматриваемой задаче в качестве саркомера принимается кардиомиоцит. Уравнения (1) имеют интегралы
xi + Х2 + Х3 = N, xi + X4 + X5 = N, xi + X2 + X7 + X8 = s = const > 0, (2)
где s - общее число молекул аденозина, образующих вместе с остатками фосфорной кислоты АТФ и АДФ клетки.
При N = const уравнения (1) имеют [5, 6] равновесное решение
x2 = W-X1(1 + /37-1), *3 = /37-1*!, *4 = (3)
ax2(N—xi) o—l —1 о —1
ж5 = aX2+XX(i > xr = «Р xi x8 = 2ax X2X4.
Величина Xi определяет давление сжатия соответственно предсердий и желудочков. Действительно, пусть f - продольная составляющая силы, создаваемая присоединенной к актину одной головкой миозина. Тогда вклад одной клетки миокарда в общую силу сжатия будет пропорционален величине Xif.
Уравнения (1) должны быть дополнены уравнениями, выражающими баланс Ca2+, т. е. фазовой переменной xq. В межклеточной среде концентрация Ca2+ существенно
i
больше, чем в цитоплазме, указывается даже, что в 10 000 раз [10]. Если мембранный потенциал равен потенциалу покоя, то кальциевые каналы в мембране клетки закрыты. При деполяризации они открываются. Можно считать, что существует порог потенциала un,ca, при достижении которого значительно увеличивается проницаемость мембраны к ионам Ca2+. Будем полагать секундный поток ионов кальция равным jca = ©(u — Mn,ca)Qca, где © - функция Хевисайда, u - активный потенциал, Qca = const > 0. Источником увеличения концентрации Ca2+ в цитоплазме кардио-миоцита служит не только внешняя среда, но и находящийся в клетке третьего типа саркоплазматический ретикулум (SR). Это запас ионов кальция в емкостях, которые открываются при присоединении к ним некоторого количества тех же ионов Ca2+. Обратно ионы Ca2+ отсасываются кальциевыми насосами, т. е. ферментами, встроенными в SR и наружную мембрану кардиомиоцита. Потому уравнения баланса кальция запишем так:
ax6x9 — biX6 — b2X6 + ©(u — Un,Ca)Qca — Q акт, (4)
dt
dx- — —ax exg + bixe.
Здесь a - вероятность одной группы ионов Ca2+ покинуть в единицу времени SR, xg -число таких групп в SR, bi и b2 - характеристики насосов, QaKT - секундный расход Ca2+ на активизацию центров актиновых волокон.
4. В отличие от аксона нейрона кардиомиоциты не имеют специальных открывающихся и закрывающихся калиевых каналов, а ионы K+ проходят через мембраны клеток по постоянно действующим проходам (порам). Поэтому для математического описания активного потенциала более других подходит модель ФитцХью-Нагумо в принятом в статье [7] виде
dV ~ d2V dJ
с—— = D^—г + F — k\V + JH, —- = -k2V, k\ = const >0, k2 = const > 0, (5) dt dx2 dt
где c - удельная емкость мембраны клетки, приходящаяся на 1 см2 ее поверхности; V = u — ир есть отклонение активного потенциала от потенциала покоя (равновесия); up « —90 мВ; D - коэффициент диффузии, умноженный на c; x - продольная координата; F - суммарная плотность тока через открываемые каналы, а также выходов ионов из SR и расхода Ca2+ на активизацию протофибрилл; JH - плотность тока ионов через поры мембраны.
Поскольку клетки миокарда являются конечными в процессе передачи активного потенциала, который на них разряжается, вызывая напряжение и сокращение предсердий и желудочков, то при построении простейшей математической модели будем опускать диффузионный член в уравнении (5). При этом частная производная по t заменяется на обыкновенную производную
Математическая модель сердечного сокращения должна быть адекватной процессам, представленным на рисунке. Вверху представлен график изменения во времени относительного потенциала для конкретного значения продольной координаты, внизу - соответствующий фрагмент электрокардиограммы (ЭКГ). Были выявлены следующие этапы возбуждения миокарда левого желудочка. Участок предварительного возбуждения AB, отвечающий передаче потенциала от волокон Пуркинье. Точка B соответствует достижению порога по и открытию натриевых каналов. В точке C открываются кальциевые каналы. В точке D закрываются натриевые каналы. Точка
График изменения во времени относительного потенциала для конкретного значения продольной координаты (вверху) и соответствующий фрагмент ЭКГ (внизу)
E отвечает началу высвобождения ионов Ca2+ после ферментной реакции в результате сжатия желудочка. Точкой C отмечается закрытие кальциевых каналов. Точка F соответствует полному изгнанию ионов Ca2+ из цитоплазмы и окончанию периода абсолютной рефрактерности, или полной невозбудимости. На ЭКГ участок RST отвечает периоду сжатия желудочка.
Как показано в книге [6], механизм открытия и закрытия каналов может быть основан на свойствах трехслойной наружной мембраны клетки.
5. Анализ уравнений (1), (4), (5) будем проводить при возможных упрощениях для различных участков возбуждения миокарда желудочка. При этом сохраняем лишь главные члены и не учитываем диффузию в уравнении (5). К точке A (см. рисунок) кардиомиоцит приходит в состоянии равновесия при u(tA) = ир, V(tA) = 0. К моменту времени tA ионы Ca2+ почти полностью удалены из цитоплазмы с помощью кальциевых насосов. Как следует из расчетов, проведенных в книге [11], напряжение при расслаблении (диастоле) ничтожно мало по сравнению с давлением сжатия (в систоле). Поэтому принимаем x1(tA) = 0. Поскольку 30% объема кардиомиоцита занимают митохондрии [2], перерабатывающие АДФ в АТФ, то будем считать, что xs(tA) = 0. Уравнения равновесия (3) вместе с третьим интегралом формулы (2) и уравнения (4) дают
xi =0, x2 = Np, x3 = x4 = 0, x5 = Np, x6 = 0, Np = const, (6)
x7 = Np — s > 0, xs =0, x9 = xg^rnax = L = const > 0.
На участке АВ происходят передача через электрический синапс положительных ионов на граничащую с волокнами Пуркинье клетку и переход этого возбуждения на контактирующие с первой другие клетки желудочного миокарда. Кардиомиоциты соединены в длинные цепи с помощью вставочных дисков, которые позволяют свободно перемещаться ионам цитоплазмы. Быстрота выравнивания потенциала в смежных клетках обсуждается в работе [4].
В точке В потенциалы почти всех кардиомиоцитов желудочка достигают первого порогового значения V(гв) = = «п^а — ир > 0, при котором открываются
натриевые каналы [6] и начинается процесс автогенерации в каждой клетке миокарда желудочка. Ограничиваясь в формуле (5) учетом главного члена ^ — ^ма, получаем V (1а) — V (1в) + — ¿в).
В точке С потенциал кардиомиоцитов достигает второго порогового значения V(гс) — Уса, Уса — мп,са — «п.ш > 0. До некоторого значения времени г1 > гс основным будет поток Са2+ извне во внутрь клетки. При этом, в силу формул (4)—(6), Х6(¿1) — Qca(tl — гс), V(¿1) — V(гс) + + ЯоаШ — гс).
Вышедшие из межклеточного пространства ионы Са2+ открывают емкости БИ, из которых выходит основная часть сократительных ионов Са2+. После точки С мостики миозина последовательно входят в ферментную реакцию, моделируемую уравнениями (1). Возбуждение миокарда будет изометрическим N — Ыр, так как клапаны со стороны как предсердия, так и аорты закрыты и объем крови в желудочке остается постоянным.
На промежутке (£1 ) главные члены уравнений (4)
dx6 dxg dt dt
дают
ax6xg, Х6 + xg = L = const > 0
X6 = LCi(Ci + e-a(t-tl))-i, Ci = (Lx— (ti) - 1)-1. (7)
С помощью формул (1), (2) и (7) получаем
dx2 2
х\ = Np — Х2, —jj- = —ax2 + ax2xz, x3 = 0, Ж4 = X2 — Ж5.
Поскольку до точки величины x7 не присоединяются к xi, то принимаем в = 0. Отсюда имеем
x5 = C2(1 + CieaL(t-tl))K-1a, C2 = Np(1 + Ci)K-1a, x7 = s - Np, xs = 0.
В качестве главного члена формулы (5) принимаем F = QNa + vxc6, где v = const > 0 - некоторый коэффициент. Для активного потенциала находим
V(t) = c-iQNa(t - ti) + c-ivx6 + C3, C3 = V(ti) - vc-ix6(ti).
В точке D происходит закрытие натриевых каналов. К моменту времени t = tp содержащиеся в SR ионы Ca2+ полностью перешли в цитоплазму. Происходит бурный процесс присоединения кальция величины QaKT за единицу времени к центрам тонких протофибрилл. Для изменения активного потенциала запишем
dV
с—— = —bQaKT, Ъ = const > 0, dt
V(t) = V(tD) — c-1bQaKT(t — to), to < t < tE.
При резком спаде потенциала по-прежнему происходит изометрическое возбуждение мышечных волокон с возрастанием xi. Повышение напряжения сжатия желудочков приводит к последовательному более жесткому прижатию миокарда к эпикарду и околосердечной сумке, что четко прослеживается на фрагменте ЭКГ (см. рисунок) в виде отрезка QR. По литературным источникам [4, 9] промежуток времени tD — tA очень короткий. Вся ломаная линия ABCD (см. рисунок) изображается в виде одного почти вертикального отрезка. Мызок DE мало заметен и весьма кратковременный. Можно считать, что до момента tE в основном еще не заканчиваются ферментативные реакции, связывающие мостики миозина с актином.
6. В точке E (и отвечающей ей точки R на ЭКГ) величина xi достигает наибольшего значения. Давление в желудочках несколько выше, чем в магистральных артериях. Открываются полулунные клапаны и, преодолевая инерцию столба крови в указанных артериях, начинается быстрое изотоническое сокращение желудочков. Происходит последовательный отход миокарда от перикарда, что выражается отрезком RS на ЭКГ. Будем считать, что в точке E величины Х2 и Х5 становятся столь малыми, что ими можно пренебрегать, т. е. x2(tE) = x3(tE) = x5(tE) = 0. Тогда x1 (tE) = Np, x4(tE) = 0, Х7(tE) = s — Np, xg(tE) = 0. Все мостики, находящиеся в перекрытии миозина и актина, оказываются задействованными. Вместе с тем имеется достаточно предварительно активизированных АТФ головок миозина и ионами Ca2+ центров актина, чтобы обеспечить быстрое увеличение N(t), t ^ tE, т. е. быстрое сокращение желудочков.
Длительность быстрого сокращения статистически отвечает времени перехода молекулы АДФ в АТФ с помощью митохондрий. При этом после отхода головок миозина от центров актина должно пройти некоторое время рефрактерности, по истечении которого головки миозина и центры актина могут снова возбуждаться. Так что увеличение N сначала происходит за счет поступления в перекрытие новых центров актина и головок миозина. Поэтому на некотором коротком промежутке должно быть к = y = X = 0. При xi = Np на основании (1) получаем
d dx3 dlnN
— InN = (Зх7, x2 = xA = 0, = Хз~~¿t--\-pNpx7,
dxr dxg
хь = x3, — = -2fíNpx7, — = 2pNpx?.
Отсюда находим
x7 = (s — Np)e-2eNp(t-tE}, lnN = lnNp + Np-1(s — Np)(í — e-2eNp(t-tE >). На основании (2)
x3 = x5 = Np, xg = s — Np — x7. Ввиду кратковременности рассматриваемого интервала можно принять
N = Npe2e(s-Np)(t-tE), x7 = (s — Np)(1 — 2pNp(t — tE)).
Как отмечалось, быстрое сокращение желудочков заканчивается, когда концентрация АТФ (переменная x7) сильно уменьшается и скорость сокращения желудочков лимитируется работой митохондрий по восстановлению АТФ на основе АДФ. Начинается медленное сокращение желудочков, для которого следует рассматривать неупрощенные уравнения (1). Качественные выводы относительно хода ферментативных реакций здесь сделать трудно. При численном интегрировании потребуются конкретные
значения входящих в указанные уравнения коэффициентов. Поэтому возникает задача идентификации коэффициентов, используя экспериментальные данные.
В момент Ье начала сжатия желудочков выход Са2+ из БИ и активизация кардио-миоцитов в основном заканчиваются. В уравнениях (4) можно получить а = 0, Qaкт = 0 и для участка ЕС'
Х6 = САв-(Ъ1+Ъ2)(1-1Е ) + Qca(Ьl + Ы-1, С4 = Х6(Ье ) - Qca(bl + Ъ2)- (8)
На участке С'Е картина несколько изменяется:
Qca = 0, Х6 = С5е-(Ъ1+Ъ2)(4-4^'), С5 = Хб(Ьс') > 0. (9)
Можно считать, что точка Е (или точка Т на ЭКГ) соответствует полному удалению ионов Са2+ из цитоплазмы кардиомиоцита. В момент времени Ьр заканчивается абсолютный рефрактерный период полной новой невозбудимости кардиомиоцитов [9].
Рассмотрим форму волны активного потенциала на участке ЕЕ. Освобождение ионов Са2+ при закрытии выходов из БИ, окончании активизации центров тонких протофибрилл и ферментативных реакций между миозином и актином, внешний вход Са2+ через кальциевые каналы при закрытых натриевых каналах приводят к резкому замедлению в момент времени Ье падения активного потенциала, несмотря на продолжающийся выход ионов К+. Возникает участок ЕС' так называемого «плато».
Для качественного описания плато, или участка сравнительно медленного убывания активного потенциала, сохраним главные члены системы (5)
¿V ^ ¿х
c— = Jк + Jca, -тг =-«гЦ =
аЬ аЬ аЬ
6
С учетом (8)
Отсюда имеем
¿2л/
С71+ к2У = иСА{Ъ1 + Ы2е-(Ь1+Ь2)(^в). (Ю)
Асоъ{Щс-ъ{г -гЕ) + а) + Ве-{Ь1+Ь2)^-*Е\ В = 1^СА{к2{Ъ1 + Ь2у2 + с)"
сова = А-\У(гЕ) - В), 5та = А^Щс-^ 1(^-\г=гЕ+0 +{Ъ1+Ъ2)ВЛк ^^
V
ч аъ
На участке С'Е сохраняются формулы (10) и (11) при использовании (9) и замене Ье на Ьс', С4 на С5.
График изменения активного потенциала V(Ь) (см. рисунок) показывает, что на отрезке [Ьс>,Ьр] должно быть ^ < 0. Тоже следует из (11) при Ьр — Ьр < (тг — а)к2 2 . На рисунке видно, что на полуинтервале [ЬЕ,Ьр) имеем < 0, т. е. выпуклость вверх. Указанное неравенство равносильно получаемым из формулы (10) условиям
V (Ь) >к-1и(Ъ1 + б2)2С4в-(Ъ1+Ъ2)(-4в), Ье < Ь < Ьс,, (12)
V(Ь) > к-1и(Ъ1 + Ъ2)2С5в-(Ъ1+Ъ2)(4-4с'), Ьс, < Ь< Ьр. (13)
Соотношения (12) и (13) можно обеспечить с помощью формулы (10) за счет выбора параметров. Точка Е является точкой перегиба. Поэтому момент Ьр может быть рассчитан из трансцендентного уравнения
V(tF) = k-iv(bi + b2)2C5e-(bl+b2)(t-tC),
где левая часть определяется формулой (11) при замене t на tF, tE на tC>, C4 на C5.
Заключительная часть заднего фронта волны активного потенциала, начинающаяся с точки F, не содержит кальциевого обмена и может быть построена по одному из описанных в статье [7] варианту.
7. В настоящей статье сделана попытка рассмотреть изменение активного потенциала клеток миокарда и их мышечного сокращения с помощью единой, хотя и очень грубой математической модели. Предлагаемая модель качественно правильно отвечает основным нюансам физиологической задачи по управлению автоволновым процессом сердечного возбуждения. При этом для разных этапов сердечного цикла константы скорости мономолекулярных и бимолекулярных реакций уравнений (1) могут быть различными. Так, до момента времени начала сокращения tE отсутствует выход молекул АДФ после окончания ферментативной реакции и следует полагать в = 0. В период быстрого сжатия желудочков ввиду рефрактерности и наличия предварительно активизированных головок миозина и центров актина было принято к = y = х = 0. Когда для поддержания требуемого напряжения xi = Np привлекаются неактивизиро-ванные головки миозина и центры актина, то это приводит к замедлению сердечного цикла. Возникает период медленного сокращения. Удалось получить чисто теоретическое обоснование существования участка «плато», на котором имеется необычная для профиля нервного импульса выпуклость вверх.
Литература
1. Словарь физиологических терминов / отв. ред. О. Г. Газенко; редкол.: И. А. Агаджанянц, О. С. Андрианов, А. А. Башкиров и др. М.: Наука, 1987. 448 с.
2. Биологический энциклопедический словарь / гл. ред. М. С. Гиляров; редкол.: А. А. Баев, Г. Г. Винберг, Г. А. Заварзин и др. М.: Сов. энциклопедия, 1989. 864 с.
3. Анатомия человека. Т. 2 / авт.: Э. И. Борзяк, В. Я. Бочаров, Л. И. Волкова и др; под ред. М. Р. Сапина. М.: Медицина, 1986. 480 с.
4. Keener J., Sneyd J. Mathematical physiology // Interdisciplinary Applied Mathematic. 1998. Vol. 8. 727 p.
5. Новоселов В. С. Статистические модели нейродинамики: препринт факультета прикл. математики—процессов управления С.-Петерб. ун-та. СПб., 2004. 64 с.
6. Новоселов В. С. Статистическая динамика: учеб. пособие. СПб.: Изд-во С.-Петерб. ун-та, 2009. 393 с.
7. Новоселов В. С. К имитационному моделированию нервного импульса // Вестн. С.-Петерб. унта. Сер. 10: Прикладная математика, информатика, процессы управления. 2011. Вып. 4. С. 73—83.
8. Новоселов В. С. К математической модели пейсмекера // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2012. Вып. 4. С. 59—64.
9. Кринский В. И., Медвинский А. Б., Панфилов А. В. Эволюция автоволновых вихрей. М.: Знание, 1986. 48 с. (Сер. Математика, кибернетика, № 8).
10. Беркинблит М. Б., Жердяев А. В., Тарасова О. Г. Задачи по биологии человека и животных: экспер. учеб. пособие. М.: МИРОС, 1995. 176 с.
11. Бегун П. И., Афонин П. Н. Моделирование в биомеханике: учеб. пособие. М.: Высшая школа, 2004. 390 с.
Статья поступила в редакцию 30 мая 2013 г.
