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

DOI: 10.15593/RZhBiomeh/2016.3.07 УДК 53.083.92

к Российский

Журнал / Биомеханики

www.biomech.ru

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТРЕХПЛОСКОСТНОЙ ДЕФОРМАЦИИ ПОЗВОНОЧНОГО СТОЛБА ЧЕЛОВЕКА

1 Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, Россия, 197101, Санкт-Петербург, проспект Кронверкский, 49, e-mail: culter88@yandex.ru

2 Восстановительный центр детской ортопедии и травматологии «Огонек», Россия, 198515, Санкт-Петербург, пос. Стрельна, Санкт-Петербургское шоссе, 101, e-mail: ogonek@zdrav.spb.ru

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

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

Одним из тяжелейших заболеваний человека прямоходящего (Homo erectus) является трехплоскостная деформация позвоночного столба (сколиоз). В теории и практике есть только два бесспорных факта о данном заболевании - оно проявляется всегда в виде трехплоскостной деформации, возникновение и развитие которой связано с процессом роста [8, 10]. И это при том, что в соответствующей литературе в качестве его причины рассматривается и обсуждается большое количество так называемых «этиологических» факторов [3, 7].

В литературе уже имеется ряд работ [5, 6, 8], в которых на основе различных математических подходов описываются анатомические параметры и характеристики деформированного позвоночного столба с целью решения самой востребованной в клинической практике задачи по объективному прогнозированию развития сколиоза. Однако в указанных работах имеется один существенный недостаток - в их

© Дудин М.Г., Балошин Ю.А., Бобер С.В., Поморцев И.Ю., 2016

Дудин Михаил Георгиевич, д.м.н., главный врач, директор центра, Санкт-Петербург

Балошин Юрий Александрович, д.т.н., профессор кафедры нанофотоники и метаматериалов,

Санкт-Петербург

Бобер Степан Васильевич, инструктор-методист лечебно-физической культуры, Санкт-Петербург Поморцев Иван Юрьевич, сотрудник центра, Санкт-Петербург

М.Г. Дудин2, Ю.А. Балошин1, С.В. Бобер2, И.Ю. Поморцев2

2

Введение

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

C точки зрения классической биомеханики нормальный (здоровый, не имеющий деформаций) позвоночный комплекс рассматривается как двухколонное образование [4], состоящее из несущей колонны (24 тела позвонков + 23 межпозвонковых диска) и колонны функциональной (позвонковые дужки, остистые и поперечные отростки, межпозвонковые, или дугоотростчатые, суставы). Первая, несущая (опорная), решает его главную задачу по поддержанию анатомического положения и веса всех органов и систем, входящих в состав туловища, а вторая, функциональная, обеспечивает движения в многозвенном позвоночном комплексе. Последовательно расположенные позвонковые дужки формируют позвоночный канал, в котором располагается спинной мозг. Кроме этого, состояние нормального позвоночного столба определяют граничные условия его существования: 1) вертикальное положение комплекса (колонн); 2) верхние и нижние концы колонн всегда на одном уровне; 3) расстояние между колоннами на любом уровне - величина постоянная.

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

Модель 1 и описывающее ее уравнение

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

Ниже (рис. 1) представлена взятая из работы [2] схема образования трехплоскостной деформации в двухколонной модели позвоночного столба. В условия задачи, описывающей сколиоз как процесс деформирования модели, были добавлены определенные требования к отрезкам. Дорсальный прямой (S, short) рассматривается как отрезок, не меняющий свой продольный размер O1O3, а

вентральный (L, long) представлен в виде кривой ABA, обладающей свойством менять продольный размер на некоторую величину AL, но при этом любая точка данной кривой равноудалена от отрезка OO на постоянное расстояние r.

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

Для решения задачи рассматривается развертка поверхности цилиндра радиусом r и высотой S, которая содержит кривую L, симметричную относительно оси O2B2, поэтому достаточно найти длину ее половины и удвоить.

Аппроксимируем кривую L полиномом третьей степени:

y(x) = ax3 + bx2 + cx + d . (1)

Рис. 1. Схема образования трехплоскостной деформации позвоночника [2] По формуле нахождения длины кривой, заданной уравнением y(x)

L = 2 jyj 1 + y'2 dx.

(2)

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

у(0) = ё = гф,

у^у'^ (3)

У(0) = С.

Подставив условия (3) в (1) и подсчитав производную, получаем

у(х) = ах3 + Ьх2 + Гф; у'(х) = Ъах2 + 2Ьх. (4)

Далее, подставив равенства из системы (3) в полином и получив решение системы уравнений, имеем значения коэффициентов полинома

16гф 12гф

а = —г1; Ь =--т-1. (5)

Я3 Я2

Производя несложную замену переменных, получаем производную в следующем виде

125

ô = гф ; 4k = ■

S

е 2 x 5=:5х,

(6)

0

= k (е Ч); = , x е

° f

* = f !] • (7)

5 '

Применив данные преобразования к формуле расчета длины кривой (2), получим

— = 5 +к (£2 —£)2 а £. (8)

0

В полученной формуле Ь - длина колонны, меняющей свой продольный размер, £ - длина колонны, не меняющей свой продольный размер; г - радиус вращения, расстояние между Ь и 8 (постоянная величина на любом горизонтальном уровне модели); ф - наибольший угол поворота, наблюдаемый при скручивании Ь вокруг 8. Также необходимо обратить внимание на переменную х, по которой производится интегрирование в формуле длины кривой; значения этой переменной являются координатами по оси абсцисс точек кривой Ь на графике развертки этой кривой

2 х

(рис. 1). Для удобства расчетов данная переменная была заменена на £ = — .

Математическое моделирование, представленное в данной работе, имеет своей целью дальнейшее развитие идей работы [2] для более качественного объяснения закономерностей образования трехплоскостной деформации позвоночного столба.

Прежде всего был произведен численный анализ формулы (8), который дал возможность в прикладных целях получать значение угла ф в зависимости от значений остальных параметров системы (см. рис. 1).

Интеграл в формуле (8), вероятнее всего, невозможно вычислить аналитически.

Однако, используя численные методы, можно анализировать этот интеграл как

—

некоторую зависимость — = Г (к) = + к (£2 — £) а£ . Если быть точнее,

5 0

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

В частности, при к е [0,40] график Г2 (к) имеет вид, представленный на

рис. 2. После аппроксимации на данном интервале значения коэффициентов получились такими: а = 0,031092491; Ъ = 1,000022401.

Общий вид графика изменения погрешности АГ2 (к) = |г 2 (к) — (аак + Ь)|

представлен на рис. 3. Как можно заметить на рис. 3, при к = 16 достигается наибольшее отклонение, составляющее примерно 15 %.

Далее, получив неплохую аппроксимацию для Г2 (к) , возвращаемся к Г (к). Предполагаем, что, если квадрат функции можно неплохо аппроксимировать линейной функцией, вероятно, саму функцию можно неплохо аппроксимировать функцией

квадратного корня: Г (к) = Vак + Ь = л/а^к + Ь = а^/к + р . Вновь применяя численный

метод на интервале к е [0, 40], получаем такие значения коэффициентов:

а = 0,1763306; р = 32,1628267. (9)

F2 (к)

40 к Рис. 2. Функция Р2 (к) в интервале к е [0, 40]

AF 2( к )

0,15

0 16 40 к

Рис. 3. Изменение погрешности АР2 (к) в интервале к е [0, 40]

При этом погрешность этой аппроксимации АР (к) = |р (к) - (ак + Ь)| не превышает 6,3 %.

Таким образом, мы аппроксимировали с погрешностью, не превышающей —

6,3 %, функцию Р(к) = — функцией а^к + р на интервале к е [0, 40] (сами значения

Я

функции на этом промежутке были в диапазоне от 1 до 1,5 - примерно этот диапазон

соотношения Ь и £ и интересен с точки зрения медицины).

—

Далее, используя равенство — = а^к + р и замены (6), после несложных

Я

преобразований можно записать значение ф в следующем виде:

ф = —^Л2 -а2Я2р2 12а-г ,

(10)

при этом а и р являются некоторыми вычисленными константами, которые, в свою очередь, зависят от отношения — и в общем случае находятся в диапазоне

а«0,175-0,185; р«30-32.

Благодаря численному анализу формулы (8) удалось также построить графики зависимости разности АЬ длины вентрального отдела Ь и дорсального отдела Я от угла поворота ф при условном расстоянии между отделами г в 2,5 см (рис. 4).

2

1

АЬ, %

Рис. 4. Зависимость разности длины АЬ вентрального отдела Ь и дорсального отдела от угла поворота ф при условном расстоянии между отделами г в 2,5 см

Идентичность картин в математической модели и в реальном позвоночном столбе дает основание предполагать идентичность механизма их инициации [3, 4]. Поэтому появляются основания считать, что, продолжая исследование математической модели, можно будет найти те закономерности в развитии реального сколиоза, которые в экспериментах и опытах на животных не обнаруживались [1].

Уравнение и модель 2

Однако в реальном позвоночном столбе изначально имеются физиологические сагиттальные изгибы - грудной кифоз (средняя величина примерно 165-170°) и поясничный лордоз (примерно 150-170°) [4], на протяжении которых можно видеть изначальное различие (± АЬ) между продольными размерами длин вышеуказанных отрезков (колонн) Ь и 5. В грудном отделе это -АЬ, а в поясничном +АЬ.

На рис. 5 представлена модель реального позвоночного комплекса в сагиттальной плоскости. Она описывается двумя параллельными кривыми 5-образной формы с постоянным расстоянием Я между ними. Вентральный отрезок Ь состоит из Ькиф+Ьлрд, а дорсальный 5 из 5киф + 5лрд. При этом повторим, что Ькиф - 5киф = -АЬкиф, а Ьлрд - 5лрд = +АЬлрд. Каждая из сагиттальных прямых («кифоз» и «лордоз») имеют радиус Якиф и ^лрд.

Таким образом, для описания процесса в реальном позвоночном столбе уравнением (1) должно произойти устранение -АЬкиф и +АЬлрд. Это явление и можно наблюдать у реальных больных. Речь идет о симптомокомплексе так называемой «плоской спины».

Переход модели с заложенными в нее физиологическими сагиттальными изгибами, представленной на рис. 5, в строго прямое состояние (см. модель на рис. 1) математически может быть описан уравнением (9) с учетом уравнения (8) как

Рис. 5. Модель реального позвоночного комплекса в сагиттальной плоскости

N

AL = £ф,

R2 + r2 +

2R ■ r2

Vr2 + r 'Чф,)

- r

(11)

Охарактеризуем кратко математическую модель, в рамках которой получена формула (11). На рис. 5 показаны две кривые Ь и Я в полярных координатах. Рассматривается разбиение этих кривых точками , р). Отрезки, которые

образованы каждой из пар (Q, р ), имеют длину г, они перпендикулярны касательным к кривым Ь и Я в точках Q и р. Предполагается, что есть возможность интерполировать кривую Я функцией г(ф), и разбиение будет производиться по достаточно малым промежуткам. Рассмотрим некоторый сегмент OPQjPj+1Qj+1 (рис. 6). По теореме косинусов

R2 = R + r2 - 2R ■ r ■ cos(90 + 0г) = R2 + r2 + 2R ■ r ■ sin 0г.

(12)

Используем простые тригонометрические соотношения

sin 0,. =

1

'1 + е1в20г'

подставим их в формулу (12) и получим

'М ф=ф = r" ct§0i ^ 0 = arcct§

ЛФ,)

R = IR2 + r2 +

2R ■ r2

' 2 + г'2(фг) '

N

Далее, Ь = ^ рф ^ АЬ = Ь - Б - получаем формулу (11).

(13)

,=0

Вид по стрелке А

Рис. 6. Пояснение к математической модели (11)

Уравнение (11) помогает решить задачу о том, на какую величину АЬ необходимо удлинить часть вентрального отрезка Ькиф в грудном кифозе и уменьшить соответствующую часть в поясничном лордозе Ьлрд, чтобы позвоночный комплекс стал прямым.

На рис. 7 представлен результат численного решения уравнения (10), устанавливающего взаимозависимость угла наклона модели позвоночного комплекса от изменения соотношения между вентральной и дорсальной колонами. На рис. 7 видно, как при постепенном увеличении продольной длины меньшей колонны модель, представленная на рис. 5, стремится в сагиттальной плоскости к условному выравниванию Ь = 5, приближаясь к вертикальному положению, а далее на небольшом отрезке значений (при достаточно малых величинах АЬ находится в пределах от -0,125 до +0,125 см при угле ф ~ 1,65 рад) достигает его (пик графика, заключенный в прямоугольник). После этого модель начинает скручиваться в вентрально-дорсальном либо дорсально-вентральном направлении, достигая максимума в крайних значениях графика, представленного на рис. 8.

Заключение

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

2

-1 0 1

АЬ, см

Рис. 7. Зависимость угла наклона ф позвоночного комплекса от увеличения продольной длины АЬ короткой дуги

Рис. 8. Динамика изменений конфигурации модели позвоночного комплекса относительно графика на рис. 4

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

Однако в реальном позвоночном столбе имеются анатомо-физиологические особенности, которые противостоят прямому применению такого вывода, поскольку

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

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

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

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

1. Гайворонский Г.И. Экспериментальный сколиоз: автореф. дис. ... д-ра мед. наук. - Л., 1982. - 36 с.

2. Дудин М.Г., Синицкий Ю.Ф. О механогенезе торсионных изменений при сколиозе // Ортопедия, травматология и протезирование. - 1981. - № 2. - С. 33-36.

3. Дудин М.Г., Михайловский М.В., Садовой М.А., Пинчук Д.Ю., Фомичев Н.Г. Идиопатический сколиоз: кто виноват и что делать? // Хирургия позвоночника. - 2014. - № 2. - С. 8-20.

4. Капанджи А.И. Позвоночник. Физиология суставов. - М.: Эксмо, 2009.

5. Чигрик Н.Н. Геометрическое моделирование многопараметрических процессов сколиотических деформаций позвоночника с целью создания системы диагностики и прогнозирования: автореф. дис. ... канд. мед. наук. - Омск, 2002.

6. Abedrabbo G., Fisette P., Absil P.A., Mahaudens P., Detrembleur C., Raison M., Banse X., Aubin C.E., Mousny M. A multibody-based approach to the computation of spine intervertebral motions in scoliotic patients // Research into Spinal Deformities 8. Studies in Health Technology and Informatics. - 2012. -Vol. 176. - P. 95-98.

7. Bagnall K. How can we achieve success in understanding the aetiology of AIS? // Study. Health Technology. Inform. - 2008. - № 135. - P. 61-74.

8. Duval-Beaupere G., Dubousset J., Queneau P., Grossiord A. Pour unetheorie unique de l'evolution des scoliosis // Presse Med. - 1970. - № 78 (25). - P. 1141-1146.

9. Kumar B., Bylski-Austrow D.I., Liu Y. Finite element model of spinal hemiepi physiodesis: Effect of contact conditions, initial conditions, and growth // Research into Spinal Deformities 8. Studies in Health Technology and Informatics. - 2012. - Vol. 176. - P. 99-103.

10. Roth M. Idiopathic scoliosis caused by short spinal cord // Acta Radiol. Diagnosis. - 1968. - Vol. 7. -P. 257-271.

MATHEMATICAL MODELLING OF 3D HUMAN SPINE COLUMN DEFORMATION

M.G. Dudin, Yu.A. Baloshin, S.V. Bober, I.Yu. Pomortsev (St. Petersburg, Russia)

The mathematical description of the genesis of a 3D deformation in double-column model, approximate to a real spine column of the person is given. This deformation is regarded as a result of two independent processes, in first of which the length surplus of one

of the columns "is absorbed" by the sagittal curves introduced in the model, and in the second the development of 3D ("scoliotic") deformation at exhaustion of this reserve begins. Each of these processes is described by its own equation. The stage-by-stage solution of these equations allows us to define dynamics of changes of a vertebral complex configuration. The presented mathematical description of process of 3D deformation of double-column model of the spine allows us to point out the key links of pathogenesis and of the "target", the influence on which can lead to success in opposition to scoliosis.

Key words: mathematical modelling, deformation, 3D deformation of the spine, mathematical description of the deformation process, double-column model.

Получено 12 апреля 2016