CC BY
63
УДК 531/534: [57+61]
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕСТАБИЛЬНОСТИ ПОЗВОНОЧНИКА И МЕТОДОВ СТАБИЛИЗАЦИИ
С.В. Орлов1, Р.Л. Седов2, Н.Д. Бобарыкин2, В.И. Аполлинариев2
Институт биомеханики позвоночника и суставов, Россия, 236022, Калининград, Советский пр., 1, e-mail: [email protected]
2 Калининградский государственный технический университет, Россия, 236022, Калининград, Советский пр., 1, e-mail: [email protected], [email protected]
Аннотация. К одному из наиболее важных вопросов биомеханики позвоночника относится нестабильность сегментов кинематической цепи позвоночника. По мнению американского физика Ильи Пригожина (1991), «нестабильность» - это состояние системы, характеризующееся неоднородностью и разновремённостью каждого из протекающих процессов и всех изменений в целом. Это определение справедливо и для механической системы позвоночного столба. В свою очередь, нестабильность позвоночного столба - это мультикаузальная проблема, которая может встречаться в любом возрасте. Причинами этой патологии могут быть как травмы с повреждением опорных структур позвоночника, так и воспалительные деструктивные процессы (остеомиелит, туберкулез), онкологические поражения позвонков, аномалии развития, дегенеративные поражения межпозвоночных дисков или связочного аппарата. Изучение природы нестабильности методом математического моделирования имеет важное прикладное значение, так как предлагает биомеханическое обоснование хирургической коррекции этой патологии, оптимизирует выбор клинических решений для типичных вариантов нестабильности. Это актуально как для врачей различных специальностей (ортопедов-травматологов, нейрохирургов, неврологов, фтизиатров, онкологов), так и для инженеров-технологов, разрабатывающих позвоночные импланты и фиксаторы. При помощи модели изучены основные формы нестабильности позвоночника при травмах, а именно повреждения сгибательного типа, разгибательного типа и ротационные. Предложены биомеханически обоснованные варианты хирургической коррекции типичных вариантов нестабильности, приведены расчёты оптимальных параметров стабилизирующих конструкций. Доказано преимущество применения демпфирующих (полуригидных) конструкций по сравнению с жёсткими.
Ключевые слова: нестабильность позвоночника, математическая модель,
оптимизация методов стабилизации позвоночника, фиксирующие устройства.
Введение
Под нестабильностью позвоночника подразумевают такое нарушение взаимодействия между телами позвонков, когда вследствие сугубо механических причин изменяются нормальные законы статики и кинетики позвоночного столба на определенном участке, что проявляется в избыточной и аномальной подвижности тел, некорректных их перемещениях, выходящих за физиологические пределы (R. Louis, © Орлов С.В., Седов Р.Л., Бобарыкин Н.Д., Аполлинариев В.И., 2010
Орлов Сергей Владимирович, к.мед.н., директор Института биомеханики позвоночника и суставов, Калининград
Седов Роман Леонидович, аспирант, ассистент кафедры прикладной математики, Калининград Бобарыкин Николай Дмитриевич, д.т.н., профессор кафедры прикладной математики, Калининград Аполлинариев Валерий Иванович, к.т.н., заведующий кафедрой прикладной математики, Калининград
09806267
1983, R. Fergusson, 1988). Это может являться причиной многочисленных сугубо медицинских проблем, часто требующих хирургических методов лечения (Ulrich A. Wagner, 1998, Thomas R. Haher, William Tallman Felmly, Michael O'Brien, 1991).
Для правильного решения важной практической задачи - хирургического лечения нестабильности позвоночника необходимым условием является оптимальный подбор корригирующих фиксирующих приспособлений, способных предотвратить некорректные перемещения тел позвонков и одновременно максимально точно протезировать упругие и динамические характеристики связочно-капсульного аппарата позвоночника, что позволяет сделать математическая модель трехпозвонкового комплекса человека. В работе [3] подробно описана математическая модель трёхпозвонкового комплекса человека при фиксации позвонков жесткими пластинами. В настоящей работе предлагается оптимизация применения жестких фиксирующих пластин и показано преимущество гибких или полуригидных фиксаторов.
В основу модели положено математическое описание динамических процессов дифференциальными уравнениями Лагранжа II рода, составленных на основе расчетной схемы трехпозвонкового комплекса, представленного дискретными сосредоточенными массами, связанными вязкоупругими элементами и обладающими определенными геометрическими параметрами [1].
Цель исследования - при помощи математической модели трёхпозвонкового комплекса человека провести расчёт и оптимизацию механических свойств применяемых фиксаторов.
За основу был принят принцип стабильности позвоночного столба, изложенный L. Rene [8], где стабильность позвоночника представлена в вертикальной, горизонтальной и аксиальной плоскостях (ротация), что обеспечивается телами позвонков с дугоотростчатыми суставами, которые связаны между собой упруго-демпфирующими элементами (межпозвоночные диски, мышечно-связочный аппарат).
Механическая система обладает следующими свойствами:
1. Механическая система является диссипативной.
2. Распределение нагрузок соответствует трехстолбовой концепции F. Denis [6].
3. Предел прочности тел позвонков и упругодемпфирующих элементов, а также их упругая деформация и плотность считались условно установленными по данным работы [1].
4. Изменение геометрических характеристик трехпозвонкового комплекса соответствовало типичным типам статико-динамических нарушений стабильности позвоночника [2, 6].
Конструктивная трёхпозвонковая модель
Расчетная схема фрагмента позвоночника человека, состоящая из трех позвонков с клиновидным средним позвонком и стабилизирующими конструкциями, представлена на рис. 1 с вариантом клиновидной деформации среднего позвонка и двумя стабилизирующими конструкциями (для передних и задних опорных комплексов).
Здесь введены следующие обозначения:
J-, Ми X - момент инерции, масса, координата i-го позвонка (i = 1, 2, 3);
Ccmj - коэффициенты жесткости j-й стабилизирующей пластины (j = 1, 2);
Copj - коэффициенты жесткости j-й опоры (j = 1, 2);
di = 25 мм; d2 = 20 мм; d3 = 5 мм; d4 = 30 мм;
d[ = 8,5 мм; d2 = 26,5 мм; d'i = 17 мм;
l1 = 32 мм; l2 = 10 мм; l3 = 25 мм; l4 = 10 мм; l5 = 50 мм;
m1 = m2 = m3 =0,1 кг;
б
1
с,
т
11
й\
2
О 12
Jl, М\, Х\, ах
XI
—
с
Ссш\
С4
О о
^, Мз, Хз, аз
Х5
Сорх
$3
IV
$4
$3
>< "Л
Х12 Хх3
с; ^ ■<. сз
2 3
О о
J2, М2, Х2, а2
Х3 Х22
х23
с < < с
с5 .„5* ,--С6
Сст2
О
Хб
Сор2
У
X
I
I
3
4
0
Рис. 1. Расчетная схема трехпозвонкового комплекса с патологией среднего позвонка и
ее двухсторонней стабилизацией
/х = J2 = J3 = 35 кг мм2; Уц = 17,2 мм; £ = 13,25 см2.
Координаты Хц и Уц центра тяжести и момент инерции J плоского позвонка определяются формулами:
Е Х,-У Е X, - У ХУ
Хц = ^------; Уц = ^------------; J = 2ТЦ(X’- + У'-)ШУ, (1)
Е у, Е х, 00
,=1 ,=1
0 ми Уц „ Ш2 У
Ш1 /~ч ^
; V \
^/2 Ч ? dl d2 Г 1 ч / dз+d4 г
01
02
Рис. 2. Распределение сил по оси х
где (X,, У,) - координаты центров тяжести элементарных фигур, составляющих плоский позвоночник, , = 1, 2, 3, 4; у - удельная поверхностная масса позвонка, кг/мм2.
Для упрощения динамической модели трехпозвонкового комплекса коэффициенты жесткости С2 и С3, а также коэффициенты жесткости С5 и С'6 приведены к одним коэффициентам С2 и С4 (рис. 2), соответственно, по следующим формулам:
С _ С^2С'2 + (<^2 + ^ X 3 . С _ <^2С2 + (<^2 + <^3)С2
2^2 + dз 2й2 + dз
(2)
На расчетной схеме (см. рис. 2) третий позвонок связывается посредством жестких элементов Сор1 и Сор2 с опорой по оси ОХ, а первый - по оси У через Су.
Для фиксации вариантов нестабильности позвоночника предусмотрено применение условных жестких плоскостных конструкций с коэффициентами жесткости Сст1 и Сст2, что позволяет моделировать как жесткие ригидные металлические системы, как и полуригидные пружинные элементы.
Математическая модель трёхпозвонкового комплекса
Разработанная математическая модель позволяет на основе вычисления внутренних нагрузок опорных комплексов каждого позвонка трехпозвонкового комплекса рассчитывать варианты переломов и нестабильности позвонков в различных зонах при их патологии. Кроме этого, можно произвести расчет смещения позвонков по оси Оу под воздействием силы 02у, что чаще всего является причиной стеноза позвоночного канала и может приводить к сдавливанию дуального мешка. Выбранная динамическая модель трехпозвонкового комплекса человека (см. рис. 2) является механической системой, для которой уравнения Лагранжа II рода имеет вид
{дТ Л —
0, к = 1, 2,..., 7, (3)
dt
дх,
дП дФ
+-----+ -
дх, дх, дх.
где Т, П - кинетическая и потенциальная энергия системы;
Ф - диссипативная функция, определяемая спинными мышцами;
0к - внешние воздействия.
В качестве обобщенных координат хк принимаются следующие координаты: х1, х2 - абсциссы центров тяжести тела 1-го позвонка и его остистого отростка; х3, х4 - абсциссы центров тяжести тела 2-го позвонка и его остистого отростка; х5, х6 - абсциссы центров тяжести тела 3-го позвонка и его остистого отростка; х7 = у - ордината центра тяжести тела 1 -го позвонка.
у
Проекции сил Q1 и Q3 на оси Х и У при наличии деформации среднего позвонка трехпозвонкового комплекса человека определяются формулами
Qlx = О2 Ql2 = О1 Ql2SІn Р1 = QlSІn2 Р1;
Qly = О1 О2 = Qlsin Р1^ Р1;
Qзx = О4 QзlSІn Рэ = Qзsin2 Рз; (4)
Qзy = Оз О4 = Qзsin Pзcos Рз;
Q2y = Qly + Qзy = СХ^т PlCOSPl + СзХ^Ш PзCOS Рз,
где Q12, Q31 - модули проекций сил на ось X, Ог- - вспомогательные точки, обозначающие границы векторов.
Кинетическая энергия механической системы трехпозвонкового комплекса человека, приведенная на рис. 2, равна
( dX )2 ( da )2
где (~^~) , (~ci~) - кваДРаты скоростей колебаний и вращений 1-го, 2-го и 3-го
позвонков относительно центров тяжестей этих позвонков.
Упругие деформации х1, х2, х3, х4, х5, х6 центрального и правого столбов 1-го,
2-го и 3-го позвонков связаны со смещениями и вращениями центров тяжестей этих
позвонков X1, Х2, Х3 и а1, а2, а3 следующими соотношениями:
X = Xi - aid1; xi = Xi + aid2; i = 1, 2, 3 . (6)
Задача решается в приближении малых смещений, то есть
Xi << Di; tg ai = sin ai = ai. (7)
Тогда X1 и a1 выражаются через х1 и х2 по формулам (для второго и третьего
позвонков аналогично)
X2 - X v X2dj + xxd2
a = —-----1; X, =-^л-----‘-^. (8)
1 D1 1 D1
С учетом соотношения (8) кинетическая энергия трехпозвонкового комплекса с патологией среднего позвонка и его двухсторонней стабилизацией, приведенного на рис. 2, запишется в следующем виде:
Т = ^ + m„ xx, + ^ ^ + m2„ x Х4 +
+ тз0 X5 x6 +
2 О
L Щ (y)
i=1
(9)
2
2
2
(т Л даА2 + 4 . (т Л тД2 + 4 . (т Л ЩёА -
ГДе ИХ = ------Д-----. КХ = Д----------’ = -----Д2----'
Потенциальная энергия механической системы трехпозвонкового комплекса человека считается равной нулю при положении статического равновесия, а отсчет деформации упругих элементов ведется от условия, когда статическая нагрузка на элемент уравновешивается упругой силой от его осадки.
В этом случае потенциальная энергия П деформации упругих элементов трехпозвонкового комплекса с патологией среднего позвонка и ее двухсторонней стабилизацией определяется следующим соотношением:
П = СоР1 х52 + Сор2х62 + Сз(Х3 -х5)2 + С4(х4 -х6)2 +
2
2
2
2
+ С1(х1 - хз)2 + С2(Х2 - х4)2 + Сст1(х3 - х5)2 + Сст2(Х4 - Х6)2 + СУ^
2
2
2
2
2
(10)
Диссипативная функция Ф трехпозвонкового комплекса человека с патологией среднего позвонка и его двухсторонней стабилизацией записывается через коэффициенты демпфирования В] как
Ф = ВоР1(х5)2 + ВоР2 (х6)2 + В3(хз - х5 )2 + В4 (х4 - ^^6)2 +
+ В1(^ - хз)2 + В2(х2 - х4)2 + Вст1(-У3 - х5)2 + Вст2(хА
)2 , ву (у)2
2
2
2
2
- +
2
(11)
Подставляя значения производных от кинетической и потенциальной энергии, а также от диссипативной функции Ф для механической системы трехпозвонкового комплекса (9)-(11) в уравнения Лагранжа II рода, получим в векторном виде
С =
М ё2 X &2 + Б- — + СХ & = О,
т1с1 т10 0 0 0 00 “
т10с1 т2 0 0 0 00
0 0 т3с3 т20 0 00
0 0 т20с3 т4 0 00 ?
0 0 0 0 т5 т30 0
0 0 0 0 т30 т6 0
0 0 0 0 0 0 М1 + М 2 +М3
С1с 10 - 1 3 0 Сст1 0 0
0 С2 0 -С2 0 0 0
1 1 1 0 (С1 + С3)с2 + Сст1 0 -Сст1 - С3 0 0
0 -С2 0 С2 + Сст2 + С4 0 -(Сст2 + С4) 0
0 0 -Сст1 - С3 0 Сор1 + Сст1 0 0
0 0 0 -Сст2 - С4 0 Сор2 + Сст2 0
*1 0 Sз 0 0 0 Су
(12)
2
2
2
2
6
B =
B1c1 0 -B1C3 0 Bcm1 0 0
0 B2 0 -B2 0 0 0
- B1c1 0 (B1 + B3)c2 + Bcm1 0 -Bcm1 - C3 0 0
0 - B2 0 B2 + Bcm2 + B4 0 -(Bcm2 + B4) 0B
0 0 -Bcm1 - B3 0 Bop1 + Bcm1 0 0
0 0 0 -Bcm2 - B4 0 Bop2 + Bcm2 0
B1 0 3 B3 ^0 0 0 0 B У
Q = [a Q2 0 0 0 0 Q-y I; х = [ x x2 x3 xA Х5 Хб у f,
где с1 = cosPj, c3 = cosP3, S1 = CjXjSinPjcosPj, S3 = C3x3sinP3cosP3, SB1 = B^sinPjcosPj,
SB3 = B3x3sin P3cos P3.
Умножая обе части векторного уравнения (12) на обратную матрицу М^1, получим следующее приведенное векторное уравнение:
d-X + M-1 -B-— + М-1 CX = М-1 Q dt2 dt
или
d X dX . , ч
—- + W------------------------------------------+ А-X = F, (13)
dt2 dt
где А = М 1 С; W = M-1 В; F = М-1-Q .
Метод численного решения
Система обыкновенных дифференциальных уравнений второго порядка с постоянными коэффициентами сводится к системе 14 обыкновенных дифференциальных уравнений первого порядка в векторной форме, которая имеет следующий вид:
dX dZ
------= Z; — + W Z = F - A X.
dt dt
(14)
При этом система дифференциальных уравнений 1-го порядка (14) решается по разностным схемам второго порядка точности по времени 0(т2):
х ,+1 - х ,
--1------ + W•Z = Г -А-Х •
т ™ -+1 -+1 -
Хм—X-=• -=0,1,..., *
т
2
при начальных условиях Х0 = Z0 = 0.
Алгоритм численного решения системы векторных уравнений (14) записывается
в виде
Z;+1 = (1 + tW) - (Z, + t(F;+1 - A - X,)); X j+1 = X j+2 (Z j+1+Z j); j = 0,1,..., k -1
(15)
Оптимизация
Цель принятия решения - устранение нестабильности позвоночного столба в очаге поражения - трёхпозвонковом комплексе человека с клиновидной деформацией среднего позвонка при помощи фиксации второго и третьего позвонков функциональными пластинами так, чтобы нагрузки на страдающий участок кинематической цепи позвоночника максимально приблизились к физиологической норме. Для этого задаётся поиск минимума отклонения нагрузки Р с патологией от нагрузки нормальной Р„ОГт с сохранением подвижности всей системы.
Критерий оптимальности:
|Р5(0, Се1л, Х1, Х3, р1, рз) -Р5пОпт\ ® Ш1п, Рб(0, Сег2, Х„ Xз, Р1, Рз) -Р6пОпт\® Ш1п,
(16)
где Аг = |р(0, Саг, Х1, Х3, Р1, Рз)-РГт\, (/ = 5,6) - компоненты двумерного вектора А
целевых функций. В силу положительности целевых функций тривиальным решением является нулевое решение. В этом случае нагрузки Рг- тождественно приблизились к константам Р"Огт, поэтому вырожденное решение является физиологически
допустимым. Целевые функции также зависят от параметров С и С3 - коэффициентов жесткости межпозвоночных дисков. Эти параметры выбираются для каждой отдельной практической задачи частным образом [5].
В качестве локальной задачи оптимизации выберем нахождение параметров жесткости левой Сс^ и правой Сс^2 полуригидных фиксирующих устройств. Таким образом, задача исследования операций в своей постановке полностью отражает структуру и динамику знаний принимающего решение субъекта о множестве допустимых решений и о критерии оптимальности [5]. Приведём результаты расчётов оптимальных параметров пластин. Векторный критерий оптимальности сведён к минимизации следующей целевой функции:
Рис. 4. Стабилизация нагрузки Р5 (сплошная линия) при оптимизации жесткости пластины 1 при смещении нагрузки на у = 0, Р"оп" - физиологическая нагрузка в норме
(пунктирная линия)
^ = а1 р5 - Р5* | + а2 р6 - Р6* | ® шт, (17)
где аг - коэффициенты предпочтения /-го критерия, г = 1, 2.
Численный метод оптимизации основан на параллельных вычислениях массивов
параметров оптимизации С&1 и Сс2. На встроенном программном модуле в пакете
ЫмЪСАП построена программа, реализующая математическую модель нестабильности
позвоночника человека. Процедура урта\ро2{х1, х2, ..., хп) выполняет нанесение
равномерной сетки на плоскость С&. Вторая процедура рогуои(х1, х2, ..., х„, Сс^, Сс^2)
выполняет численное решение разрешающей системы дифференциальных уравнений в
каждом узле сетки плоскости С&. Задавались следующие исходные данные [10]:
С1 = 3,26-103 Н/мм; С2 = 0,92-103 Н/мм; С3 = 0,46-103 Н/мм;
С4 = 3,26-103 Н/мм; С5 = 0,92-103 Н/мм; С6 = 0,46-103 Н/мм; Су = 5-103 Н/мм;
Сор1 = Сор2 = 3,26-103 Н/мм; Q = 800 Н (внешняя сила приложена к центру тяжести
позвонка при у = 21 мм); р1 = 0,175; р3 = 0,262 (углы деформации среднего позвонка);
ц = 1 мм (предел подвижности позвоночника); число временных слоев N = 1000, шаг
-2
интегрирования по времени т = 10 с [5].
Как видно из рис. 4, поиск оптимальных жесткостей пластин сопровождается приближением механической системы к нормальному состоянию, близкому к показателям модели при отсутствии деформаций (у здорового человека) [5].
В результате вычислений получен интервал допустимых жесткостей пластин, применяемых для стабилизации позвоночника человека:
350 < Сстх < 10000 Н/мм,
10 < Сст2 < 2010 Н/мм.
Заключение
По данным математического моделирования, в условиях нестабильного положения позвоночника, связанного с разрушением межпозвонкового диска или тела позвонка, в том числе с изменением его упруго-прочностных свойств или геометрических параметров, оптимальным вариантом его стабилизации является применение фиксирующих конструкций упругого типа (например транспедикулярных
систем с пружинными штангами из NI-TI-сплава или системы «Dynesis»). Это позволяет сохранить распределение эпюр нагрузок на тела позвонков в пределах значений, близких к природным - 0,7 / 0,3 [3].
Выводы
1. При нестабильности позвоночника более физиологично применение динамических фиксирующих систем.
2. Применение динамических фиксирующих систем при межтеловой стабилизации является профилактикой раннего износа смежных дисков в системе [9].
Список литературы
1. Громов А.П. Биомеханика травмы. - М.: Медицина, 1979. - С .17 9-210.
2. Давыдов Е.А., Загородний Н.В., Древаль О.Н. Основные направления разработки и использования конструкций с эффектом памяти формы для динамической стабилизации позвоночника // Современные технологии в травматологии и ортопедии: материалы III Междунар. конгр. - М., 2006. - С. 241.
3. Орлов С.В., Бобарыкин Н.Д., Латышев К.С. Математическая модель стабильности трёхпозвонкового комплекса // Математическое моделирование. - 2006. - Т. 18, № 10. - С. 55-70.
4. Орлов С.В., Каныкин А.Ю., Москалёв В.П., Щедрёнок В.В., Седов Р.Л. Математический расчёт прочности позвоночного столба при хирургическом лечении нестабильных переломов позвоночника // Вестник хирургии имени И.И. Грекова. - 2009. - Т. 168, № 2. - С. 61-64.
5. Седов Р.Л., Орлов С.В., Бобарыкин Н.Д. О расчёте параметров динамических стабилизирующих конструкций на основе математической модели трёхпозвонкового комплекса человека // Математическое моделирование. - 2010. - Т. 22, № 2. - С. 113-123.
6. Denis F. Spinal instability as defined by the three column spine concept in acute spinal trauma // Clin. Orthop. - 1984. - Vol. 189. - P. 65.
7. Fergusson R., Tencer A., Woodard P., Allen A. Biomechanical comparison of spinal fracture models and the stabilizing effects of posterior instrumentations // Spine. - 1988. - Vol. 13. - P. 453.
8. Reno L. Surgery of the Spine. - Berlin-Heidelberg: Springer-Verlag, 1983. - P. 55-58.
9. Haher T.R., Felmly W.T., O’Brien M. Thoracic and lumbar fractures: diagnosis and management // Spinal Surgery. - 1991. - Vol. 2. - P. 857-910.
10. Wagner U.A. Unter mitarbeit von Ottmar Schmitt, Hans-Martin Schmidt, Wallny T. Atlas der Pedikelschraubenimplantate. - Georg Thieme Verlag Stuttgart, 1998. - P.1-2.
MATHEMATICAL MODELING INSTABILITY OF THE SPINE AND STABILIZATION METHOD
S.V. Orlov, R.L. Sedov, N.D. Bobarikin, V.I. Apollinariev (Kaliningrad, Russia)
The one of most important issues in biomechanics is the instability of segments of the kinematic chain of the backbone. According to American physicist Ilya Prigogine (1991), "instability" is a state of the system characterized by heterogeneity and the difference in each of the processes taking place and all the changes in general. This definition is valid for the mechanical system of the spine. In turn, the instability of the spine is a many reasons problem and can occur at any age. The causes of this pathology may be as injuries with damage to supporting structures of the spine, so the destructive inflammatory processes (osteomyelitis, tuberculosis), cancer, vertebral lesions, developmental anomalies, degenerative lesions of intervertebral discs or ligaments. The study of the instability by the method of mathematical modelling is of great practical importance, since offers a biomechanical substantiation for surgical correction of this pathology, optimizes the selection of clinical decisions for the typical variants of instability. This is important for doctors of various specialties (orthopedists-traumatologists, neurosurgeons, neurologists, phthisiologists, oncologists), and for the engineers developing the spinal implants and fixators. Using the model, the basic
forms of instability at spinal trauma, namely damage of the flexor type, extensor type and rotary type are studied. Biomechanically-based options are offered for surgical correction of the typical variants of instability. Calculations of optimum parameters of the stabilizing structures are described. The advantage of damping (semi rigid) structures in comparison with rigid is proved.
Key words: spinal instability, mathematical model, optimization of methods of spinal stabilization, fixing devices.
Получено 31 мая 2010
