ПРИКЛАДНАЯ МАТЕМАТИКА И МЕХАНИКА
Том 79. Вып. 2, 2015
УДК 539.3 : 51.76
© 2015 г. Л. Б. Маслов
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВОССТАНОВЛЕНИЯ МЕХАНИЧЕСКИХ СВОЙСТВ КОСТНОЙ МОЗОЛИ
Представлена математическая модель и вычислительный алгоритм регенерации костной ткани, управляемой законом дифференциации клеток и действием внешнего механического стимула периодического характера. В основе расчета восстановления упругих свойств костной ткани лежит обобщенная динамическая модель изменяющейся пороупругой сплошной среды и метод конечных элементов в трехмерной постановке. Разработанное программное обеспечение дает возможность исследовать процессы восстановления поврежденных костных элементов опорно-двигательного аппарата человека при наличии стационарной динамической нагрузки и теоретически обосновать выбор оптимального периодического воздействия на поврежденные ткани с целью их скорейшего и устойчивого заживления.
В математическом виде задача о репаративной регенерации (репарации, или структурной перестройки) костной ткани была сформулирована относительно недавно. Предложена теория [1], согласно которой действие механических сил на процесс дифференциации клеток и структурной перестройки кости осуществляется посредством поля упругой деформации. Гипотеза состояла в том, что характер регенерации ткани определяется комбинацией первого и второго инварианта тензора деформаций. В частности, был рассмотрен процесс возникновения фиброзной, фиброзно-хрящевой и хрящевой ткани из недифференцированного вещества мезенхимальной ткани, возникающей после первоначального заживления в зоне перелома кости. Было выдвинуто предположение, что сдвиговые компоненты тензоров упругости и деформаций являются специфическим стимулом для развития коллагеновых волокон, т.е. фиброзной ткани, в то время как гидростатические составляющие ответственны за формирование хрящевой ткани.
В качестве регулирующих механических стимулов были введены положительное гидростатическое напряжение и знакопеременное октаэдрическое касательное напряжение, а их линейная комбинация трактовалась как "индекс остеогенеза" [2, 3]. Краевая задача линейной теории упругости изотропной однородной среды решалась в статической постановке итерационно с использованием метода конечных элементов. Гипотетически предполагалось, что большие значения индекса остеогенеза обуславливают преобразование хрящевой ткани в костную, а низкие значения индекса поддерживают развитие хряща, задерживая формирование костной ткани. Числовые значения параметров, входящих в индекс остеогенеза, определялись из сравнения анатомического строения типичного сустава скелета и соответствующего модельного решения. С помощью разработанного численного алгоритма в плоской постановке краевой задачи было рассчитано возникновение зон оссификации в процессе восстановления кости в области перелома и заживления кости вокруг ортопедического им-плантата [4, 5].
Подходы, представленные в публикациях Картера и Джиори с соавторами [2—5], дают реалистичное предсказание о зоне оссификации моделируемого органа. Однако используемая статическая постановка задачи не позволяет, начиная с некоторой на-
чальной конфигурации, исследовать развитие области оссификации. Это приводит к необходимости введения времени в уравнения и рассмотрения изменений механических свойств ткани в процессе решения. Кроме того, статическая постановка существенно сужает возможности изучения влияния временных характеристик прикладываемой нагрузки на структурную перестройку ткани.
Были предложены количественные выражения пороговых величин, входящих в закон развития кости из мезенхимальной ткани [6]. Согласно выдвинутой гипотезе интра-мембранная оссификация, которая имеет место в основном в плоских костях и на плоских поверхностях трубчатых костей, возникает, если локальные деформации меньше 5%. Эндохондральное (внутрихрящевое) развитие, обеспечивающее рост трубчатых костей в длину, возникает под действием сжимающего гидростатического давления выше 0.15 МПа и локальных деформациях меньше 15%. Другие виды напряженно-деформированного состояния органа стимулируют формирование фиброзной соединительной ткани или смешанной фиброзно-хрящевой. С помощью представленного управляющего правила и коммерческого конечно-элементного комплекса ANSYS была решена задача о моделировании процесса восстановления упругих свойств костной ткани в зоне перелома в осесимметричной постановке.
В предыдущих публикациях биологическая ткань рассматривалась как упругая однородная среда. Однако известно, что костное вещество — высоко структурированный биологический композитный материал, образованный матрицей из коллагена и кристаллов гидроксиапатита и сложной иерархической системой пор, заполненных сосудистой и внутритканевой жидкостью. Поскольку основная масса специфических костных клеток располагается на стенках канальцев и поверхностях лагун, входящих в систему микропор, то предполагается, что возмущения, вносимые внешней механической нагрузкой в установившееся движение жидкости в транспортной системе кости могут обеспечивать передачу управляющих сигналов между клетками кости в процессе ее структурной перестройки [7, 8]. Перечисленные ограничения были в значительной степени преодолены [9, 10] при рассмотрении задачи в динамической постановке в случае двухфазной изотропной модели среды. Альтернативно было предложено новое управляющее правило, определяющее фенотип ткани, образующейся в текущей точке среды при механической стимуляции. Вместо "индекса остео-генеза" был введен безразмерный "механо-регулирующий индекс" (МРИ) — линейная комбинация максимального значения октаэдрической сдвиговой деформации упругого каркаса двухфазной среды и максимального значения скорости потока внутритканевой жидкости в порах с двумя эмпирическими коэффициентами. С помощью коммерческого конечно-элементного комплекса MARC, имеющего возможности анализа насыщенных грунтов в статической пороупругой постановке, была решена осесимметричная задача о моделировании процесса формирования костной ткани вокруг имплантата под действием циклической нагрузки, приложенной к протезу. Согласно алгоритму итерационного приближения упругие модули и проницаемость каждого элемента исследуемого слоя ткани, прилегающего к поверхности имплантата, на каждой итерации изменялись в соответствии со значением МРИ до тех пор, пока решение не сходилось к некоторому устойчивому состоянию. В результате было получено конечное распределение вида ткани по основным фенотипам: фиброзная, хрящевая и костная, соответствующее, как было указано в статьях, ранее проведенным клиническим экспериментам.
Рассмотренная модель интенсивно применялась и в последующих работах [11 — 14]. Основное усовершенствование состояло в том, что в механобиологическую схему было включено уравнение диффузии, описывающее миграцию и распределение "клеток-предшественников", которые формируются на периостальной поверхности кости из внешних к костной мозоли мягких тканей и костного мозга. Для всех тканей была
принята статическая модель линейного пороупругого тела, твердая и жидкая фаза которого сохраняли свойство сжимаемости во время всего расчета. С помощью рассмотренной модели и различных коммерческих конечно-элементных комплексов (ABAQUS, ABAQUS Inc., USA и DIANA, TNO, The Netherlands) были рассмотрены задачи о трансформации живой ткани вокруг биосовместимых имплантатов, заживления перелома кости и дистракционного остеогенеза.
Разработанная концепция при незначительных модификациях применялась впоследствии и другими авторами для решения различных практических задач ортопедии и биоинженерии тканей [15, 16]. С теоретической точки зрения интересно также отметить статьи [17, 18], в которых с помощью фундаментальных законов и соотношений механики сплошных сред рассматривались обобщенные модели роста живой ткани.
К основному недостатку рассмотренных работ можно отнести то, что моделировалось действие только статической или циклической нагрузки низкой частоты (1 Гц). Несмотря на декларируемый динамический подход, фактически задача структурной перестройки решалась в квазистатической постановке, при использовании которой инерционными эффектами пренебрегалось. Однако экспериментально показано, что резонансные режимы относительно высокой частоты могут оказывать более существенный эффект для стимулирования восстановления механических свойств, чем низкочастотная периодическая нагрузка [19]. Поэтому в настоящей статье предложен вычислительный алгоритм регенерации кости, основанный на полной динамической модели анизотропной пороупругой среды с учетом всех инерционных слагаемых, характеризующих движение упругого каркаса и внутритканевой жидкости.
1. Общий алгоритм. Блок-схема алгоритма математического моделирования репара-тивной регенерации костной ткани представлена на фиг. 1. В начальный момент времени ткань в исследуемой области структурной перестройки относится к гранулированному типу, клетки которой недифференцированны. Данный вид ткани образуется при быстром заживлении ран большого размера и содержит большое количество мелких кровеносных сосудов, позволяющих активным клеткам-предшественникам мигрировать внутрь костной мозоли. Исходные значения пороупругих модулей модели костной мозоли соответствуют этому типу ткани.
Расчетная схема состоит из двух параллельных алгоритмов, основанных на единой конечно-элементной модели исследуемого биомеханического объекта, которым может быть часть костно-мышечной системы или фрагмент кости, включающие костную мозоль или внутрикостный имплантат. Первый алгоритм реализует вычислительную схему на основе уравнения диффузии и позволяет рассчитать концентрацию активных клеток-предшественников, способных к дифференциации и производству других типов клеток. В результате на каждом шаге по времени определяется пространственное распределение концентрации активных клеток, мигрирующих из тканей, окружающих область репарации кости.
Вторая часть схемы структурной перестройки реализует динамические уравнения пороупругого тела и предназначена для расчета основных пе
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.