Моделирование сократительной активности мышечной ткани Текст научной статьи по специальности «Математика»

УДК 531/534: [57+61]

Российский Журнал

www.biomech.ru

МОДЕЛИРОВАНИЕ СОКРАТИТЕЛЬНОЙ АКТИВНОСТИ МЫШЕЧНОЙ ТКАНИ

Д.А. Черноус, С.В. Шилько

Институт механики металлополимерных систем им. В.А. Белого Национальной академии наук Беларуси, Беларусь, 246050, Гомель, ул. Кирова, 32а, e-mail: shilko_mpri@mail.ru

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

Ключевые слова: мышечная ткань, сократительная активность, нити саркомера, изотоническое и изометрическое нагружение, феноменологическая модель.

Введение: О проблеме адекватного описания сократительной активности мышц

Так как эффективность функционирования различных органов во многом определяется работой мышечных тканей, составляющих около 50% массы тела человека, адекватное описание сократительной активности мышц является актуальной задачей биомеханики. Практическое значение подобных исследований заключается в том, что модели деформирования биотканей начинают использоваться при оптимизации методик лечения, реабилитации и спортивных тренировок, позволяя исключить повреждения, повысить тонус мышц, снизить утомляемость и улучшить координацию движений. Актуальность моделирования сократительной активности мышечной ткани обусловлена также попытками создания искусственных мышц на основе активных и адаптивных материалов [1] с использованием различных физических принципов генерации усилий, включая магнито- и электрострикцию, влагопоглощение полимерными гелями, эффект «памяти формы» на основе фазовых переходов и т.п.

Ранее авторами [2] была предложена феноменологическая модель мышцы, состоящая из последовательно соединенных пружины, сократительного элемента и элемента Фойгта (рис. 1). Для характеризации сократительного элемента использовано эмпирическое уравнение Хилла [3] и изложена методика идентификации параметров модели по экспериментально определяемым зависимостям длины мышцы от времени при изотоническом сокращении. Однако следует признать, что предложенное в [2] последовательное соединение сократительного и вязкоупругих элементов не вполне

© Черноус Д.А., Шилько С.В., 2006

09806267

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

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

Описание модели

Вышеназванные недостатки феноменологического описания [2] устранены в модели, состоящей из двух параллельно соединенных частей (рис. 2). Первая (верхняя по схеме) часть образована вязким (^), инерционным (к\) и двумя упругими (с, е\) элементами. Эта часть модели описывает деформирование мышечной ткани в пассивном состоянии, когда она ведет себя подобно материалу высокомолекулярной природы с вязкоупругими свойствами в широком диапазоне деформаций. Поэтому для моделирования данного состояния можно применить модель, предложенную в [5] для полимерных материалов. Отличие правой части феноменологической модели мышцы (рис. 2) от модели [5] заключается в отсутствии второго «незапараллеленного» демпфера.

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

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

Р - Р

и, =-ЪР^-. (1)

Р + а

Здесь Р0 - максимальная сила, развиваемая мышцей при изометрическом возбуждении; а, Ь - экспериментально определяемые параметры; Ыц - удлинение сократительного элемента; Ря - растягивающая сила, действующая на сократительный элемент. Точка над символом означает производную по времени от соответствующей величины.

ч

Рис. 1. Феноменологическая модель мышцы с «последовательным» сократительным

элементом

ч

Рис. 2. Феноменологическая модель мышцы с «параллельным» сократительным элементом

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

Наиболее распространенное и проработанное к настоящему времени описание взаимодействия нитей саркомера основано на гипотезе «скользящих нитей» [6, 7]. В простейшей постановке такое описание сводится к системе дифференциальных уравнений для числа тянущих и тормозящих поперечных мостиков

V V

П = к^(а- п - т) — п; т = — п - к2т . (2)

5 5

Здесь V - скорость относительного скольжения нитей; а - общее число мостиков на половине длины саркомера в области переплетения нитей; 5 - расстояние между положениями равновесия тянущего и тормозящего мостика; к\, к2 - временные константы образования и разрыва мостика, соответственно. Многочисленные уточнения и усовершенствования данной модели приводят к существенному усложнению математического аппарата и увеличению числа исходных параметров [7, 8], но не обеспечивают приемлемого соответствия экспериментальным данным [9] и затрудняют переход от анализа отдельного саркомера как элементарной

сократительной единицы мышцы к анализу мышечной ткани. Далее для характеризации сократительного элемента в рамках предлагаемой феноменологической модели (рис. 2) будем использовать дифференциальные уравнения (2).

Если мышечное волокно содержит N саркомеров, скорость относительного скольжения V выражается через скорость удлинения сократительного элемента

и

феноменологической модели V =-------—. Число N определяется как отношение длины

2 N

мышечного волокна —0 к длине саркомера в естественном состоянии /0: N = —. Для

10

различных мышц длина саркомера 10 различается незначительно и составляет около 2,2 мкм [4].

С учетом деформации нитей саркомера скорость удлинения сократительного

элемента связана со скоростью удлинения мышцы следующим образом: й5 = и - Р,

С0

где Р - сила, действующая на упругий элемент с жесткостью с0. Усилие, развиваемое сократительным элементом, является результатом действия тянущих и тормозящих мостиков: Р5 = 2N(п - т) , где/- сила, генерируемая каждым тянущим мостиком.

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

Р = Р5 + ки5 = 2 Ы/ (п - т) + к

(.. Р^

и-------

Сп

. Подставив выражение для скорости V в систему

'0 J

(2), получим дифференциальные уравнения для силы Р и чисел п, т

П = к1 (а - п - т) +

п

Г рЛ

2 N5

и — •

Со у

п

т = —к2 т —

( р\

2 N5

и—

Со у

(3)

Р = — [2Ы/ (п — т) — Р + ки\. к

Общее число а поперечных мостиков, попадающих в «рабочую» область перекрытия актиновой и миозиновой нитей, зависит от текущей длины мышцы, а следовательно - и от удлинения и. Для исходной длины мышцы (и = 0) число мостиков максимально (атах). Малое сокращение мышцы не сопровождается изменением а. В процессе растяжения параметр а практически линейно уменьшается с ростом и. Однако при значительных сокращениях (более 30% естественной длины) происходит перекрытие актиновых нитей, а миозиновые нити вступают в контакт со стенками саркомера [10]. В рамках настоящего исследования данное явление не учитывается и принята следующая зависимость для а(и) с условием в виде неравенств

Гатах - Си, и > 0; а (и) = \ тах ’ (4)

а™_ и < 0.

Здесь С - коэффициент пропорциональности.

56

Идентификация параметров модели

Система дифференциальных уравнений (3) с учетом зависимости (4) описывает деформирование мышцы в активном состоянии. Параметры модели к1, к2, к, с0, / 5, атах и С определяются при сопоставлении результатов расчета с экспериментальными данными. В работах [7, 9] выдвинута гипотеза о том, что генерируемая одним поперечным мостиком сила / и расстояние 5 между положениями равновесия тянущего и тормозящего мостика не зависят от типа мышцы и составляют соответственно

/ = 3 • 10-12 Н; 5 = 8 • 10~9 м. (5)

Принимая данную гипотезу, в дальнейшем будем считать параметры / и 5 известными.

Экспериментальное исследование сократительной активности мышцы проводят в изометрическом и изотоническом режимах. В изотоническом режиме к мышце в пассивном состоянии крепится груз заданной массы т. После электростимуляции мышцы регистрируется смещение у груза. Смещение у связано с удлинением мышцы и следующим образом: у = и0 - и . Здесь и0 - удлинение пассивной мышцы под действием силы тяжести mg. На рис. 3 представлены экспериментальные зависимости у(^) из работы [3] для портняжной мышцы лягушки при температуре 0 0С. На данных диаграммах можно выделить следующие характерные этапы:

1. Начальный этап (^ < 0,1 с), когда сокращение мышцы практически отсутствует и происходит рост силы в активной части мышцы. Длительность данного этапа возрастает с увеличением массы груза т.

2. Стационарное сокращение мышцы (0,1 < ^ < 0,33), которое характеризуется практически постоянной скоростью движения груза V. Этап начинается с того момента, когда сила в сократительном элементе уравновесит силу тяжести груза Р = mg. При стационарном сокращении сила, действующая на груз со стороны мышцы, не изменяется со временем: Р = Р = п = т = 0 . Тогда из системы (3) можно определить скорость сжатия мышцы V при постоянной нагрузке Р

V = -и| = 2 т^к 2 —2М/а тах - Р—. (6)

|Р*=п=0 1 2 Р(к1 + к2) + 2Ы/ак1

Данное выражение совпадает с уравнением Хилла (1), если параметры эмпирического соотношения Р0, а, Ь связаны с параметрами модели «скользящих нитей» к1, к2, атах, 5, / следующим образом:

Р0=2Ща, а = Р0 к , Ь = 22 . (7)

0 0 к1 + к2 к1 + к2

Если величины / и 5 известны, то из системы (7) можно при заданных параметрах уравнения Хилла Р0, а, Ь найти значения параметров к1, к2, атах:

—, к2 = -ЬР^, к1 =к2—^~ . (8)

2/ 2Ш5 Р0 - а

Параметры Р0, а, Ь определяются по экспериментальной зависимости ^Р). Как было показано в работе [2], для нахождения Р0, а, Ь достаточно измерить установившуюся скорость изотонического сокращения при трех различных нагрузках. В соответствии с зависимостями, представленными на рис. 3, для портняжной мышцы

сх

тах

0

0

0,1

0,2

0,3

0,4

0,5

0,6

/, с

Рис. 3. Временная зависимость изотонического укорочения портняжной мышцы лягушки. Цифры у кривых соответствуют нагрузке в граммах. Пунктирные кривые -экспериментальные зависимости из [3]; сплошные - результаты использования

математической модели (3)

лягушки (Ь0 = 7 см) Р0 = 0,42 Н; а = 0,27 Н; Ь = 3,2 см/с. Используя соотношения (8), получим атах = 2-10 ; к1 = 160 с- ; к2 = 88,9 с- . Для описания первого и второго этапов изотонического сокращения необходимо также знать значения параметров к, с0, С.

3. Этап быстрого уменьшения скорости сокращения до нуля (остановка груза). Уменьшение скорости сокращения при изотоническом режиме обусловлено прекращением работы поперечных мостиков. Прекращение актомиозинового взаимодействия происходит в тот момент, когда вся энергия, переданная мышце при электростимуляции, расходуется на увеличение внутренней энергии мышцы, ее деформирование и работу по перемещению груза. При описании данного этапа в феноменологической модели следует устранить сократительный элемент. Тогда поведение мышцы определяется упругим (с0) и инерционным (к) элементами, а третье уравнение системы (3) примет вид: кР + с0Р = с0ки. Используя динамическое

уравнение движения груза, прикрепленного к мышце (ти = mg - Р ), для удлинения и получим дифференциальное уравнение 4-го порядка

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

с0 кт к

(9)

ускорением у = -и = -g------. Имея экспериментальную зависимость у({) на третьем

к + т

участке сокращения, можно найти параметр к феноменологической модели. Для определения к необходимо измерить промежуток времени At, в течение которого скорость сокращения уменьшится от V до нуля. Тогда

к = mgA^L. (10)

V

Исходя из экспериментальной зависимости смещения груза от времени для портняжной мышцы лягушки при массе груза m = 9 г (рис. 3), получено к = 150 г. При известном значении к можно достаточно точно (рис. 3) описать третий этап изотонического сокращения.

4. Возвращение груза в исходное (до электростимуляции) состояние. Описание данного этапа с использованием того же дифференциального уравнения (9), что и для третьего этапа сокращения, возможно только при значительных нагрузках (массе груза m). При малой массе (менее 5 г) груз возвращается в исходное состояние значительно быстрее, чем это следует из расчетов по математической модели (3). Данное несоответствие расчетных оценок экспериментальным данным обусловлено, по-видимому, двумя обстоятельствами:

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

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

Параметры феноменологической модели с0 и С удобнее определять по результатам испытаний в изометрических условиях. Испытания в изометрическом режиме заключаются в получении и анализе временной зависимости силы P, генерируемой мышцей после электростимуляции при фиксированной длине [11]. На экспериментальной диаграмме P(t) (рис. 4) можно выделить характерные стадии: быстрый рост силы до максимального значения Pmax и медленную релаксацию после прекращения взаимодействия между нитями саркомера (определяемую в основном свойствами пассивной части мышцы и в настоящей работе не исследуемую).

Рассмотрим подробно стадию роста силы до максимального значения Pmax. Значение Pmax зависит от длины мышцы L, при которой осуществляется изометрическое стимулирование. Наибольшее значение Pmax соответствует естественной длине мышцы и совпадает с параметром P0 уравнения Хилла

PmaxIт-т = P) = 2N f amax . Если до электростимуляции в изометрических условиях

L—L0

мышцу растягивают на некоторую длину AL, то реализуемая сила Pmax оказывается меньше параметра P0. Сила Pmax практически линейно уменьшается с ростом смещения AL. Для описания этапа роста силы при изометрическом сокращении в системе уравнений (3) следует положить: u — AL — const, u — u - 0 . Тогда из системы (3) можно

определить значение Pmax, при котором прекращается увеличение силы (P — P — 0 ) и количества тянущих мостиков ( n — m — 0 )

Pnax — 2Nf a(AL) — 2N f (amx - C AL) — P0 - 2NfC AL . (11)

Таким образом, зная величину P0 и измеряя максимальную силу Pmax при различных удлинениях AL, определим параметры функции a(u) (4). На рис. 4 показана

зависимость Р(Т) для портняжной мышцы лягушки, предварительно растянутой на 3 см. Так как в этом случае Ртах = 0,12 Н, получим С = р>^ ~ 5' 10-5 м-1.

Коэффициент жесткости с0 выбирается таким образом, чтобы расчетная функция Р(Т) наиболее точно описывала экспериментальную зависимость. Предлагается следующая методика определения параметра с0.

а) Путем численного дифференцирования экспериментальной зависимости Р(Т)

определяется третья производная по времени от силы Р в начальный момент времени • - - Р

Рэкс(0) = 8—1. Здесь Р\ - первое ненулевое значение силы, регистрируемое при Т1

изометрических испытаниях; Т\ - момент времени, соответствующий силе Р1.

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

Р (0) = С0 *12N1 а(АЬ) = С0 к Рдах . к к

с) Выражение для начальной третьей производной от силы приравнивается значению Рэкс (0), полученному при обработке экспериментальной зависимости: Р с

8 “Г = ~Т к1 Ртах . Тогда Тл к

1

кР1

к Р Т п,11 тах '1

0=8. „ 1 3. (12)

Анализируя экспериментальную зависимость Р(Т) на рис. 4 и используя ранее найденные значения параметров феноменологической модели, получим, что для портняжной мышцы лягушки параметр с0 составляет « 84 Н/м.

Описание механического поведения мышцы

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

П = к, (а(м0) - п - т)-Р—п, т = -к2т л---------------Р—п,

и 0 ' 2Шс0 2 2N5c0

(13)

Р = 1 [2Щ(п - т) - Р]. к

Начальные условия для уравнений (13) имеют вид: п(0) = т(0) = Р(0) = Р(0) = 0 . Представленная на рис. 4 расчетная зависимость Р(Т) при изометрической активации является решением системы (13) при и0 = АЬ = 3 см. Можно отметить приемлемую точность описания экспериментальной зависимости генерируемой мышцей силы от времени (рис. 4). Наличие «всплеска» на расчетной диаграмме Р(Т) (Т « 0,12 с)

г, с

Рис. 4. Временная зависимость силы, генерируемой мышцей при изометрической активации. Пунктирная кривая — экспериментальная зависимость из [10]; сплошная -результат использования математической модели (3)

обусловлено пренебрежением в рамках модели вязкостью нитей саркомера. В изометрических условиях сила Р, определяемая системой (13), возрастает до Ртах (11).

Первый этап изотонического сокращения заканчивается, когда сила Р, также определяемая из системы (13), уравновесит силу тяжести груза Р = mg. При этом рост силы прекращается и начинается движение груза. Система (3) для данного процесса имеет вид

п = к, (а(и) - п - т)н——п, т = -к2т — п,

и ;2Ж 2 N5 (14)

и =1 [2Щ (п - т) - mg]. к

Начальные условия для системы (14) следующие:

и (0) = и0, и (0) = 0, п(0) = п1, т(0) = т1. Здесь п1 (т1) - количество тянущих

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

Заключение

В результате использования разработанной феноменологической модели мышцы и соотношений теории «скользящих нитей» установлено, что введение в модель инерционного и упругого элементов (учитывающих инертность и деформативность нитей саркомера, соответственно) позволяет интерпретировать и

адекватно описать сократительную активность мышцы в изотонических и изометрических условиях.

Список литературы

1. Шилько, С.В. Механика адаптивных композитов и биоматериалов / С.В. Шилько, Ю.М. Плескачевский // Материалы, технологии, инструменты.- 2003.- Т. 8, № 4.- C. 5-16.

2. Черноус, Д.А. Феноменологическая модель мышцы как активного композита / Д.А. Черноус, С.В. Шилько // Механика композитных материалов.- 2006.- Т. 42.- № 5.

3. Hill, A.V. First and last experiments in muscle mechanics / A.V. Hill. - Cambridge, 1970.

4. Бегун, П.И. Биомеханика: учебник для вузов / П.И. Бегун, Ю.А. Шукейло. - СПб.: Политехника, 2000.

5. Черноус, Д.А. Модель резонансной дисперсии полимеров / Д.А. Черноус // Материалы, технологии, инструменты. - 2004. - Т. 9, № 3. - С. 22-26.

6. Wong, A.Y.K. Mechanics of cardiac muscle, based on Huxlev’s model: mathematical simulation of isometric contraction / A.Y.K. Wong // Journal of Biomechanics. - 1971. - Vol. 4, No. 6. - P. 529-540.

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

8. Wong, A.Y.K. Some proposals in cardiac muscle mechanics and energetics / A.Y.K. Wong // Bull. Math. Biol. - 1973. - Vol. 35, No. 3. - Р. 375-399.

9. Шестаков, Д.А. Математическая модель механических свойств волокон скелетной мышцы с учетом растяжимости актиновых нитей / Д.А. Шестаков, А.К. Цатурян // Биофизика. - 1998. - Т. 43, № 2. -С. 329-334.

10. Лавендел, Э.Э. Моделирование процесса сокращения сердечной мышцы / Э.Э. Лавендел, Л.И. Мачабели, И.О. Типанс // Механика композитных материалов. - 1981. - № 6. - С. 1088-1092.

11. Дубровский, В.И. Биомеханика: учебник для средних и высших учебных заведений /

В.И. Дубровский, В.Н. Федорова. - М.: ВЛАДОС-ПРЕСС, 2003.

MODELLING OF CONTRACTIVE ACTIVITY OF THE MUSCLE TISSUE

D.A. Chernous, S.V. Shilko (Gomel, Belarus)

A modified phenomenological model of the muscle that accounts for the inertial and deformation properties of sarcomere fibers and avoids certain inconsistencies and inaccuracies found in the previous descriptions is proposed. A computation method for deformation parameters of the muscle tissues in the active state (phase of contraction) based on the model developed and known relations of the theory of sliding fibers is set forth. The model is compared with available experimental data on deformation of the frog tailor’s muscle in the isotonic and isometric regimes.

Key words: muscle tissue, contractive actvity, sarcomere fibers, isotonic and isometric loading, phenomenological model.

Получено 15 августа 2006