Моделирование мышечных усилий при движении стопы, вызванном снятием внешней нагрузки

RUSSIAN JOURNAL OF BIOMECHANICS

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

(ЭМГ) – биоэлектрической активности мышц [2]. Описание мышечного сокращения в ответ на изменение потенциала действия мотонейронов в мышце представляется достаточно сложной задачей, так как она требует рассмотрения совокупности различных процессов как механических, так и биохимических.

Классическая работа Хилла [3] и ее различные модификации [4–9] не позволяют решить эту задачу в

: 0009-0001-1063-8484

: 0000-0003-2839-4817

Эта статья доступна в соответствии с условиями лицензии Creative Commons Attribution-NonCommercial 4.0 International

License (CC BY-NC 4.0)

This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0)

должной мере, так как в рамках этих моделей отсутствует четкость описания сокращения мышцы в ответ на нервные импульсы мотонейронов. λ – модель, предложенная Фельдманом [10], позволяет учитывать работу ЦНС в сокращении мышц. Однако предлагаемые подходы основаны на результатах статических экспериментов и носят феноменологический характер, связь их параметров с особенностями функционирования мышечной ткани неочевидна. Подходы, учитывающие сложную динамику актин-миозинового взаимодействия [11; 12], рассматривают сокращение на уровне одного саркомера, что не позволяет напрямую применять данные подходы к описанию сокращения целой мышцы. Более подробный обзор исследований, посвященных моделированию мышечного сокращения, представлен в работе [13].

Модели [11; 12] хоть и описывают сокращение только на уровне саркомера, имеют физиологически обоснованные параметры, полученные из экспериментальных данных [11; 12], что, несомненно, является их преимуществом. Поэтому, модель сокращения мышцы может быть основана на модели сокращения одного саркомера. Однако высокий порядок дифференциальных уравнений модели затрудняет их использование.

В работе [13] на основе [11; 12], описывающих сокращение саркомера сердечной мышцы, построена модель сокращения скелетной мышцы, использующая уравнения пониженного порядка. Она учитывает механизмы сокращения мышечных волокон и представляет мышцу как комплекс активных двигательных единиц. Входом является средняя скорость притока ионов кальция в мышечное волокно. Предполагается, что она пропорциональна потенциалу действия мотонейрона. При этом модель пониженного порядка пренебрегает описанием части быстрых переходных электрохимических процессов в мышечном волокне. Основным ограничением является требование к изменению потенциала, который должен быть «медленно» меняющейся гладкой функцией: характерные времена его изменения должны быть одного порядка с временем рассматриваемого движения. Поэтому использование данной модели сокращения оправданно в условиях гладкого тетануса.

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

Показано, что модель работает корректно, если на ее вход поступает суммарный эквивалентный сигнал для всей мышцы. В работе [14] подход [13] был применен для описания реакции мышцы на одиночное импульсное воздействие. Такой вид сокращения не типичен для функционального сокращения мышц, но он позволяет оценивать различные свойства мышечных волокон.

Целью данной работы является применение модели [13] для описания мышечных усилий в процессе реального сокращения мышц человека. Поскольку измерение усилия отдельно взятой мышцы в процессе реального движения затруднительно, в работе обсуждается описание движения стопы, которое вызвано мгновенным снятием нагрузки. Вместо измерения усилия мышц голени напрямую с помощью методов видеоанализа оценивается движение стопы. В вызванном движении участвуют мышцы-антагонисты голени, биоэлектрическая активность которых может быть измерена в процессе движения. На основе ЭМГ мышц голени формируются входные сигналы для модели [13]. Неизвестными параметрами выступают коэффициенты линейной связи между ЭМГ и входом в модель сокращения [13]. Для определения недостающих параметров модели проводится серия исследований, в рамках которых регистрируется ЭМГ мышц-антагонистов голени, и эта информация используется для «восстановления» движения стопы. Найденные неизвестные параметры далее используются для оценки состоятельности модели в описании других подобных движений стопы.

Материалы и методы

Обследование проводилось следующим образом. Испытуемый располагался на кушетке в положении лежа на животе. Стопа c жестко прикрепленной к ней подошвенной платформой свободно свисала и могла совершать вращение в голеностопном суставе. С помощью нити стопу оттягивали на угол порядка 20о; силу натяжения измеряли динамометром. В некоторый момент времени нить перерезали, что приводило к свободному движению стопы и изменению угла в голеностопном суставе.

Для оценки изменения угла использовалась система видеоанализа Qualisys . Маркеры размещались на стопе и дополнительной платформе (рис. 1). Одновременно на мышцы голени – камбаловидную, латеральную и медиальную головки икроножной мышцы, а также переднюю большеберцовую мышцу – накладывались поверхностные электроды для регистрации ЭМГ. Процедура с обрезанием нити повторялась несколько раз для одного испытуемого.

Для дальнейшей обработки ЭМГ-сигнал сглаживался с использованием окна Ханна шириной 150 мс. На рис. 2 и 3 приведены примеры изменения угла в голеностопном суставе по данным видеоанализа и соответствующие ЭМГ мышц голени.

В данном примере амплитуда движения стопы составляет приблизительно 15˚, а характерное время движения составляет приблизительно 0,3 с, что позволяет применить модель 1-го порядка [14] для описания процесса сокращения мышц в данном движении.

Рис. 1. Система «платформа-стопа» и маркеры системы видеоанализа, прикрепленные к ней.

Черные сферы – маркеры видеоанализа.

Синяя линяя – натянутая нить

Время, с

Рис. 2. Изменения угла в голеностопном суставе по данным видеоанализа

Время, с

Рис. 4. Система «платформа-стопа». O – центр голеностопного сустава. S – центр масс системы стопа-платформа. Ось Ox направлена вдоль ноги. Ось Oy направлена вдоль силы тяжести. F _1, F _2 – силы, создаваемые передней и задней группой мышц голени, H _1, H _2 – плечи данных сил

Рис. 3. Миограммы для передней большеберцовой (голубая линяя), латеральной икроножной (оранжевая линяя), медиальной икроножной (желтая линяя), камбаловидной мышц (фиолетовая линяя)

Математическая модель движения стопы

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

Систему стопа-платформа будем описывать физическим маятником с массой m и моментом инерции относительно центра голеностопного сустава J . Голень считаем неподвижной, а голеностопный сустав моделируем цилиндрическим шарниром. К стопе прикреплена пара мышц антагонистов голени. Центр масс системы S находится на расстоянии is = |OS| от центра голеностопного сустава. Тогда уравнение движения этой системы запишем следующим образом:

I          п )

J s V 1 = m s gl s sin I V 1 - - I ⁺ M S — k v V 1 ,

Ms = cos p2F2H2 - cos pFH, где угол ^ соответствует углу в голеностопном суставе, p соответствует углу перистости мышц-антагонистов голени, которые определены в [15], M – момент сил, создаваемый мышцами-антагонистами в голеностопном суставе, F, F – это силы, создаваемые передней и задней группой мышц голени соответственно, которые будем определять с помощью подхода [16], H, H – это плечи соответствующих сил, которые определены в [16], k – это параметр вязкости в голеностопном суставе.

Для полного описания модели (1) требуется определить параметры системы ms, Js, ls . Масса системы ms складывается из массы стопы mt и платформы m2.

Масса платформы  m2 = 0,45. Масса стопы оценивается    для    каждого    человека    по антропометрическим данным согласно зависимости, представленной в работе [16]. Момент инерции системы относительно центра голеностопного сустава складывается из моментов инерции стопы J и платформы J2 относительно центра голеностопного сустава Js = J + J2.

Оценка момента инерции платформы   J₂

относительно центра голеностопного сустава определяется по характерному времени свободных колебаний      платформы      и      составляет

0,0354 кг^м2. Оценка момента инерции стопы относительно центра голеностопного сустава приведена в работе [17]. Для оценки расстояния центра масс от центра голеностопного сустава l требуется оценить  dx  и d2  - радиус-векторы от центра голеностопного сустава до центров масс платформы и стопы соответственно:

,   m,d, + md ds =              , ms                       (2)

l s = d s |-

W (x ) =

Ф(5) = <

G (5) =

00 (5) =

1 - k,

1 — к,2

f ^1,  ₅2

*

l

' S 0 ,1

^l s 0 ,2

(5.-5)2,

- b + с 5 ²,

‘ I ⁵

b +--,

[   5.-8

с

' 1 + e^yS ’

a (x) = (1 + knx) 1,

Х( c )= 1 +

У

^где F pi

s 0 ,2

Bk B

(C + кв )2 J ’

пассивная сила

i-й мышцы, обусловленная

упругостью мышечного комплекса ( i = 1 соответствует передней группе мышц, а i = 2 - задней), к_г - пассивная

жесткость i-й мышцы, L_; - длина i-й мышцы, определенная в [18], L_Oi - длина i-й мышцы в нерастянутом состоянии, F_a. - активная сила i-й

Координаты d оцениваются по данным из работы [17]. Координаты d оцениваются по координатам маркеров видеоанализа для платформы (маркеры 4, 5, 6). Усилие, создаваемое мышцами-антагонистами голени, будем описывать моделью первого порядка из [13; 14]:

F i = F p , i + F a , i ,

F ,i= k i ( L i - L ,o ) ,

F a , , = S i EN xb N m hW | T L - F +6 i ) n , y¹ ’ sarc J

ni   ф(5 i) + G (5,-) ’

0 = — v_si A .    s ’ '

1 5 i G ( 5 ) ; V s ,i

d dt

0i =0o (5i),

с, (1 - Ai)m - a(п,0i)(Ai)m = 0, Г dc =       Qm    ci- c*)

dt x(ci)[ ‘        (c + kx)

мышцы, создаваемая актин-миозиновым комплексом, S; - площадь «поперечного сечения» i-й мышцы (для антагонистов голени существуют характерные значения поперечников [16]), N - число саркомеров в поперечнике, E - жесткость миозиновых мостиков, Nxb - количество головок миозина, приходящихся на половину длины миозиновой нити, N - количество миозиновых нитей на единицу поперечного сечения саркомера, n, - средняя доля миозиновых мостиков в i-й мышце, присоединенных к актину, 5 i - среднее по ансамблю безразмерное смещение миозиновых мостиков после их присоединении к актину в i-ой мышце, 0i - доля мостиков, находящихся в состоянии генерации силы в i-й мышце, h - перемещение головки миозина в состояние генерации силы, w (x) - безразмерная длина зоны перекрытия в половине длины саркомера, k , k , l , l     - параметры кусочно-линейной зависимости длины зоны перекрытия от длины саркомера, 5 - максимальное значение 5 во время растяжения, b, c, у, A, - константы, Ai - средняя вероятность присоединения мостиков в i-й мышце, с, - средняя безразмерная концентрация кальция в i-ой мышце, m и кп - параметр аппроксимации Хилла

зависимости «сила - концентрация кальция» и степень зависимости связывания Ca ²⁺ c тропонином от количества присоединенных к актину миозиновых головок соответственно, B - полная безразмерная концентрация кальциевого буфера в мышечных клетках, k_B - нормированная постоянная равновесия буфера, X - скорость откачки кальция в саркоплазматический ретикулум, с. - нормированная концентрация кальция в расслабленной мышце, а₊ - скорость перехода тропомиозина в открытое состояние после того, как тропонин связался с кальцием, k_x - константа, Q. ( t ) - средняя скорость притока ионов кальция в i-й мышце, N_sarc - число саркомеров вдоль мышцы. Значения параметров описывают работу актин-миозинового комплекса, считаются одинаковыми для всех скелетных мышц и приведены в работе [13].

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

Второе уравнение системы (3) описывает линейную упругую пассивную силу мышцы. Третье уравнение системы (3) задает связь между активной силой мышцы и процессом актин-миозинового взаимодействия. Активное усилие пропорционально доле присоединенных мостиков п, , безразмерному смещению мостиков 6 , и 0_; доли мостиков, находящихся в состоянии генерации силы. Уравнение 4 системы (3) описывает динамику присоединения и отсоединения миозиновых мостиков [19]. Уравнение 5 системы (3) описывает деформацию миозиновых мостиков в результате их присоединения к актину и смещения [20].

Уравнение 6 системы (3) описывает динамику изменения числа миозиновых мостиков, которые после присоединения перешли в состояние генерации силы [19]. Уравнение 7 системы (3) регулирует взаимодействие между белками актина и миозина в ответ на изменение концентрации ионов кальция [21]. Уравнение 8 системы (3) представляет собой простейшее уравнение баланса ионов кальция в мышечных клетках [22]. Обоснование перехода от системы дифференциальных уравнений, описывающих сокращение саркомера, к одному дифференциальному уравнению и набору алгебраических уравнений, описывающих сокращение мышцы для времен сокращения свыше 0,5 с., приведено в [13]. Все параметры модели определены на основе различных экспериментальных данных [11; 12], за исключением величины Q которая является входом для модели.

Модель движения (1), (3) описывает движение стопы в ответ на изменение средней концентрации кальция в мышцах голени. Поскольку концентрация кальция в мышце увеличивается в ответ на изменение потенциала действия на мотонейронах мышцы, между потенциалом действия и концентрацией кальция имеется сложная связь. Миограмма отражает изменение потенциала действия в группе двигательных единиц, потенциал которых отражают измерения миографа. Обозначим V потенциал из миограммы для i-й мышцы, тогда скорость притока ионов кальция в i-й мышце можем оценить:

Q i = rV                  (4)

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

Определение неизвестных параметров модели для оценки угла в голеностопном суставе по ЭМГ мышц голени

Применим модель движения стопы (1), (3), (4) для определения неизвестных параметров модели и последующего восстановления угла в голеностопном суставе по данным ЭМГ мышц-антагонистов голени.

Обозначим Q_va ( t ) угол в голеностопном суставе (угол между осью Оу и радиус-вектором OS), восстановленный по координатам маркеров видеоанализа на участке движения стопы [ 0, T ] , где T обозначает момент времени, когда стопа пришла в положение равновесия, а в качестве начальной точки по времени возьмем момент обрыва нити.

В данной задаче будем восстанавливать относительное изменение угла в голеностопном суставе, то есть нас будет интересовать изменение угла в  суставе. Поэтому, не нарушая общности рассуждений, будем считать, что 0^ (0) = 0 .

До обрыва нити мышцы голени находятся в напряженном состоянии, так как на стопу действует сила натяжения нити. Это создает дополнительный момент в голеностопном суставе. Уравнение равновесия системы запишем в виде:

m s gl s sin ^¥ 1 ( 0 ) - ПJ ⁺ M S ( ^ 1 ( 0 ) ) = M T ,

M S = cos P 2 F 2 H 2 ( ¥ 1 ( 0 ) ) - cos P 1 F 1 H 1 ( ¥ 1 ( 0 ) ) ,  (5)

F = F p ,. + F a ,i ,

F p , i = k i ( L i ( ¥ . ( 0 ) ) - L i ,0 ) ,

T ₀                                   (9)

Smin = min S ( X ) , где S(X) описывает среднеквадратичную ошибку оценки. Задача минимизации (9) решалась при ограничениях на вектор переменных:

где M обозначает момент сил, вызванный натяжением нити. В конечный момент времени T стопа будет свободно свисать, и тогда уравнение равновесия будет иметь простой вид:

m s gl s sin ^¥ 1 ( T ) ^- |j ⁺ M S = 0,

M S = cos P 2 F 2 ^H 2 ( ¥ 1 ( T ) ) ^{- cos} P 1 ^F 1 H 1 ( ¥ 1 ( T ) ) ,  (6)

F i = F p ,i + F_a , i ,

F p , i = k i ( L i ( ¥ 1 ( T ) ) - L. 0 ) ,

Уравнения (5), (6) используем для нахождения пассивных     жесткостей     мышц-антагонистов голени k, k. Для этого требуется определить начальное и конечное положение стопы ¥, (0), ¥i (T):

¥1 (0) = П + !9va (T)l ,

²                             (7)

Соотношения (5)–(7)  позволяют определить параметры пассивных жесткостей мышцы k, k.

Обозначим Q_m ( t ) оценку угла в голеностопном суставе, полученную с использованием миограмм мышц голени. По построению модели 1, 3 данный угол будет зависеть от набора неизвестных параметров:

0 m ( t ) = 0 m ( t , r 1 , r 2 , ^k v , J S ) .                  ⁽⁸⁾

В набор искомых параметров X = ( t , r , r ₂, k_v , j_s ) добавлен момент инерции системы, поскольку его предварительная оценка неточна и заметно влияет на описание динамических процессов, а также параметр вязкости в голеностопном суставе, поскольку напрямую его измерить не представляется возможным, а различные его оценки из литературы варьируются достаточно сильно в зависимости от человека и движения, которое он совершает [23–25]. Для нахождения искомых параметров поставим задачу минимизации отклонения оценки угла в голеностопном суставе, полученной с использованием ЭМГ мышц голени и значениями этих углов из данных видеоанализа:

10 4 <  r <  10 7 ,

0,5 <  k_v <  3,

1 сВ ^, Нмс,

2 ⁽ J 1 ⁺ J 2 ) <  J S <  ^{3 (} J 1 ⁺ J 2 ) .

Задача (9), (10) решалась в MATLAB методом fmincon [25–27].

Результаты

С целью анализа применимости приближенной модели [13] проведено экспериментальное обследование пяти испытуемых, для каждого из которых регистрировалась серия из 3–5 движений стопы, выполненных по описанной выше методике.

Так как набор соответствует одному человеку и крепление электродов миографа при этом не изменялось, для каждого набора движений определялось единственное значение вектора искомых параметров X и величины жесткостей k , k . Для каждого набора данных выбиралась самая «активная» мышца из трех задних мышц голени, и ее миограмма использовалась для описания сокращения задней «мышцы» модели, так как интенсивности других миограмм существенно ниже.

На рис. 5 приведен пример сравнения двух оценок углов: по миограммам мышц голени, полученным в результате решения модельных уравнений, и по данным видеоанализа. Как видно из графика, использование упрощенной модели позволяет приближенно восстановить ориентацию движения стопы и скорость ее изменения. В табл. 1 приведены значения S для каждого из 5 наборов движений и значения искомых параметров, A J_s - отклонение найденного значения момента инерции от аналитического, 5 S_min - отношение S_m к максимальной амплитуде изменения угла. Среднее квадратическое значение ошибки определения угла составило 11–12 %. Коэффициенты усиления миограмм получились примерно одинаковые по порядку для каждого набора движений, т.е. для каждого человека, оценки момента инерции системы отличались от априорного значения J в среднем на 16 %. Параметры вязкости голеностопного сустава, найденные в результате решения задачи минимизации,

Таблица 1

СКО и значения искомых параметров для всех наборов движений

Набор движений	u min ,	8 S min , %	r 1 , 105 · с-1В-1	r 2 , 105 · с-1В-1	k , Нмс	A Js , %
I	1,2	8,1	50,1	2,51	1,06	10
II	2,2	12,0	12,6	3,16	1,10	0
III	3,4	11,8	31,6	1,26	0,74	47
IV	1,8	11,9	31,6	2,00	1,10	5
V	2,3	12,9	25,1	4,00	0,92	-20

Таблица 2

СКО для всех наборов движений в случае использования половины серии движений для нахождения неизвестных параметров

Набор движений	\ . 0 u min ,	8 S min , %	r 1 , 105 · с-1В-1	r 2 , 105 · с-1В-1	k , Нмс	AJs , %
I	1,3	10,0	45,8	2,61	1,06	10
II	2,6	13,4	13,2	3,32	1,10	0
III	4,0	14,9	62,7	1,29	0,67	56
IV	2,0	11,7	29,1	1,87	1,10	5
V	3,8	20,8	28,0	4,01	0,92	-20

укладываются в диапазон значений, представленный в литературе для различных движений человека [23–25].

В моделировании выше показано применение для восстановления траектории движения стопы по миограммам мышц, когда вектор неизвестных параметров X находился для всех движений. Можно рассмотреть задачу нахождения вектора неизвестных параметров X по ограниченному набору движений, т.е. из всей серии движений использовать, например, только половину движений для нахождения Xmin = arg min S ( X ) , а для второй половины X восстанавливать движение по найденному решению

^X min

.

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

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

Обсуждение

Проведенное обследование показывает принципиальную возможность восстановления движения стопы по ЭМГ мышц голени. Опосредованно это подтверждает возможность восстановления мышечных усилий по ЭМГ с помощью модели [13] для

Рис. 5. Пример восстановления угла в голеностопном суставе по ЭМГ (синяя сплошная линия) в сравнении с углом в голеностопном суставе, восстановленным из данных видеоанализа (красная прерывистая линия)

реальных движений. Преимуществом данного подхода является то, что для каждой мышцы требуется определить всего один параметр связи между ЭМГ и входом модели [13], чтобы по ЭМГ восстанавливать усилие мышцы в процессе движения. В литературе существуют различные подходы к восстановлению мышечных усилий по ЭМГ [2], которые основаны либо на моделях Хилла [28], либо на современных методах анализа данных с привлечением алгоритмов машинного обучения [29–31]. Первые подходы требуют определения различного рода физиологических параметров и начальных условий, вторые подходы требовательны к объему данных для тренировки и валидации. Модель [13] допускает проведение небольшой серии обследований для определения связи между ЭМГ и средней скоростью притока ионов кальция в мышце.

Заключение

В работе использовалась упрощенная модель сокращения скелетной мышцы [13] для описания движения стопы, вызванного мгновенным снятием внешней нагрузки. Для движения стопы было проведено экспериментальное обследование, в ходе которого регистрировались электромиограммы мышц голени и угол в голеностопном суставе. Решение задачи минимизации невязки между изменением угла в голеностопном суставе, восстановленным по данным ЭМГ с использованием приближенной модели мышцы, и результатами видеоанализа позволило определить параметры модели, а также оценить мышечные усилия и угол в голеностопном суставе по данным электрической активности мышц. Средняя ошибка восстановления угла составила 11–12 %.

Отметим, что предлагаемая упрощенная модель отличается от широкоизвестных упругих моделей,