CC BY
47
Научный журнал КубГАУ, №99(05), 2014 года
1
УДК 51-76
УСТОЙЧИВОСТЬ СТАЦИОНАРНЫХ СОСТОЯНИЙ КИНЕТИКИ ЛЕЙКОПОЭЗА
Тумаев Евгений Николаевич д.ф.-м.н., профессор
Кубанский государственный университет, Россия, 350040, Краснодар, Ставропольская, 149, [email protected]
Шарай Иван Александрович аспирант
Кубанский государственный университет, Россия, 350040, Краснодар, Ставропольская, 149, [email protected]
В статье приведены результаты исследования устойчивости модели нейтрофиломоноцитопоэза. При помощи критерия Рауса-Гурвица вычислено, что приведенная система дифференциальных уравнений, описывающих созревание клеток, является асимптотически устойчивой.
Определены пороговые значения параметров модели, при которых система становится неустойчивой
Ключевые слова: НЕЙТРОФИЛОПОЭЗ, МОНОЦИТОПОЭЗ, КОСТНЫЙ МОЗГ, УСТОЙЧИВОСТЬ
UDC 51-76
STABILITY OF STATIONARY CONDITIONS OF KINETICS OF LEYKOPOEZ
Tumayev Evgeny Nikolaevich
Dr.Sci.Phys.-Math., professor
Kuban State University, Krasnodar, Stavropolskaja
str., 149, Russia
Sharay Ivan Aleksandrovich postgraduate student
Kuban State University, Krasnodar, Stavropolskaja
str., 149, Russia
The results of the research of stability of the model of neutrophilomonocytegenesis are shown in the article. With the criterion of Routh-Hurwitz it's calculated that the system of the differential equations of cells growing is asymptotically steady. Threshold values of parameters of model at which the system becomes unstable are defined
Keywords: NEUTROPHILOGENESIS, MONOCYTEGENESIS, MARROW, STABILITY
Введение
Отличительной особенностью лейкопоэза служит то, что он представляет собой главную защитную систему организма. Один из важнейших показателей для оценки работы лейкопоэза - кинетика кроветворения и кроверазрушения.
Для описания данных процессов могут использоваться дифференциальные уравнения, описывающие последовательные переходы клеток из одной фазы созревания в другую. Так, все клетки-нейтрофилы проходят следующие стадии: единая колониеобразующая единица
гранулоцитарно-макрофагального рядов (КОЕ-ГМ), колониеобразующая единица гранулоцитов (КОЕ-Г), нейтрофильный миелобласт (НМб), нейтрофильный промиелоцит (НПм), нейтрофильный миелоцит (НМ), нейтрофильный метамиелоцит (НМм), палочкоядерный нейтрофил (Пн), сегментоядерный нейтрофил в костном мозге (Снкм), сегментоядерный
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
2
нейтрофил в крови (Снк), сегментоядерный нейтрофил в тканях (Снт). В свою очередь для моноцитов выделяют стадии: КОЕ-ГМ,
колониеобразующая единица моноцитов (КОЕ-М), монобласт (Мб), промоноцит (Пм), моноцит в крови (Мн) и макрофаг (Мф) [1-3].
Полученная таким образом математическая модель будет зависеть от большого числа параметров, влияющих на кроветворение. Модели нейтрофилопоэза и монопоэза в отдельности были опубликованы ранее в работах [4-5]. Однако в этих работах не рассматривался вопрос об устойчивости системы дифференциальных уравнений.
В таком случае представляется важной проверка устойчивости модели, так как даже малые колебания нормальных условий могут привести к серьезным изменениям в производстве нейтрофилов и моноцитов. Помимо этого, вопрос о важности исследования устойчивости стационарных состояний для биологических систем неоднократно обсуждался другими авторами, в частности [6-7].
В данной работе проверяется устойчивость одной из ветвей лейкопоэза - системы производства нейтрофилов и моноцитов, имеющих единого предка КОЕ-ГМ.
Цель исследования
Исследование устойчивости стационарного состояния нейтрофиломоноцитопоэза.
Методы исследования
Устойчивость стационарного состояния модели проверялась при помощи критерия Рауса-Гурвица.
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
3
Расчеты кинетики переменных и решение дифференциальных уравнений выполнены при помощи математического комплекса MATHCAD 14.
Результаты исследования
В уравнениях (1-15) представлена исследуемая модель производства нейтрофилов и моноцитов, основанная на уравнениях, опубликованных в более ранних работах [4-5]:
dn0
dt
10 + (n0 - G00 )
n0
G00 ^
r0 -(1 -70)(n0 -G0o)ko,
1
(1)
+ (n н1
dnH1
dt
(
G 0 н1)
1
= сн (1 -У0 )(n0 -
Пн1 - G0н1
K Гн1
K н1 )
G00 )kн0 + (1 -Ун1 )(пн1
G0 н1 )k н1,
(2)
dn
н2
dt
= (1 - a)(1 - Гн1 )(n н1 - G0 н1 )kн1 +
+ (n н2 - G0 н2 )
L nн2 - G0
1
н2
K
н2
гн2 - (1 - Ун2 )(nн2 - G0н2 )kн2
(3)
dnн
dt
= (1 - Ун{і-1) )nн(і-1) - G0н(і-1) kн(i-1)
+
+ (пні - G0ні )
( nн - G0 н ^
K
(4-5)
гні - (1 - Уні )(пні - G0ні )kні, i = 3,4
ні )
dn3l
dt
(1 Ун4 )(пн4 G0н4 )kн4 nн5^н5,
(6)
dnHj
dt
- nн(і—1)kн(і-1)
n rnkні,
і - 6 - 9.
(7-10)
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
4
dnMl
dt
= см(l ~Yo )(п0
G00)k мі +
+ (пMl - G0Ml )
L Пмі - G0Ml
1 -
мі
K
мі
V -(l -Умі )(п мі - G0 мі )kMl,
(11)
dn
м2
dt
= (l - a)(l - Умі )(пмі - G0мі )kмі +
+ (пм2 - G0м2 )
пм2 - G0м2
K
м2
(12)
гм2 - (1 - Ум2 )(пм2 - G0м2 )kм2
1
dn
м3
dt
= (l -Ум2 )(пм2 - G0 м2 )kм2 +
+ (nм3 - G0м3 )
L пм3 - G0
1 -
м3
K
м3
Гм3 (1 Ум3)(пм3 G0 м3)k м3,
(13)
d^4
dt
(l - ум3 )(пм3 - G0 м3 )kм3 - пм4kм4,
(14)
dnм5
dt
— пм4kм4
пм5^м5 •
(15)
где I0 - поток клеток КОЭ-ГЭММ.
G00 - число клеток КОЭ-ГМ в фазе обратимого покоя на кг.
G0Hi - число клеток нейтрофилов в фазе обратимого покоя на кг. G0Ш - число клеток моноцитов в фазе обратимого покоя на кг. а - потеря клеток вследствие неэффективного гемопоэза. у0 - коэффициент потери клеток КОЭ-ГМ путем апоптоза. уні - коэффициент потери клеток нейтрофилов путем апоптоза. уш- - коэффициент потери клеток моноцитов путем апоптоза. сн - доля клеток КОЭ-ГМ дифференцирующихся в нейтрофилы. см - доля клеток КОЭ-ГМ дифференцирующихся в моноциты. k0 - скорость потока КОЭ-ГМ в следующую стадию. kHi - скорость потока нейтрофилов в следующую стадию.
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
5
kMi - скорость потока моноцитов в следующую стадию.
п0 - общее число клеток КОЭ-ГМ на кг массы тела.
пні - общее число клеток нейтрофилов на стадии созревания на кг.
пмі - общее число клеток моноцитов на стадии созревания на кг.
r0 - скорость роста КОЭ-ГМ.
гні - скорость роста нейтрофилов на стадии.
гмі - скорость роста моноцитов на стадии.
К0 - поддерживающая ёмкость среды для КОЭ-ГМ.
Кні - поддерживающая ёмкость среды для нейтрофилов на стадии.
Кмі - поддерживающая ёмкость среды для моноцитов на стадии.
В таблицах 1 -3 приведены численные значения параметров модели в норме, рассчитанные по аналогии с [4-5].
Таблица 1 - Данные по состоянию КОЕ-ГМ в норме.
П0 3,83 • 105 Г0, ч-1 0,19
О О 3,26 • 105 К0 7,46 • 104
k0, ч-1 0,046 І0, ч-1 39,99
сн 0,88 см 0,12
Y0 0,03
Таблица 2 - Данные по состоянию нейтрофилов в норме.
Пні и О km , ч'1 Гні , ч-1 Кні Уні
КОЭ-Г 6,4 • 106 5,4 • 106 0,15 0,64 О Os 0,03
НМб 1,2 • 107 9,8 • 106 0,50 1,03 Js> “о о о.
НПм 00 О 00 О 0,05 0,07 00 О іп
НМ 2,0 • 109 О 0,19 0,52 00 О
НМм 2,9 • 109 - 0,02 - - -
Пн Os О іп - 0,01 - - -
Снкм 3,8 • 109 - 0,01 - - -
Снк 00 О - 0,12 - - -
Снт Os О 1/"Г - 0,01 - - -
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
6
Таблица 3 - Данные по состоянию моноцитов в норме.
пмі О S кмі, ч-1 Г мі, ч'1 Кмі Умі
КОЭ-М 1,8 • 106 о о\ 0,08 0,03 3,4 • 105 0,03
Мб 3,8 • 107 3,2 • 107 0,08 0,29 “о о о.
Пм 6,0 • 108 ОО О 1/"Г 0,02 0,03 ОО О оу
Мн С~~- О CD ОО - 0,02 - - -
Мф О - 0,001 - - -
Стационарное состояние лейкопоэза характеризуется в среднем постоянной численностью на разных стадиях созревания. Математически стационарному состоянию отвечает равенство нулю всех производных в системе кинетических уравнений (1-15). Полученные таким образом алгебраические уравнения определяют стационарные численности клеток п0ст, пшст и пшст. Для исследования устойчивости стационарных состояний изучалось поведение решений системы уравнений (1-15) вблизи значений, отвечающих стационарным численностям, для чего полагалось (16):
п
исл
0
- n0 + У0,
У0 < n0
Пшл - + Ун, Ун < пн, і - 1 - 9. (16)
пТ -пмт + Умі, Умі < пмт, і -1 -5.
В результате получена линеаризованная система уравнений (17-31):
dt
- Уо
гр(1 - 2 По G 0о) - к°(1 -Го) К
(17)
d^ -
dt
У о [сн k°(1 - Г°)] + Ун1
Ги(1 - 2
пн1
- G°H1 К.1
) - к н1(1 -Гн1) ,
(18)
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
7
dy н2 dt
У ні [k н1(1 -a)(1 -Гн1)] + У н2
^н2(1 - 2
пн 2
- G0н2
K н2
) -^н2(1 -?н2) ,
(19)
dt
Ун(і-1) Vен(і-1)
[k н„-1,(1 -Гн(1 -1))] + У н
Гн (1 - 2 п“ г00 ■') - k н, (1 -Гм)
K н,
і = 3,4.
(20-21)
dyH5
dt
Ун 4 [k н4 (1 -Гн4)]- Ун5kн5,
(22)
dyH,
dt
_ yн(і-1^н(і-1)
Унік
ні ’
і = 6 - 9.
(23-26)
dyM1 =
dt
y0 [см k0(1 - Г0)] + Ум1
Гм1(1 - 2
”M1 - 0<Ц - *м1(1 -Гм1)
м1
(27)
dy м2
dt
У м1 [k M1(1 -a)(1 -Гм1)] + y м2
Гм2(1 - 2
пм2 - 00м2
K м2
) - kM2(1 -gM2) ,
(28)
dyM3
dt
yм2 [kм2 (1 -Гм2)] + y м3
Гмз(1 - 2-
nM3 - 00
K
—) -kмз(1 -Гмз)
м3
(29)
dyM4
dt
yм3 [kM3 (1 gM3)] yM4kM4
(30)
Фм5
dt
yM4kM4 yM5kM5
(31)
Далее для построения матрицы Гурвица были введены коэффициенты С00, Сшу и Сщ для КОЕ-ГМ, нейтрофилов и моноцитов соответственно. Данные коэффициенты можно разделить на две части относительно места возникновения: первая группа представляет собой
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
8
влияние, оказываемое на устойчивость предыдущей стадией, а вторая -влияние собственных параметров. Преобразованная система уравнений (17-31) представлена в уравнениях (32-46):
dyp
dt
+ C00 yo = 0
dyHi
dt
+ Сн01у0 + Сн11 ун1
= 0
dyH,
dt
+ Сн(і-1)іун(і-1) + Снііуні = 0,
і = 2 - 9.
dyM1
dt
+ См01 y0 + См11 yM1 _ 0
4уж
dt
+ См(і-1)іум(і-1) + Сііумі = 0,
і = 2 - 5.
(32)
(33)
(34-41)
(42)
(43-46)
В таблицах 4-5 приведены рассчитанные значения для С00, Сшу и Сшу.
Таблица 4 - Значения коэффициентов матрицы Гурвица для
нейтрофилов.
С00 0,145
Сн11 0,498 Сн01 -0,039
Сн22 0,669 Сн12 -0,110
Сн33 0,043 Сн23 -0,480
Сн44 0,362 Сн34 -0,046
Сн55 0,018 Сн45 -0,182
Сн66 0,012 Сн56 -0,018
Сн77 0,014 Сн67 -0,012
Сн88 0,121 Сн78 -0,014
Сн99 0,010 Сн89 -0,121
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
9
Таблица 5 - Значения коэффициентов матрицы Гурвица для
моноцитов.
С00 0,145
См11 0,249 См01 -0,005
См22 0,213 См12 -0,055
См33 0,019 См23 -0,081
См44 0,023 См34 -0,020
См55 0,001 См45 -0,023
Как видно из таблиц, все коэффициенты, представляющие влияние предыдущей фракции на устойчивость, имеют отрицательные значения, в то время как коэффициенты данной стадии больше нуля.
Ниже приведены построенные матрицы Гурвица (47-48):
Сн =
( с с00 0 0 0 0 0 0 0 0 0 N
Сн01 0,11 0 0 0 0 0 0 0 0
0 Сн12 сн22 0 0 0 0 0 0 0
0 0 Сн23 Сн33 0 0 0 0 0 0
0 0 0 Сн34 сн44 0 0 0 0 0
0 0 0 0 Сн45 Сн55 0 0 0 0
0 0 0 0 0 Сн56 сн66 0 0 0
0 0 0 0 0 0 Сн67 Сн77 0 0
0 0 0 0 0 0 0 Сн78 Сн88 0
v 0 0 0 0 0 0 0 0 Сн89 Сн99 у
( C00 0 0 0 0 0 "
См01 см11 0 0 0 0
См =
м
с
м12
0
0
0
с,
0
м22 См23 См33
с
м34
0
0
0
м44
'м45
с
0 0 0
м55 J
(47)
(48)
В таблице 6 приведены главные диагональные миноры матриц и Лмг, рассчитываемые по формуле (49):
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
10
C00 0 0 ... 0
C01 C11 0 ... 0
0 C12 C22 .0 (49)
0 0 0 ••• C
Как видно из (49), коэффициенты, представляющие влияние, оказываемое предыдущей, более молодой фракцией, не играют никакой роли в расчете главных диагональных миноров. Из этого следует вывод о том, устойчивость стационарного решения каждого из уравнений, описывающих гемопоэз, зависит от значения параметров, характеризующих только данную стадию созревания.
Таблица 6 - Главные диагональные миноры матриц Гурвица для
нейтрофилов и моноцитов.
Л0 0,145
Аи1 0,072 Ам1 0,036
Ан2 0,048 Ам2 7,69 ■ 10-3
Ан3 2,077 ■ 10-3 Ам3 1,461 ■ 10-4
Ан4 7,520 ■ 10-4 Ам4 3,361 ■ 10-6
Ан5 1,354 ■ 10-5 Ам5 3,361 ■ 10-9
Ан6 1,624 ■ 10-7 - -
Ан7 2,274 ■ 10-9 - -
Ан8 2,751 ■ 10-10 - -
Ан9 2,751 ■ 10-12 - -
Согласно критерию устойчивости Рауса-Гурвица [8], для того, чтобы все корни характеристического уравнения имели отрицательные действительные части, необходимо и достаточно, чтобы все главные диагональные миноры матрицы Гурвица были положительны при условии С00 > 0, то есть выполняются условия (50):
Ану > 0 Аму > °.
(50)
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
11
Таким образом, стационарное состояние приведенной в начале системы уравнений (1-15) асимптотически устойчиво.
Из (49) видно, что потеря устойчивости системы обуславливается переходом хотя бы одного из коэффициентов Cij в зону отрицательных значений, что возможно только в уравнениях, в которых присутствует деление клеток. Для определения границ устойчивости выведены соотношения (51-58), показывающие условия, при которых происходит нарушение:
п0 - G00 <-у
K Г1 - ^'
к0
у
nrn - G0rn < ‘^
K ■ ( 1-g ^
н7 1 - k . ' Ш
2
V
к
i = 1 - 4,
ні У
пмг - G0m < ‘"Ш
K ( 1 -g ^
мі ^ j А /мі
1 - kмi
2
V
к
і = 1 - 3.
мі у
(51)
(52-55)
(56-58)
Из (51-58) следует, что изменения скорости перехода клеток в следующую стадию, скорости роста и естественной гибели не могут привести к нарушению устойчивости. В таблице 7 приведены пороговые значения остальных параметров, при которых устойчивость теряется.
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
12
Таблица 4 - Пороговые значения параметров модели.
Норма Нарушение
П0 3.83 • 105 < 3.5 • 10y
Пн1 6.4 • 106 < 5.9 • 106
Пн2 1.1 • 107 < 1.0 • 107
Пн3 4.85 • 108 < 4.3 • 108
пн4 1.9 • 109 < 1.8 • 109
Пм1 1.7 • 106 < 1.6 • 106
Пм2 3.7 • 10y7 < 3.4 • 107
пм3 6.0 • 108 < 5.3 • 108
G00 3.2 • 105 > 3.54 • 105
G0^ 5.44 • 106 > 5.9 • 106
G0н2 9.83 • 106 > 1.08 • 107
G0нз 4.12 • 108 > 4.62 • 108
^1- и О 1.6 • 109 > 1.8 • 109
G0мl 1.4 • 106 > 1.6 • 106
G0м2 3.1 • 107 > 3.4 • 107
G0мз 5.1 • 108 > 5.7 • 108
K0 7.4 • 104 > 1.4 • 105
Кн1 1.2 • 106 > 2.5 • 106
Кн2 2.93 • 106 > 6.5 • 106
Кн3 1.5 • 108 > 4.8 • 108
Кн4 4.3 • 108 > 9.0 • 108
Км1 3.4 • 105 > 6.8 • 105
Км2 7.7 • 106 > 1.5 • 107
Км3 1.8 • 108 > 5.9 • 108
Из таблицы 4 видно, что для нарушения требуется практически двукратное увеличение поддерживающей ёмкости среды во всех случаях, составляющее по расчетам 233,2 ± 51,9%.
В связи с этим, основными причинами потери устойчивости являются изменения количества клеток в популяциях n0, G00, пні, О0ні, пмі и О0мі. Расчеты показывают, что пороговое значение в данном случае составляет в среднем 8,6 ± 1,2%.
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
13
Заключение
В статье приведены результаты исследования устойчивости модели нейтрофиломоноцитопоэза.
При помощи критерия Рауса-Гурвица вычислено, что приведенная система дифференциальных уравнений, описывающих созревание клеток, является асимптотически устойчивой. Также установлено, что
устойчивость каждого уравнения для отдельных популяций созревающих клеток зависит исключительно от значения параметров, характеризующих только данную стадию.
Определены пороговые значения параметров модели, при которых система становится неустойчивой. К потере устойчивости приводят лишь значительное увеличение поддерживающей ёмкости среды (в среднем
233,2 ± 51,9%). При этом, неустойчивость вызывают гораздо меньшие колебания численности клеток в n0, G00, пні, G0ni, пмі и GO^ (в среднем 8,6 ± 1,2%), что позволяет сделать вывод о решающей роли этих параметров в устойчивости состояния лейкопоэза.
Список литературы
1. Ciesla В. Hematology in Practice. — Philadelphia : F. A. Davis Company, 2007. — 348 p.
2. Свищенко В. В., Гольдберг Е. Д. Математическое моделирование кинетики эритропоэза. — Томск : Изд-во Том. ун-та, 1995. — 94 с.
3. Шиффман Ф. Дж. Патофизиология крови. — СПб. : БИНОМ - Невский диалект,
2000. — 448 с.
4. Шарай И. А., Тумаев Е. Н. Математическое моделирование кинетики
нейтрофилопоэза // Материалы XI научно-практической конференции молодых ученых и студентов юга России "Медицинская наука и здравоохранение". — Краснодар, 2013. С. 5-17.
5. Шарай И. А., Тумаев Е. Н. Математическое моделирование кинетики
моноцитопоэза // Труды X Всероссийской научной конференции молодых ученых и студентов "Современное состояние и приоритеты развития фундаментальных наук в регионах". — Краснодар : Просвещение-Юг, 2013. С. 58-60.
6. Луковенков А. В. Устойчивость стационарных состояний в кинетике многостадийных химических и биохимических процессов : автореф. дис. ... канд. хим. наук : 02.00.04. — М, 2012. — 24 с.
http://ej.kubagro.ru/2014/05/pdf/34.pdf
Научный журнал КубГАУ, №99(05), 2014 года
14
7. Липунова Е. А., Скоркина М. Ю. Физиология крови. — Белгород : Изд-во БелГУ, 2007. — 324 с.
8. Рабинович М.И., Трубецков Д.И. Введение в теорию колебаний и волн. — НИЦ «Регулярная и хаотическая динамика», 2000. — 560 с.
References
1. Ciesla V. Hematology in Practice. — Philadelphia : F. A. Davis Company, 2007. — 348 p.
2. Svishhenko V. V., Gol'dberg E. D. Matematicheskoe modelirovanie kinetiki jeritropojeza. — Tomsk : Izd-vo Tom. un-ta, 1995. — 94 s.
3. Shiffman F. Dzh. Patofiziologija krovi. — SPb. : BINOM - Nevskij dialekt, 2000. — 448 s.
4. Sharaj I. A., Tumaev E. N. Matematicheskoe modelirovanie kinetiki nejtrofilopojeza // Materialy XI nauchno-prakticheskoj konferencii molodyh uchenyh i studentov juga Rossii "Medicinskaja nauka i zdravoohranenie". — Krasnodar, 2013. S. 5-17.
5. Sharaj I. A., Tumaev E. N. Matematicheskoe modelirovanie kinetiki monocitopojeza // Trudy X Vserossijskoj nauchnoj konferencii molodyh uchenyh i studentov "Sovremennoe sostojanie i prioritety razvitija fundamental'nyh nauk v regionah". — Krasnodar : Prosveshhenie-Jug, 2013. S. 58-60.
6. Lukovenkov A. V. Ustojchivost' stacionarnyh sostojanij v kinetike mnogostadijnyh himicheskih i biohimicheskih processov : avtoref. dis. ... kand. him. nauk : 02.00.04. — M,
2012. — 24 s.
7. Lipunova E. A., Skorkina M. Ju. Fiziologija krovi. — Belgorod : Izd-vo BelGU,
2007. — 324 s.
8. Rabinovich M.I., Trubeckov D.I. Vvedenie v teoriju kolebanij i voln. — NIC «Reguljarnaja i haoticheskaja dinamika», 2000. — 560 s.
http://ej.kubagro.ru/2014/05/pdf/34.pdf
