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

МОДЕЛИРОВАНИЕ ДИНАМИКИ НЕЛИНЕЙНОГО СОКРАЩЕНИЯ САРКОМЕРА

С.А. Охотников (Уральский государственный университет им. А.М. Горького) Научный руководитель - к.ф.-м.н., профессор Г.П. Быстрай (Уральский государственный университет им. А.М. Горького)

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

Динамика нелинейного сокращения саркомера

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

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

Подставляя (1) в динамическое уравнение продольного сокращения относительно оси растяжения 1 ёг ф ш

получаем термодинамическое нелинейное уравнение, описывающее сокращения сар-комера, в виде

шг = -Ф(1оХг3 -10кг2 +1ог-сте), (2)

ш

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

ds dG ^ ц „ys ц Qks 1 2 ...

dS =-paG • G =bir-V+2 2-а's- (3)

Стационарная скорость укорочения. Вводя для напряженного саркомера скорость укорочения и внутренние напряжения в виде ёг

и = —, ст = ст. = 10 г, Л г

представим, после деления правой и левой частей на итах, уравнение (2) в виде

= - A

Г „Л3 Г „Л2 г „л

и m

а

а „

+ B

а

- C

V°Q J

а

V°Q J

max \ Q J

где соответствующие константы равны

+ D , (4)

A =

ФХа0

2

M 0 umax

B =

ф£а0

; С =-фа^- D = _Фа*.

ит

ит

M о Umax "max "max

При а / ао = О и / umax = 1, тогда D=1. В результате получаем теоретическую зависимость скорости укорочения от а / а о (рис. 1, кривая 2: A=7.3; B=8.15; С=3.95), которая соответствует экспериментальным результатам [3].

и

и max 0,8 - ч

0,6

0,4 - • х

0,2

0 - • \ •

0,25 0,5 0,75 а / а

0

Рис. 1. Зависимость относительной скорости и / итах от относительной силы (напряжения) а / ад. Кривые 1 и 2 отличаются различными значениями констант и показывают

рабочие диапазоны действия нелинейной модели. Точки соответствуют экспериментальным результатам [3]

Редукция к классическим моделям

Модель А.Хилла. Полагая в уравнении (4) A=0 и B=0, получаем и

= -С

и m

+ D.

Чао J

Делая замену в последнем выражении следующим образом

а P С — = —, С = а-

b ( P Л b P

b 1 + ^1, d = —P

а

aa

и m

a J U max a

где a - малый параметр, получаем уравнение А. Хилла [4] и b

и m

b P0 + a

и и P + a

max max

- + -

Коэффициент D был выбран из следующих соображений: а/а0 = 0 и / umax = 1,

тогда D=1. Это условие удовлетворяет выражению для максимальной скорости стационарного укорочения

и = bP0.

max

a

Модель Б.С. Эббота и Д.Р. Уилки. Принимая в выражении (4) коэффициент A=0, получаем

и = вТлГ - с(-а1+d .

Ча0 J

Чао J

Для перехода к модели Б.С. Эббота и Д.Р. Уилки необходимо сделать замену

0

a F 2

— = —, B = a

C = a-

b f F,

a 0

aa

и m

и

max V

1 + _L I, D = ■ a J

b F]

U max a

С помощью этой замены получаем уравнение Б. С. Эббота и Д.Р. Уилки [5]

и

b Fl — F

U max U max a + F

Таким образом, данная модель может считаться обобщением модели А.Хилла, так как содержит еще и квадратичный член по напряжению.

Соответствие II началу термодинамики. Уравнение (2) должно быть совместимо с условием положительности производства энтропии

dtS dt

1 1

1

= J(в)8 = Цоs2I---ks + —%s2 !> о

2 3

4'

где константа до = д0/Т0 имеет размерность Дж/К, а параметр Т0 - температура, при которой происходит сокращение саркомера. Выражение для производства энтропии /dt может быть представлено выражением:

G*■ = 1 s*2(*2 -4s0s* + б)> О,

(5)

где в* =в/8с, 8с = 1/3%, в0 = 80/ес = к/3%8с. Далее необходимо рассматривать решения неравенства вида в*2 - 4в0в* + 6 > 0, входящего в (5), в зависимости от дискриминанта Б = Ь2 - 4ас ; а = 1, Ь = -4в 0, с = 6. Для его выполнения необходимо, чтобы па-

раметры, входящие в него, удовлетворяли следующим неравенствам:

* 2 3 , 2 9 s 0 < —, или k < — % .

2

2

(6)

Делая замену переменной в (5) в* = п + в*, получаем

О * =1 п4 +1 а У + а, п + О0 * > 0, 4 2 0

* ^ / *2 1 N ^ * /-4*3 /Т * ^ *2 \Г\ *2 I / Л Т-\

здесь а =-3(в0 -1), аж = 3в0 -2в0 , О0 = 3в0 (2-в0 у4, а п- параметр порядка. Выражение (6) не содержит внешнее поле а е = а = 0. Положительности производства энтропии отвечает условие а * > -3/2 (см. рис. 2). Отсюда следует, что функции О*1 > 0,

О*1 < 0 являются функциями Ляпунова. Приведем аналог теоремы Пригожина для нелинейных систем, которая впервые сформулирована в [6].

С*'

Рис. 2. Производство энтропии в канонической форме. 1 - s0 =1.08, 2 - s0=1.155,

*

3 - s0=1.204. Нулевые значения производства энтропии соответствуют s = 0 в выра-

жении n = s* — s0. Второй экстремум Gсоответствует стационарному состоянию

b

Теорема. Временная эволюция в нелинейной термодинамической системе при заданных постоянных граничных условиях происходит так, что производство энтропии О ' стремится убывать йО *

йг

< 0

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

Приведенная знакопеременная потенциальная функция, равная относительной (безразмерной) скорости изменения энтропии системы также примет вид:

йЗ *—( * 7 1 4 1 * 2 т * >0 Л Л Л |>0

О = —г = °е + = - П4 + - «л2 + ь п>°, О *е =-(О* + а*п)| <°.

Здесь обратимые потоки энтропии О*е также могут принимать разные знаки, поэтому

потенциальная функция О * может иметь любой знак.

Хаотическая динамика параметра порядка. Используя переход от релаксационных уравнений к уравнениям второго порядка и учитывая эффект последействия [6], получаем однородное каноническое уравнение второго порядка для величины деформации п в первом приближении:

тг п+Г(г) п+п3 + а*п = °0*С08 ®г,

(7)

здесь коэффициент затухания и амплитуда внешней силы равны соответственно Г(г) = 1 -т(3п2 + а') > 0, а 0* = а*0(1 + (иг)).

100 г

Рис. 3. Поведение приведенного функционала производства энтропии Ок*7 (а, б) в условиях воздействия периодической внешней силы Ь* =1.9, ш = 2.35, тг = 1.1, т = 0.216, е0* = 1.204, /о=0.34е. Начальные условия для трех переменных: 0.1,

0.02, 0.01. Линией указано среднее значение функционала во времени Ок *7 = 0.381;

продолжительность всей истории движения Л=/(п)=100, шаг разбиения Д/=0.01; в) вид функционала от параметра порядка, определенного в предшествующий момент

времени с задержкой Д= 20 расчетных точек

Рассмотрим поведение производства энтропии при наличии хаоса (рис. 3). Расчет производился по функционалу

-—у 1 4 1*2 * *

О* = -п- + 2а п п* + Оо,

в который подставлялись значения решений уравнения (7) (к - индекс шага расчета); здесь

3

* о/ *2 14 * о* о*3^->* *2/О *2 ч

а = -3(в 0 - ^ = 3в 0 - 2В0 , О0 = - В0 (2-е 0 ).

На рис. 3, в представлено возмущенное состояние саркомера с временной задержкой Д= 20. В отличие от невозмущенного состояния (рис. 3, а) данный график характеризуется замкнутыми областями, отвечающими за определенные невозмущенные стационарные состояния. Центры этих областей являются глобальными или локальными минимумами.

Согласно второй теореме об устойчивости функционала [7], производство энтропии и его производную можно оценить некоторыми числами сверху. Таким образом, перебирая различные значения параметра е* можно получить различные средние значения функционала Ок*', которые будут ограничены числами снизу (0 - равновесное состояние) и сверху Ок .

Пульсации температуры

В данной работе считалось, что существенное изменение средней температуры происходит на расстояниях в0 (основной масштаб сокращения), на которых меняется

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

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

Пульсации температуры по Обухову. В теории турбулентности вводят энергетический, инерционный и диссипативный интервалы пульсаций [8, 9]. Инерционный интервал берет энергию из больших вихрей и снабжает энергией малые пульсации. Инерционный интервал является в то же время конвективным - выравнивание температур в нем происходит путем механического перемешивания различно нагретых «частиц» без участия истинной теплопроводности; свойства температурных пульсаций в этом интервале не зависят и от крупномасштабного движения. В случае саркомера, находящегося в растворе, считается, что число Прандтля Рг= х / V =1 и поле температур подобно полю скоростей деформации. Последнее означает, что коэффициенты температуропроводности и кинематической вязкости для инерционного интервала равны

X = v = sy•nУ, /1; в^ и пУ- приведенные пульсации пространственного масшта-

ба и скорости пульсаций, отнесенной к средней скорости сокращения. Определим относительное повышение температуры ДТ * = ДТ / Т0 в инерционном интервале по Обухову [8], где Тз =309.6 К - температура внутренней среды организма. Тогда скорость

диссипации энергии за счет теплопроводности со значением х = V = ву • пУ дается вы-

ражением [9]

ф* =х

'лт*>2

у Лх у

= (8УПУ )

дт

*2

.У 2

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

дт =

Ф*

п

(8)

здесь ф - скорость диссипации энергии за счет теплопроводности: ф* = Ет / ерТ0 = Ет 1, Ет - приведенная средняя энергия диссипации, ср - удельная теплоемкость при постоянном давлении. Представим Ет у в виде полинома четвертой степени от пространственного масштаба :

Ету= 4Гву )4 + 2а(ву )2 + ЬвУ ,

где а и Ь - некоторые приведенные параметры. В формуле (8) неизвестно значение ф* = Ет''. Для его нахождения воспользуемся выражением

1 "

=NI Ет' С п ),

1 у п=1

где N - общее число точек, на которое делится функция Ет по оси еу, 8У п -значение 8У в точке п. Для нахождения приведенной скорости пульсаций пУ используем

=

Лв^

Ж

. Формула (8) дает достаточно реальные значения повышения температуры 4-

16оС для саркомера в растворе. На рис. 4 приведена зависимость приведенных пульсаций температура от приведенного значения пространственного масштаба в7 .

Расчет средней диссипации энергии. Расчет производился для раствора хлорида кальция СаС12 [10] (ср = 40.88 Дж/(моль-К) [11]) по формуле Ет = Ет 1 срТ0. В [3] для

портняжной мышцы лягушки при сокращении на 2 мкм энергия рассеивания равна 0.13 Дж/г. Данное значение соответствует расчетной пульсации температуры при сокращении саркомера в 4.7 К.

дт

0 05

О 0196

а

Е 7

-0.5

0.5

б

Рис. 4. Зависимость приведенной температуры Дт (а) и приведенной диссипации энергии Ет7 (б) от пространственного масштаба пульсаций а7 при а= -0.07, Ь=0,

Ет7 =3.309-10-4 , Дт = 4.7 К

у

г

г

У

8

У

8

Спектр мощности пульсаций параметра порядка п- Для модели построены нормированные спектры мощности пульсаций параметра порядка п методом Фурье-преобразования (рис. 5). Из рис. 5 видно, что теоретически полученный спектр является сплошным, что свидетельствует о существовании многомасштабной структуры поля скорости деформации саркомера [12]. Именно многомасштабность и является важнейшим признаком развитой турбулентности, приводя к возбуждению гигантского числа степеней свободы [13]. Полученный спектр соответствует экспериментальным данным, описанным в книге [14]. В области больших частот (ш>4-104 сек-1) спектральная плот-

—2 — 1 ность изменяется по закону ~ ш . С уменьшением частоты (104 - 4-104 сек ) последняя резко возрастает, что соответствует закону ~ ш-7 турбулентных пульсаций Гейзен-берга в диссипативном (вязком) интервале [15]. В достаточно узкой области частот

(103-104 сек-1) спектр можно сравнить с колмогоровским Sn~ ш-5/3, который всегда наблюдался в развитой турбулентности [16, 17]. В области еще меньших частот

(ш<1-103 сек-1) в узком интервале имеет место фликкер - шум Sn~ ш-1 [13]. В области

инерционного интервала происходит переход от спектра Sn~ ш-7 к колмогоровскому

-5 /3

спектру Sn~ ш .

Sn

O.Ol -

1 10

1 10 -

1 10 -

1 10 " -

1 10 -

1 10 -

1 10

1 -10

œ = 2. б

10

100

1 10

1 10

1 10

1 -10

œ, с

Рис. 5. Зависимость спектра пульсаций параметра порядка п (Зп) от частоты®

при следующих значениях параметров а * = -1.5 , и = 2.6 (частота внешних гармонических колебаний), т = 0.216 , п(0)=0.1, о0* = 1.8. На спектре проявляется пик при ю = 2.6

Цикл реакций, проходящих при сокращении саркомера

Миозин обладает ферментативной активностью. Он катализирует гидролиз АТФ. Активные центры у миозина расположены в его «головках» [18]. Рассматривается следующая модель: саркомер находится в растворе, к которому добавляют АТФ. В этой системе идут химические реакции, однако внешних периодических воздействий нет. Цикл реакций с участием АТФ и саркомера представляет следующую схему: X 0 + 2 X! о 2 X 2, (1)

2X2 о 2X3 + X4, (2)

2X3 о X5, (3)

X5O X6, (4)

X6 о X7 + 2X9, (5)

X6 о Xо + 2X8 + 2X9, (6)

X7 о Xо + 2X8, (7)

где X0 = AM - актинмиозиновый комплекс, X1 = ATP, X2 = M.ATP - комплекс миозина и молекулы АТФ, X3 = M.ADP.Pi - комплекс миозина, молекулы АТФ и фосфора, X4 = H + - ион водорода, X5 = AM.ADP.Pi, X6 = AM * .ADP.Pi, X7 = AM ". ADP, X8 = ADP, X9 = Pi - фосфор. Данная схема построена на основании старых моделей

[19-21] c добавлением реакции (6), которая может идти параллельно с (5). В данном цикле учтено также участие 2 «головок» миозина в реакциях (1)-(3) и (5)-(7), поэтому в кинетических уравнениях возникнут квадратичные члены. Реакции (1), (2) и (6), описанные в работах [19, 20], протекают не мгновенно. Процесс (3) усовершенствован для 2 «головок» миозина. Представляемая модель содержит 7 реакций, в которых участвуют 10 веществ. Этими реакциями описываются следующие процессы:

(1) - присоединение АТФ к «головкам» миозина с образованием АТФ-миозинового комплекса;

(2) - гидролиз M .ATP с образованием M. ADP.Pi комплекса и ионов H + . Ионы водорода в дальнейшем уходят в водную среду;

(3) - образование из комплекса M.ADP.Pt вещества AM.ADP.Pt. В ходе данной реакции происходит продвижению головок миозина к актину;

(4) - образование энергетически активной конформации миозина AM * .ADP.Pi ;

(5) - изменение конформации легкой части миозина, которое происходит при распаде AM * .ADP.Pi на AM * .ADP и Pi («головка» миозина производит тянущее усилие);

(6) - описывает распад AM *. ADP.P на соответствующие продукты (наряду с (5)). При

распаде происходит мгновенное выделение энергии, а на механическом уровне «головки» миозина производят тянущее усилия;

(7) - распад AM * .ADP с выделением ADP (выделяется энергия).

Эти реакции служат основой для написания кинетических уравнений. Кинетические уравнения. Используя закон действующих масс и принцип макроскопической обратимости [22-24], можно записать кинетические уравнения для реакций (1)-(7), приведенные к безразмерному виду и определяющие изменения концентраций веществ xi (i = 0,.. .,9) со временем.

dx , 2 2 ,22 2,2

— = k1(x2) - k1 x0(x1) + k 7 x6 - k 7 x0(x8) (x9) + k6 x7(x8) - k6 x0 (x8 ) , dt

dx1 2 2 —— = k1(x2) - k1 x0(x1) , dt

dx2 2 2 2 2

— = -k1(x2) + k1 x0(x1) - k2(x2) + k2x4(x3) , dt

dx3 2 2 2

— = k2 (x2) - k2 x4(x3) + k3 x5 - k3 x4(x3 ) , dt

dx

dt

4 = k2 (x2) k2 x4(x3) + k3 x5 k3 x4(x3 ) '

dx

dt dx,

— ^4 x5 I k4 x^ k3 x5 I k"^ x4 ( x^ )

dt dx

k4 x5 k4 x6 k5 x6 1 k5 x7 (x9 ) k7 x6 1 k7 x0 (x8 ) (x9 )

dt

dt dx

— k5 x6 k5 x7(X9) k6 x7(x8 ) + k6 x0(x8) J

dt

— k6 x7 (x8 ) k6 x0 (x8 ) ,

— k 5 x6 k 5 x7(x9) -k6 x7(x8) + k6 x0(x8) + k7 x6 - k7 x0 (x8 ) (x9 ) J

где кi и к' - константы (/ = 0,...,9), г=г/г0,а х7. приведены в безразмерном виде. Таким

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

10 О

1 1 1

2

а 1 3 1

x,

20

40

t

1

1 1 б

1 1

0 20 40 t

Рис. 6. Зависимость концентраций веществ X0 — AM (а), X, — ATP (б)

от времени /ъ=1.6 мс

xe

x~i

Рис. 7. Зависимость концентраций веществ Х6 = АМ*.АБР.р (а), Х7 = АМ*.АБР (б), Х8 = АБР (в), Х9 = Р{ (г) от времени £///0, /0=1.6 мс. 1 - временной ход реакции без самовозбуждения, 2 - самовозбуждение, 3 - выход на стационарное значение

Численные решения. Решая полученную систему численными методами, были полу-

3

t

чены следующие результаты. Принимая во внимание, что весь процесс составляет 80 мс [3, 4], получаем время t0=1.6 мс. Конформация легкой части миозина происходит при распаде AM * .ADP.Pi на AM * .ADP и Pi или на составляющие это вещество части. Полный цикл завершается после выделения ADP. На рис. 6, а, можно выделить 3 участка, отвечающих механике сокращения саркомера: 1 - процесс сокращения саркомера; 2 - релаксация саркомера к стационарному состоянию; 3 - процесс так называемого "дрожания" саркомера. Не весь АТФ расходуется на конформацию головок миозина и на полимеризацию актина (см. рис. 6, б). Концентрация ионов водорода возрастает до определенного момента времени. Часть ионов связывается с АТФ [3, 18], а часть - с ионами магния, что ведет к релаксации саркомера к невозбужденному состоянию (рис. 6, а, участок 2). В модели В.И. Дещеревского [25] отсутствует участок 3, представленный на рис. 6, а участки 1 и 2 соответствуют тянущим и тормозящим усилиям при сокращении саркомера. На рис. 7 представлены графики самовозбуждения веществ X6 = AM * .ADP.P , X7 = AM * .ADP, X8 = ADP, X9 = Pi. История процессов на рис. 7

состоит из 3 стадий: 1 стадия - предварительное нелинейное изменение соответствующих концентраций (73.6 мс); 2 стадия - начало самовозбуждения, содержащие высокочастотные пульсации и периодические движения (4.8 мс); 3 стадия - развитие неустойчивых низкочастотных пульсаций, которые приводят к неустойчивости процесса (1.6 мс). Как видно из рисунков, самовозбуждение в системе происходит после выделения ADP . В свою очередь, это ведет к началу процесса релаксации саркомера в растворе.

Литература

1. Быстрай Г.П., Богинич А.В. Термодинамика многоядерных клеток: системное моделирование самоорганизующегося саркомера с хаотической динамикой параметра порядка // Вестник кибернетики. - 2007. - № 6. - ИПС СО РАН - С. 77-91.

2. Быстрай Г.П., Макаров Л.В., Шилин Г.Ф. Неравновесная термодинамика процессов горного производства. - М. Недра, 1991. - 119 с.

3. Рубин А.Б. Биофизика. Т. 2. - М.: Наука, 2000. - 467 с.

4. Хилл А. Механика мышечного сокращения. - М.: Мир, 1972. - 182 с.

5. Wilkie D.R. The mechanical properties of muscle // British medical bulletin. - 1956. -V.12.

6. Быстрай Г.П. Термодинамика открытых систем. Учебное пособие. - Екатеринбург: Изд-во Урал. университета (гриф УМО). 2007. - 120 с.

7. Ким А.В. Ко второму методу Ляпунова для систем с последействием // Диф. уравнения. - 1985. - Т.21. - №3. - С. 385-391.

8. Обухов А.М. Структура поля температуры в турбулентном потоке // Изв. АН СССР. - Геогр. и геофизика. - 1949. - Т.13. - С. 58-69.

9. Ландау Л.Д., Лифшиц Е.М. Гидродинамика. - М.: Наука, 1988. - 733 с.

10. Волькенштейн М.В. Физика мышечного сокращения // УФН. - 1970. - Т. 100. -Вып.4. - С. 681-717.

11. Кикоин И.К. Таблица физических величин. - М.: Атомиздат, 1976. - 1008 с.

12. Шустер Г. Детерминированный хаос. Введение. - М.: Мир, 1988. - 240 с.

13. Мирабель А.П. О роли нелокальных взаимодействий в формировании спектра пассивной примеси в двумерной турбулентности // Этюды о турбулентности. - М.: Наука, 1994.

14. Колесниченко А.В., Маров М.Я. Турбулентность многокомпонентных сред. - М.: МАИК Наука, 1998. - 336с.

15. Хинце И.О. Турбулентность. - М.: Физматгиз, 1963.

16. Meksyn D. New methods in laminar boundary layer theory. - Pergamon Press, London, 1961.

17. Белоцерковский О.М., Андрющенко В.А. Динамика пространственных вихревых течений в неоднородной атмосфере. - М.: Янус-К, 2000.

18. Румянцев Е.В., Антина Е.В. Химические основы жизни. - М.: Химия, КолоС. 2007.-560с.

19. Siththanandan V.B., Donnelly J.L., Ferenczi M.A. Effect of strain on actomyosin kinetics in isometric muscle fibers // Biophysical Journal. 2006. - V.90. - Р. 3653 - 3665.

20. Ranatunga K.W., Coupland M.E., Pinniger G.J., Roots H., Offer G.W. Force generation examined by laser temperature - jumps in shortening and lengthening mammalian (rabbit psoas) muscle fibres // J. Physiol. - 2007. - V.585. - №1. - Р. 263-277.

21. Bendall J. Muscles, molecules and movement. Heinemann, Lnd., 1969.

22. Murray J.D. Mathematical biology. Springer - Verlag Berlin Heidelberg. 1989. - p. 767.

23. Николис Г. Пригожин И. Самоорганизация в неравновесных системах. - М., 1979. -512 с.

24. Жуховский А.А., Белащенко Д.К., Бокштейн Б.С., Григорян В.А., Григорьев Г.А., Гугля В.Г. Физико-химические основы металлургических процессов. - М.: Металлургия, 1973. - 392с.

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