Математическое моделирование движения, управляемого мышцами, действующими в противоположных направлениях Текст научной статьи по специальности «Математика»

УДК 517.958:57+531/534:57

К. А. Кликунова, В. П. Трегубое

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ, УПРАВЛЯЕМОГО МЫШЦАМИ, ДЕЙСТВУЮЩИМИ В ПРОТИВОПОЛОЖНЫХ НАПРАВЛЕНИЯХ

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

А Б

Рис. 1. Схема работы миозинового мостика при изменении длины мышцы (А)

и ее укорочении (Б)

Y - число свободных миозиновых мостиков; n, m, и p - число тянущих, тормозящих и супертянущих мостиков в слое мышцы длиной в половину саркомера и площадью поперечного сечения 1 см2; ki, , кз - скорости перехода мостиков из одного состояния в другое; v - скорость сокращения саркомера; 5 - смещение головки мостика в тянущем состоянии.

Кроме того, в отличие от лабораторных условий, моделированию которых авторами предшествующих работ уделялось основное внимание (см., например, [2-6]), необходимо было учесть зависимость общего числа миозиновых мостиков от степени растянутости сократительной компоненты мышцы. В результате существенным образом была усовершенствована модель, разработанная в работе [7]. Хорошо известно, однако,

Кликунова Ксения Алексеевна - аспирант кафедры моделирования электромеханических и компьютерных систем факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Научный руководитель: доктор физико-математических наук,

проф. В. П. Трегубов. Количество опубликованных работ: 6. Научное направление: биомеханика. E-mail: kliksa@gmail.com.

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

© К. А. Кликунова, В. П. Трегубов, 2010

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

Механическая модель мышцы. В качестве модели скелетной мышцы использовалась механическая система из четырех деформируемых элементов (рис. 2). Модель состоит из сократительной компоненты, поведение которой моделируется, как было описано ранее [1], а также параллельной упругой компоненты и задемпфированной последовательной упругой компоненты (рис. 2). Параллельная упругая компонента обладает жесткостью с1. Ее физическим эквивалентом является упругая оболочка мышцы. Задемпфированная последовательная упругая компонента, представляющая сухожилие, обладает жесткостью с2 и коэффициентом демпфирования Ь. Инерционный элемент М подразумевает под собой массу мышцы вместе с массой прикрепленного к мышце груза.

Рис. 2. Механическая модель мышцы

I - длина мышцы; 1\ - длина сократительной и параллельной упругой компоненты;

12 - длина последовательной упругой компоненты.

При моделировании сократительной компоненты нужно иметь в виду, что при некотором изменении ее длины развиваемое усилие не зависит от степени растяжения, за исключением случая предварительно растянутой мышцы. Это было продемонстрировано в работе [1].

При построении механической модели были учтены два важных свойства мышцы. Во-первых, мышца не оказывает сопротивления на сжатие, т. е. при 11 < 101 С1 = 0 и при 12 < 102 с2 = 0, где 101,102 - соответствующие размеры покоящейся мышцы 10 = 101 + 1о2. Во-вторых, пассивное сопротивление растяжению, которое моделируется параллельной упругой компонентой, начинается лишь после 30%-ного ее растяжения. В результате с1 остается равным 0 при 11 < 1.31о1. При записи уравнений, описывающих

механическое поведение предложенной модели при наличии внешней силы Г, в условиях постоянного максимального возбуждения мышцы, переходим от длин компонент к укорочению контактильной компоненты Х1 = Д/1 и укорочению всей мышцы х = Д/. Таким образом, эти уравнения имеют вид

VI >

йп

сМ

йт

сМ

йр

сМ

йм

Л

йм1

сМ

йх1

Л

йх

Л

0:

М1

к^ао/ЗБх! - п- то - р) - ——-п,

21\д

VI

---- = ------п — кот,

" 2 N6 ’

—кзР,

п1 (п — т + 2р) — С1Х1 — пБГ

(1)

М

йм йм2

М1,

V.

где

М1 < 0:

^ = к1(а0(33х1 -п-т-р) + аъ 21\д

йт

-Ж = -к2,п-

йр

сМ

М1

п — к3р,

2Nд

йм п1 (п — т + 2р) — С1Х1 — пБГ

йъ

йм1 йм йм2

М

йъ йъ йхл

йъ

йх

Л

йъ

= М1,

= V.

М2

r]f(n - т + 2р) + (с-1 + С2)ж1 - с2ж Ъ '

(2)

В системах (1), (2) ао = 10 - общее число мостиков в слое мышцы длиной в поло-

-1-е

см,

вину саркомера и площадью поперечного сечения 1 см2, / = 3 • 10 дин, д = 10 к1 = 50с-1 - константа скорости перехода мостиков из свободного в тянущее состояние, 2N - число саркомеров в волокне, М - масса мышцы и присоединенной нагрузки, Г - внешняя сила на 1 см 2 поперечного сечения, Б - площадь поперечного сечения волокна, п - число волокон в мышце, в = 1.4 • 104 см-1, х = Д/, Х1 = Д/1, М1 - скорость сокращения конкрактильной компоненты, V - скорость сокращения целой мышцы. Начальными данными для системы (1) является состояние покоя х1 =0, х = 0, М1 = 0,

V = 0 при максимальном возбуждении мышцы п = 1013, т = 0, р = 0, т. е. все мостики в начальный момент времени тянущие. При смене направления сокращения, т. е. при смене знака скорости, переходим в решению системы (2). Начальными данными для системы (2) будут конечные значения системы (1).

На рис. 3, 4 представлены результаты работы модели соответственно при постоянной и циклической внешней нагрузке, а на рис. 5 - результат эксперимента, проведенного на изолированной мышце [3].

Рис. 3. Результаты численного эксперимента для внешней силы Е = 0, 75Р0 (Р0 = 3 • 106 дин - максимальная сила, вырабатываемая мышцей)

Р/Р0

0 50 100 150 200 250 300 350 400 450 500

І, мс

Рис. 4. Результаты численного эксперимента для циклической внешней силы 1 - сила, вырабатываемая мышцей; 2 - внешняя сила Е = 0, 5Ро (вт24£ + 1).

Р/Ро 1.00 г

0.90

t, мс

Рис. 5. Результат эксперимента, проведенного на изолированной мышце [3]

Модель управления мышцей. Как уже отмечалось, уравнения (1) и (2) описывают поведение мышцы при максимальном ее возбуждении, которое в лабораторных условиях вызывается переменным током. Для того чтобы воспроизводить движения, сходные с реальными, необходимо построить модель управления со стороны нервной системы. Известно, что мышца представляет собой совокупность так называемых двигательных единиц [9]. Двигательная единица является группой мышечных волокон, возбуждаемых одним нервным окончанием. Волокна, принадлежащие одной двигательной единице, распределены по всему поперечному сечению мышцы. Уровень возбуждения одного мышечного волокна зависит от частоты следования нервных импульсов и, которая может изменяться в пределах от 20 до 60 Гц. При выполнении циклических установившихся движений эта частота остается практически постоянной.

Другим параметром управления уровнем возбуждения мышцы является амплитуда нервных импульсов А. В соответствии с принципом Хейнемана [10], чем крупнее двигательная единица, т. е. чем больше мышечных волокон в нее входит, тем больше пороговое значение амплитуды нервных импульсов нужно преодолеть, чтобы эта двигательная единица включилась и приняла участие в сокращении мышцы. В результате параметр А определяет число двигательных единиц, вовлеченных в процесс сокращения мышцы. Поскольку амплитуда А возрастает от нуля до заданного нервной системой значения в течение некоторого времени, то и сила, развиваемая мышцей, увеличивается аналогичным образом.

Итак, управление скелетной мышцей осуществляется с помощью двух параметров: частоты следования нервных импульсов и и амплитуды нервных импульсов А, которая меняется во времени. На уровне микроструктуры мышцы можно считать, что частота следования нервных импульсов определяет количество миозиновых мостиков, которые могут принимать участие в сокращении.

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

Кроме того, введем в рассмотрение функцию Л(£), воспроизводящую управление числом двигательных единиц в мышце со стороны нервной системы с помощью изменения амплитуды нервных импульсов. Поскольку от амплитуды зависит число задействованных двигательных единиц, то с ее ростом будет увеличиваться сила, вырабатываемая мышцей. Это приведет к следующим изменениям в системе уравнений (1):

dn ав^Б VI

л=ы^Х1-п-т-р)-Ш"'

dv Л(Ь)ц/(п — т + 2р) — с\х\ — Р

= М ’

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

Рис. 6. Модели лучезапястного (А) и локтевого (Б) суставов

Для каждой мышцы используем описанную выше механическую и математическую модели. При сокращении одной из мышц внешней силой для нее является сила, вырабатываемая мышцей-антагонистом, однако, исходя из экспериментальных данных [11], мышцы-антагонисты работают поочередно, т. е. вторая мышца не активна, пока сокращается первая. К тому же при колебаниях в суставе с угловой амплитудой порядка 10° мышца, работающая как антагонист, не оказывает сопротивления на растяжение. Таким образом, уравнение движения механической модели сустава имеет вид

72

I= Щг]А(г)(п - т + 2р) - С1Х1).

Поведение второй мышцы описывается полностью аналогичной системой дифференциальных уравнений.

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

Рис. 7. Схематичное изображение электромиограммы системы мышц-антагонистов (1-11) при периодических качаниях в лучезапястном суставе человека.

Рис. 8. Результаты расчета колебательных движений модели лучезапястного сустава с моментом инерции I = 3 кг-м2 (функция сигнала Л(£) = вт30£, т = 0.12 с, а = 0.05 с)

В результате сила, развиваемая каждой мышцей, будет зависить от двух параметров: частоты и, которая остается постоянной в пределах одного движения, и амплитуды сигнала нервного импульса A(t). Для простоты можно считать, что сигнал задается функцией A(t) = A sin ut.

Как параметры управления мышцами-антагонистами, участвующими в колебательном движении сустава, как показывает эксперимент [10], целесообразно рассматривать интервал активности каждой мышцы т и интервал взаимного «молчания» - а. Эти параметры служат в нашей модели для определения моментов переключения с одной системы уравнений на другую. Изменение параметров т и а позволило моделировать колебательные движения в суставе с различными частотами колебаний.

В качестве примера на рис. 8 приведены результаты численного эксперимента для колебаний модели лучезапястного сустава с частотой 5 Гц. Для сравнения на рис. 9 показаны экспериментальные записи электромиограммы двух мышц-антагонистов. Начальными данными является состояние покоя (все мостики в разомкнутом состоянии n = m = p = 0 и v = v\ =0) растянутой на 0.1 см первой мышцы (х = 0.1 см, х\ = 0.05 см). После окончания активности первой мышцы в работу вступает мышца-антагонист с такими же начальными данными.

Рис. 9. Экспериментальные кривые для периодических качаний в лучезапястном суставе с угловой амплитудой 10о и частотой порядка 5 Гц [11]

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

Литература

1. Кликунова К. А., Трегубое В. П. Математическое моделирование переходных режимов мышечного сокращения // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2008. Вып. 3. С. 56-62.

2. Дещеревский В. И. Математические модели мышечного сокращения. М.: Наука, 1977. 160 с.

3. Deshcherevsky V. I. A kinetic theory of sriated muscle contraction // Biorheology. 1971. Vol. 7. P. 147-170.

4. Huxley A. F. Muscle structure and theories of contraction // Progr. Biophys. and Biophys. Chem. 1957. Vol. 7. P. 255-318.

5. Murase M. The dynamics of cellular motility. New York: Jonh Wiley & Sons, 1992. 358 p.

6. Zahalak G. I. A distribution-moment approximation for kinetic theories of muscular contraction // Mathem. Biosc. 1981. Vol. 52. P. 89-114.

7. Tregoubov V. P. Development of the muscle kinetic theory and cyclic movements // Acta of Bioengineering and Biomechanics. 2003. Vol. 5. Suppl. 1. P. 512-519.

8. Хилл А. Механика мышечного сокращения. Старые и новые опыты / пер. с англ. Ю. А. Шаронова; под ред. Г. М. Франка. М.: Мир, 1972. 184 с.

9. Козаров Д., Шапков Ю. Т. Двигательные единицы скелетных мышц человека. Л.: Наука, 1983. 252 с.

10. Henneman E. Functional organization of motoneuron pools: the size principle // Proc. Intern. Congr. Paris. 1977. Vol. 12. Р. 50.

11. Айзерман М. А., Андреева Е. А., Кандель Э. И., Тененбаум Л. А. Механизмы управления мышечной активностью: Норма и патология. М.: Наука, 1974. 165 с.

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

Статья принята к печати 1 апреля 2010 г.