CC BY
154
УДК 531/534:[57+61]
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ СТРУКТУРНОЙ ПЕРЕСТРОЙКИ КОСТНОЙ ТКАНИ
Л.Б. Маслов
Кафедра теоретической и прикладной механики Ивановского государственного энергетического университета имени В.И. Ленина, Россия, 153003, Иваново, ул. Рабфаковская, 34, e-mail:
Аннотация. Рассмотрены основные научные подходы к решению фундаментальной задачи биомеханики живых тканей, состоящей в разработке математической теории структурной перестройки твердых биологических тканей, детерминированной процессом дифференциации клеток и управляемой внешним силовым полем. Представлены обобщенная динамическая модель изменяющейся пороупругой сплошной среды и математический алгоритм, концептуально описывающий процесс структурной перестройки костной ткани под действием внешнего механического стимула периодического характера. Модель предполагает, что индуцируемые потоки жидкости в системе микропор костного вещества, возникающие в результате деформаций костного матрикса, наряду с самими упругими деформациями играют существенную роль механического регулятора, ответственного за репаративную регенерацию ткани. Для численного анализа построенной математической модели предложены численные алгоритмы, на основе которых может быть разработано программное обеспечение для компьютерного моделирования процесса костной перестройки. Математическая модель дает возможность исследовать процессы восстановления поврежденных костных элементов опорно-двигательного аппарата человека при наличии динамической нагрузки и теоретически обосновать выбор оптимального периодического воздействия на поврежденные ткани с целью их скорейшего и устойчивого заживления.
Ключевые слова: костная ткань, структурная перестройка, механический стимул, пороупругость, колебания, математическое моделирование.
Введение
Известно, что живые ткани в процессе своего роста и развития существенным образом реагируют на внешнее силовое поле, в котором они функционируют. Механический фактор оказывает стимулирующее и регулирующее воздействие на специфические клетки тканей, что приводит к запуску и развитию процессов трансформации органа в макроскопическом масштабе [29].
Развитие скелета и его функционирование в течение жизни обеспечивается и опосредуется процессами роста, формообразования, ремоделирования и репарации костной ткани. Первые два процесса генетически детерминированы и свойственны живой ткани при развитии до взрослого состояния, но при этом реализуются наряду с генетической программой под действием биомеханических и физико-химических факторов [5]. Физиологические алгоритмы репарации и ремоделирования действуют на
© Маслов Л.Б., 2013
Маслов Леонид Борисович, д.ф.-м.н., доцент, заведующий кафедрой теоретической и прикладной механики, Иваново
протяжении всей жизни организма и служат основой преобразований костной ткани в результате процесса сложной интеграции механических стимулов на тканевом и клеточном уровнях при действии внешнего силового поля.
Структурная перестройка неразвитой мягкой субстанции в плотную костную ткань происходит в результате дифференциации костных клеток, например, во время восстановления целостности кости после перелома и вживления скелетных имплантатов в твердое вещество костной ткани, что приводит к запуску процесса репаративной регенерации кости в зоне соприкосновения с поверхностью инородного предмета или между костными отломками.
Начиная с работы Вольфа 1892 года [30] к настоящему времени накоплен значительный опыт эмпирических исследований. Так, например, медицинские данные свидетельствуют, что недостаток механической стимуляции (микрогравитация или длительный период малоподвижности) приводит к потере костной массы в трабекулярном веществе кости, а также к увеличению интракортикальной порозности костей опорно-двигательного аппарата человека. Напротив, имеются данные, полученные при исследовании атлетов, занимающихся специфическими видами спорта, что увеличение физической нагрузки ведет к увеличению костной массы и толщины слоя компактного вещества кости. Существует гипотеза, имеющая клинические подтверждения, что после первоначальной дифференциации мягкого вещества костной мозоли для успешного завершения процесса остеосинтеза переломов длинных трубчатых костей необходима дозированная физическая нагрузка на конечность [22]. Дефицит же механической нагрузки приводит к торможению процесса оссификации желеобразной костной мозоли и образованию из нее преимущественно фиброзной или хрящевой ткани. Экспериментально показано, что решающим для инициации ростовых процессов является не пиковое значение деформации ткани, а определенное значение нижней границы деформации, при превышении которого процессы новообразования ткани превалируют над резорбцией [11]. Установлено также, что определенное значение имеют и другие характеристики механического воздействия, такие как распределение напряжений и потоков внутритканевой жидкости в порах костного вещества, длительность и частота силового воздействия. Высказывается предположение, что динамические, в частности периодические, нагрузки более эффективны для стимуляции репаративной регенерации кости, чем статические [5, 29].
Несмотря на наличие общей концепции морфогенеза и экспериментальных подтверждений адаптации живых тканей к изменению их функционального состояния, в настоящее время общепринятой универсальной математической теории, адекватно и всесторонне описывающей механизмы структурной перестройки живых тканей и механохимических регуляций на системном и локальном тканевом уровнях, не разработано.
В статье представлен краткий обзор подходов к математическому моделированию процесса структурной перестройки костной ткани на основе дифференциации клеток и действия силовой нагрузки как важного регулирующего фактора, дано описание математической модели и алгоритма, учитывающего динамические эффекты относительно высокой частоты.
Обзор теоретических подходов
Задача о функциональной адаптации кости к изменению функциональной нагрузки фактически восходит к пионерской работе Вольфа [30] и является наиболее изученной. Существует достаточно большое количество публикаций, посвященных экспериментальным данным и моделям адаптирующейся ткани, приведены результаты компьютерного анализа, качественно верно описывающие процессы адаптации губчатой или кортикальной ткани к изменившимся силовым условиям [1, 26].
Однако задача дифференциации клеток и структурной перестройки костной ткани в математическом виде была сформулирована относительно недавно. В работе Пауэлса 1941 года [25] была предложена теория, согласно которой влияние механических сил на процесс дифференциации осуществляется посредством поля упругой деформации ткани. Математически гипотеза состояла в том, что комбинация первого и второго инвариантов тензора деформаций определяет характер структурной перестройки ткани. В частности, был рассмотрен процесс возникновения фиброзной, фиброзно-хрящевой и хрящевой ткани из недифференцированного вещества мезенхимальной ткани, возникающей после первоначального заживления в зоне перелома кости.
Современное изучение влияния внешнего механического стимула на структурную перестройку биологических тканей, обусловленную клеточной дифференциацией, методами математического моделирования и механики сплошных сред начинается с уже упоминавшейся работы Пауэлса [25] и его последующих публикаций. Впервые была высказана гипотеза, что сдвиговые компоненты тензоров упругости и деформаций являются специфическим стимулом для развития коллагеновых волокон, т.е. фиброзной ткани, в то время как гидростатические составляющие ответственны за формирование хрящевой ткани.
Впоследствии концепция Пауэлса была применена при разработке различных биомеханических моделей дифференциации тканей на основе линейной теории упругости изотропной однородной среды. Краевая задача в статической постановке решалась итерационно с использованием метода конечных элементов. В работах [12, 14] в качестве регулирующих механических стимулов были введены положительное гидростатическое напряжение и знакопеременное октаэдрическое касательное напряжение для каждого вида нагрузки. Их комбинация трактовалась как «индекс остеогенеза»:
где п - количество циклов нагружения для каждого вида нагрузки; а -гидростатическое напряжение; к - эмпирическая константа; т - октаэдрическое касательное напряжение.
Гипотетически предполагалось, что большие значения 01 обусловливают преобразование хрящевой ткани в костную - так называемое эндохондральное формирование кости, а низкие значения индекса поддерживают развитие хряща, задерживая формирование костной ткани. Из сравнения модельного решения и анатомического строения типичного сустава опорно-двигательного аппарата были сделаны выводы, что наилучшее соответствие дают значения к = 0,3...1. Однако конкретных пороговых значений индекса остеогенеза представлено не было. С помощью разработанного численного алгоритма в плоской постановке краевой задачи было рассчитано возникновение зон оссификации в процессе восстановления кости в области перелома и заживления кости вокруг ортопедического имплантата
Подходы, представленные в публикациях [12-14], дают реалистичное положение зоны оссификации моделируемого органа. Однако используемая статическая постановка задачи не позволяет исследовать развитие области оссификации, начиная с некоторой начальной конфигурации. Это приводит к необходимости введения времени в уравнения и рассмотрения изменений материальных свойств ткани в процессе решения. Кроме того, статическая постановка существенно сужает возможности изучения влияния временных характеристик прикладываемой нагрузки на структурную перестройку ткани.
(1)
[13, 17].
Авторы работы [15] предложили количественные выражения пороговых величин, входящих в закон развития кости из мезенхимальной ткани. Согласно их гипотезе интрамембранная оссификация, которая имеет место в основном в плоских костях и на плоских поверхностях трубчатых костей, возникает, если локальные деформации меньше 5%. Эндохондральное (внутрихрящевое) развитие, обеспечивающее рост трубчатых костей в длину, возникает под действием сжимающего гидростатического давления выше 0,15 МПа и локальных деформаций меньше 15%. Другие виды напряженно-деформированного состояния органа стимулируют формирование фиброзной соединительной ткани или смешанной фиброзно-хрящевой. С помощью представленного управляющего правила и коммерческого конечно-элементного комплекса ANSYS (ANSYSInc., USA) была решена задача о моделировании процесса восстановления костной ткани в зоне перелома в чисто упругой осесимметричной постановке.
Предыдущие публикации рассматривали биологическую ткань как классическую сплошную среду. Однако известно, что костное вещество является высоко структурированным биологическим композитным материалом, образованным матрицей из коллагена и кристаллов гидроксиапатита, и сложной иерархической системой пор, заполненных сосудистой и внутритканевой жидкостью. Поскольку основная масса специфических костных клеток располагается на стенках канальцев и поверхностях лагун, входящих в систему микропор, то предполагается, что возмущения, вносимые внешней механической нагрузкой в установившееся движение жидкости в транспортной системе кости, могут обеспечивать передачу управляющих сигналов между клетками кости в процессе ее структурной перестройки [6]. Перечисленные ограничения были в значительной степени преодолены в работах [19, 28], в которых задача рассматривалась в динамической постановке и была принята двухфазная, хотя и изотропная, модель среды. Альтернативно было предложено новое управляющее правило для описания процесса перестройки неразвитой соединительной ткани в плотную хрящевую или костную ткань. Для этой цели был введен безразмерный механорегулирующий индекс, который определяет, ткань какого фенотипа образуется в текущей точке среды в ответ на механическую стимуляцию:
m=-+q, (2)
a b
где s - максимальное значение октаэдрической сдвиговой деформации упругого каркаса двухфазной среды; q - максимальное значение скорости потока внутритканевой жидкости в порах; а, b - эмпирические константы, а = 0,0375 и b = 3 мкм/с.
На основе имеющихся экспериментальных данных авторы [19] использовали следующее управляющее правило с двумя эмпирическими пороговыми значениями механорегулирующего индекса (2): 1) М > 3 означает формирование фиброзной ткани с модулем упругости Е = 2 МПа и проницаемостью k = 0,01 мм4/Нс; 2) 1 < М < 3 -формируется хрящевая ткань (Е = 10 МПа, k = 0,005 мм4/Нс); 3) М < 1 - формируется костная ткань (Е = 4590 МПа, k = 0,37 мм4/Нс). С помощью коммерческого конечноэлементного комплекса MARC (Palo Alto, USA), имеющего возможности анализа грунтов в статической пороупругой постановке, была решена осесимметричная задача о моделировании процесса формирования костной ткани вокруг имплантата под действием циклической нагрузки, приложенной к протезу [19, 28]. Согласно алгоритму итерационного приближения упругие модули и проницаемость каждого элемента исследуемого слоя ткани, прилегающего к поверхности имплантата, на каждой итерации изменялись в соответствии со значением механорегулирующего индекса до тех пор, пока решение не сходилось к некоторому устойчивому состоянию.
В результате было получено конечное распределение вида ткани по основным фенотипам: фиброзная, хрящевая и костная, соответствующее, как было указано в статье, ранее проведенным клиническим экспериментам.
Рассмотренная модель интенсивно применялась и в последующих работах отмеченных исследовательских групп [20, 21, 23, 27]. Основное усовершенствование состояло в том, что в механорегулирующую схему было включено диффузионное уравнение, описывающее миграцию и распределение так называемых клеток-предшественников, которые формируются на периостальной поверхности кости из внешних к костной мозоли мягких тканей и костного мозга:
dV = JAv, (3)
dt
где у - объемная концентрация клеток; J - коэффициент диффузии, м2/день.
Уравнение (3) решалось независимо от остальных дифференциальных уравнений сплошной среды. Предполагалось, что клетки в каждом конечном элементе могут дифференцироваться в клетки основных типов тканей: фибробласты,
хондроциты и остеобласты в зависимости от среднего механического состояния элемента в текущий расчетный день. Каждый тип клеток производит твердую фазу (упругую матрицу) двухфазного материала с определенными модулями упругости. Конечный модуль упругости сплошной среды рассчитывался на основе полученных частных значений на каждой итерации. Резорбция кости моделировалась путем простой деактивации соответствующего элемента. Костный мозг и оригинальная кортикальная ткань были исключены из зоны структурной перестройки. Для всех тканей была принята статическая модель линейного пороупругого тела, твердая и жидкая фаза которого сохраняли свойство сжимаемости во время всего расчета. С помощью рассмотренной модели и различных коммерческих конечно-элементных комплексов (ABAQUS, ABAQUSInc., USA и DIANA, TNO, The Netherlands) были рассмотрены задачи трансформации живой ткани вокруг биосовместимых имплантатов, заживления перелома кости и дистракционного остеогенеза.
Разработанная в статьях [19, 23, 28] концепция при незначительных
модификациях применялась впоследствии и другими авторами для решения различных практических задач ортопедии и биоинженерии тканей [10, 24]. Интересно отметить статьи [7, 8], в которых с помощью фундаментальных законов и соотношений механики сплошных сред рассматривались обобщенные модели роста живой ткани.
К основному недостатку перечисленных работ можно отнести то, что моделировалось действие только статической или циклической нагрузки низкой частоты (1 Гц). Поэтому, несмотря на декларируемый динамический подход, фактически задача структурной перестройки решалась в квазистатической постановке, при использовании которой инерционные эффекты не учитывались. Однако экспериментально показано, что резонансные режимы относительно высокой частоты могут оказывать более существенный эффект при стимулировании живых тканей, чем низкочастотная периодическая нагрузка [18]. Поэтому далее в статье будет рассмотрена полная динамическая модель анизотропной пороупругой среды с учетом всех инерционных слагаемых, характеризующих движение упругого каркаса и внутритканевой жидкости.
Модель пороупругой среды
Рассмотрим основные уравнения механики пористого материала, насыщенного жидкостью или газом. Пороупругая сплошная среда является моделью гетерогенного материала, одна из фаз которого представляет собой упругую пористую матрицу, или
скелетон, а вторая - жидкость, заполняющую систему пор. Определяющие соотношения и дифференциальные уравнения, описывающие малые упругие перемещения эффективной двухфазной среды, были сформулированы Био на основе феноменологического подхода и впоследствии получили обобщение в работах Нигматулина [4, 9].
Для описания деформирования основных фенотипов биологических тканей, таких как компактное и губчатое вещество кости, хрящевая и соединительная ткани, используются определяющие соотношения двухфазного материала, записанные относительно осредненных по представительному элементу среды перемещений твердой и и жидкой и фаз:
о
(1)
(2)
(и, и ) = С ••£ ( и ) + О©( и ), (и, и ) = О ••£ (и ) + Я©(и ),
(4)
где о(1), 5(2) - осредненные по всему объему представительного элемента тензор и шаровая часть тензора напряжений в материале каждой фазы, 5(2) = Е • о(2)/3 ; £, © -тензор деформаций и объемная деформация твердой и жидкой фаз; С - тензор упругих модулей твердой фазы; О - тензор коэффициентов взаимности, определяющих влияние деформаций одной из фаз на возникающие напряжения в другой фазе; Я -гидростатическая константа, имеющая смысл модуля объемного сжатия жидкой фазы; Е - единичный тензор.
В предположении упругой модели твердой фазы и модели идеальной сжимаемой жидкости фазовые уравнения (4) преобразуются в определяющие соотношения пороупругой среды в «и-р» переменных (перемещение упругого каркаса - давление поровой жидкости). Кинематической переменной жидкой фазы
является вектор относительного перемещения w = Ф(и—и) или дивергенция этого вектора ^ = ~^^, имеющая физический смысл относительного изменения объемного содержания жидкости в порах. Для случая анизотропии упругих и гидростатических свойств определяющие соотношения примут следующий вид:
о (и Р) = °*г(и) — АР = С*г - £ (и) — АР,
С(и, р) = а • £ (и) + ф2Я—1 р,
где о - полный тензор напряжений; о^ - тензор напряжений в точках твердой фазы,
вызываемый только упругими деформациями; Саг - тензор упругих модулей среды
в дренированном состоянии; А - тензор коэффициентов эффективных напряжений Био; р - давление жидкости в поровых каналах.
В результате перехода к новым переменным возникают эффективные материальные константы, представленные тензором упругих модулей в дренированном состоянии и тензором коэффициентов Био:
С* = С — ООЯ л
А = ф(Е + ОЯ~:). ( )
Уравнения динамики эффективной среды могут быть получены из уравнений баланса количества движения, записанных для каждой фазы, с учетом силового и инерционного взаимодействия между фазами:
У-°(1) +(1— ф) Г, + Я = р„и + р12 и, (7)
V •о(2) +ф^ — Я = Р121! + Р22'Ц^,
где Г,, - плотности объемных сил, действующих в твердой и жидкой фазах;
Я - сила межфазного взаимодействия; р11, р12, р22 - частичные фазовые плотности, выражаемые через истинные плотности упругой матрицы р, и заполняющей поры жидкости р/, пористость ф и параметр искривленности поровых каналов т,
р11 = ( -ф) —Pl2, р12 =(1 —т)фр/ , р22 =тфр/ .
Сила Я для схемы Рахматулина силового взаимодействия и совместного деформирования фаз складывается из равновесной и диссипативной составляющих [4]:
К = К0 + К= - РVф+B • , (8)
где В - тензор коэффициентов вязкого трения, выражаемый через тензор
проницаемости и вязкость жидкости: В = ф2К—1, К = к/п/ .
С учетом выражения тензора напряжений в идеальной сжимаемой жидкости и соотношения (8) уравнения (7) после сложения и преобразования примут вид
^о (и, р ) + Гу =ри + р/ ,
1 1 (9)
—Vp + Г/ =р/и + тф р/-(у + фК • ,
где Гу - полная объемная сила; р - полная плотность среды.
Динамические уравнения (9) содержат дополнительную переменную w, которая выражается с помощью комплексного преобразования Лапласа:
^ = -К (, )-^р + р/, 2 и -/), (10)
где К (,) - приведенная комплексная гидравлическая проницаемость пористой среды,
К (, ) = ( 5 Е + тф-1р ^ 2К )-1 • ,К.
Окончательная система уравнений примет вид краевой задачи относительно изображений искомых функций и и р [3]:
-2 (Е-р/Г(5))-(А-Г(5))^р = -у + Г(5)•?/,
(11)
, -1 V• (К (, ) ■VI)) - ф2 Я-р -(А - Г (, )) • £ = ,_1к /,
где ]Г - тензор, характеризующий проницаемость среды и взаимодействие фаз, Г = р/^К; к/ - плотность внутренних источников жидкости, к/ =V • ( К • Гf ) .
Полученная система уравнений (11) с учетом определяющего уравнения (5) представляет собой замкнутую связанную систему уравнений динамики относительно изображений вектора перемещений упругого формообразующего скелета и давления жидкости в порах. При замене параметра Лапласа на комплексное выражение 5 = /ш получим уравнения вынужденных колебаний пороупругой среды под действием силы, изменяющейся по гармоническому закону. Не умаляя общности, можно принять естественное в постановке задачи об исследованиях вынужденных колебаний пористой упругой конструкции условие отсутствия объемных сил, действующих на точки жидкой фазы (Г/ = 0), что приводит к упрощению правой части уравнений (11).
Записанная система дифференциальных уравнений и сформулированная на их основе краевая задача пороупругости, описывающая установившиеся гармонические колебания пороупругого анизотропного тела, учитывают все основные виды межфазного взаимодействия (силовой, физический и инерционный), что отличает их от других известных в литературе соотношений.
В случае смешанной «и-р» формулировки задачи пороупругости (11) дискретная система содержит изображения искомых функций и параметр преобразования Лапласа. Если нагрузки являются произвольными функциями от времени, то для нахождения динамического отклика системы на силовое воздействие требуется использовать обратное преобразование Лапласа и специальные численные алгоритмы. Однако если внешние нагрузки изменяются по гармоническому закону, то отклик линейной системы также представляет собой гармонические во времени функции, а параметр преобразования Лапласа 5 = /ш , что дает следующую систему матричных уравнений [3]:
(К( - ш2М - Ь (/ш)) и -(Н1 + Й2 (/ш)) Р = Жу + ^,
Г ч (12)
-(й1 + Н2 (/ш)) и + (-О + /ш-1С (/ш)) Р = -/ш-1О*,
где КЛг, М - стандартные глобальные матрицы жесткости и массы; Ь -дополнительная глобальная матрица массы; Н1, Н 2 - глобальные матрицы взаимного
влияния распределения давления жидкости в порах и деформации скелета; О, С -глобальные матрицы насыщения и проницаемости; и, Р - глобальные векторы комплексных амплитудных значений перемещений и давлений в узлах конечноэлементной сетки.
Дифференциации и диффузия клеток
Биологической компонентой математической модели структурной перестройки костной ткани являются процесс проникновения активных, способных к видоизменению клеток в зону формообразования и процесс преобразования клеток в другие фенотипы. Характерным примером перестройки ткани в результате дифференциации клеток и воздействия механического фактора служит явление оссификации костной мозоли, возникающей в зоне перелома длинных трубчатых костей скелета.
Известно, что перелом разрушает местное кровоснабжение, вызывая кровотечение, кислородное голодание, гибель клеток и асептическое воспаление. Образующаяся гематома может быть первоисточником сигнальных молекул, которые обладают способностью запускать внутриклеточные процессы дифференциации и формирования костных клеток. На следующем этапе недифференцированные клетки начинают производить новые клетки, которые преобразуются в сосуды, фибробласты, межклеточный материал и поддерживающие клетки. Все вместе они формируют гранулированную ткань.
Дальнейшая пролиферация, дифференциация и организация клеток приводит к образованию в гранулированной ткани новых хондроцитов и остеобластов в течение так называемого мезенхимного процесса [5, 29]. Предполагается, что клетки-предшественники (мезенхимные стволовые клетки) могут проникать в области дифференцирования кости, хряща, сухожилия, связки, мышцы и соединительной ткани. Выбрав свой путь, клетки начинают синтезировать внутриклеточный органический матрикс ткани, который впоследствии минерализуется. В процессе сращения перелома минерализация длится несколько недель и завершается образованием переломной мозоли.
Внешняя костная мозоль может быть разделена на две компоненты: твердую мозоль с внутримембранным окостенением и мягкую мозоль, в которой продолжается внутрихрящевое окостенение. Самые важные стадии процесса сращения происходят в надкостнице в первые несколько дней сращения. Соединительно-мембранная кость формируется из недифференцированных клеток надкостницы. В течение одной-двух недель образуется множество удлиненных разросшихся хондроцитов, которые претерпевают митоз и деление. Далее количественный рост клеток уменьшается и гипертрофированные хондроциты становятся клетками главного типа в костной мозоли. После кальцинирования хряща в нем разрастается сеть кровеносных сосудов и он заполняется остеобластами. Данная стадия завершается образованием незрелой костной ткани внутри мозоли, состоящей из молодых клеток здорового типа.
После первичной оссификации костной мозоли при условии оптимальной механической нагрузки начинаются явления реконструкции кости в зоне перелома. Первично минерализованный хрящ полностью превращается в незрелую костную ткань путем дополнительной минерализации. С течением времени происходит замена части незрелой костной ткани пластинчатой костной тканью, мозоль между костными отломками замещается новыми вторичными остистыми отростками пластинчатой ткани в направлении приложенной нагрузки, полость мозоли заполняется костным мозгом, формируется медуллярная полость. На этом основной этап сращения кости считается завершенным. Однако еще в течение нескольких лет процесс реконструкции продолжается: участок большой грубой кости со шрамом уменьшается и окружающая костная мозоль постепенно исчезает.
Для математической формализации описанных процессов существующие модели структурной перестройки костной ткани предполагают, что доставка активных клеток-предшественников, способных к формоизменению, осуществляется непрерывно в течение заданного срока реконструкции по закону диффузии. Источником недифференцированных клеток являются окружающие кость мягкие ткани, костный мозг, надкостница. Уравнение диффузии представляет собой классическое уравнение Фурье, имеющее вид (3) в изотропном случае при постоянном значении коэффициента диффузии. Применительно к анизотропным средам и, как и ранее, при отсутствии объемных источников оно может быть записано следующим образом:
где I - тензор коэффициентов диффузии.
С точки зрения компьютерной реализации важно отметить, что дифференциальное уравнение в частных производных (13) аналогично второму уравнению системы (11) без перекрестного слагаемого, если перейти от изображений функций к оригиналам. Это позволяет существенно упростить разработку численного алгоритма решения комплекса уравнений при использовании одного типа конечных элементов и, следовательно, единой конечно-элементной сетки для описания поля перемещений и распределения концентрации клеток в объеме тела. Применяя метод взвешенных невязок к уравнению (13) и процедуру конечно-элементной дискретизации, получим систему дифференциальных уравнений первого порядка в матричном виде:
где О, С - глобальные матрицы системы конечно-элементных уравнений, аналогичные матрицам насыщения Б и проницаемости С в системе (12); Т, Т - глобальные векторы узловых значений скорости изменений концентрации
(13)
Б ¥ + О ¥ = О,
(14)
и концентрации клеток; О - глобальный узловой вектор заданных потоков вещества (клеток) через поверхность рассматриваемого объема тела.
Систему обыкновенных дифференциальных уравнений (14) можно проинтегрировать численно с помощью разностной схемы первого порядка (методом Эйлера), что дает матричное уравнение для определения глобального узлового вектора концентрации клеток в любой момент времени:
БТ(к+1) =(б-МС)Т(к) +МО, М = /(к+1) -/(к), к = 1,2, ... . (15)
Нахождение вектора Т из уравнения (15) на каждом шаге интегрирования может быть осуществлено любым известным способом решения больших систем алгебраических уравнений с разреженной матрицей.
Алгоритм структурной перестройки
Общий алгоритм математического моделирования структурной перестройки костной ткани представлен на рис. 1. В начальный момент времени ткань
в исследуемой области структурной перестройки относится к гранулированному типу, клетки ее недифференцированны. Данный вид ткани образуется при быстром заживлении ран большого размера и содержит большое количество мелких кровеносных сосудов, позволяющих активным клеткам-предшественникам мигрировать внутрь костной мозоли. Исходные значения пороупругих модулей конечно-элементной модели исследуемой области кости соответствуют ткани гранулированного типа.
Расчетная схема состоит из двух параллельных алгоритмов, основанных на единой конечно-элементной модели исследуемого биомеханического объекта (целостной костно-мышечной системы или фрагмента кости, включающей костную мозоль или внутрикостный имплантат). Первый алгоритм реализует вычислительную схему (15) на основе уравнения диффузии (13) и позволяет рассчитать концентрацию активных клеток-предшественников, способных к дифференциации и производству других типов клеток. В результате на каждом шаге по времени определяется пространственное распределение концентрации активных клеток, мигрирующих из тканей, окружающих область репарации кости.
Вторая часть схемы структурной перестройки реализует уравнения (12) и предназначена для расчета основных переменных задачи - перемещений и деформаций упругого скелета ткани (твердой фазы пороупругой среды), давления и потоков внутритканевой жидкости. Силы определенной частоты прикладываются пошагово по заданному временному закону, что позволяет исследовать влияние различных схем нагружения на механизм структурной перестройки кости.
После расчета ключевых переменных задачи в каждой точке пороупругой среды по критерию (2) рассчитывается значение механобиологического индекса М, определяющего, образование какого фенотипа ткани стимулирует данная нагрузка. В математическую модель структурной перестройки, кроме исходной гранулированной, включены четыре вида дифференцированной ткани - фиброзная, хрящевая, незрелая костная и зрелая костная, каждый из которых характеризуется своими значениями материальных констант (плотности, модули упругости, коэффициенты проницаемости и т.д.). В зависимости от значения индекса М в каждом узле конечно-элементной модели или во всем конечном элементе пороупругие модули сплошной среды изменяются.
Рис. 1. Блок-схема алгоритма структурной перестройки твердых биологических тканей
Поскольку значения материальных констант четырех основных типов тканей сильно различаются, а резкое изменение механических свойств регенерирующей ткани не соответствует реальному процессу заживления кости, то для повышения точности
алгоритма введен блок сглаживания пороупругих модулей по нескольким
последовательным временным шагам т, как правило, пяти-десяти. Для каждой последующей итерации после первых т шагов для любой константы пороупругой среды применяется следующая формула:
1 п+т-1 Г
Еп) = - I Е'к), п = 1,2, ..., пш,х, - (16)
т п М
где Е(п) - осредненные значения / -й материальной константы по т последовательным шагам; Е/к) - расчетные значения /-й материальной константы на временном шаге к, полученные из конечно-элементного анализа; Т - конечное значение времени работы алгоритма структурной перестройки ткани; А- - значение шага по времени.
Концентрация клеток в исследуемой зоне кости меняется на каждом шаге, но физико-механические свойства материала меняются не сразу, поскольку область, занимаемая конечным элементом, может состоять частично из дифференцированной и гранулированной тканей. Поэтому для сглаживания двойственной природы материала и усреднения свойств элемента применяют правило смешения. Новые значения материальных констант каждого конечного элемента, описывающего область репарации кости, рассчитываются на основе полученных осредненных значений и величины концентрации активных клеток в любой точке среды по закону смесей:
£;» = _^!_Щп) + фтах -У(и) E(gran), п = i 2, ., nmax, (17)
ф max ф max
где E;п) - окончательные значения упругих и гидростатических модулей
перестраивающейся ткани на n-м шаге в любой точке области репарации (для любого конечного элемента или узла конечно-элементной сетки); ф( п) - концентрация клеток в данной точке среды на n-м шаге; ymax - максимально достижимая концентрация клеток в данной точке среды, которую можно считать равной концентрации клеток здоровой кости; E(isran) - значения пороупругих модулей гранулированной ткани
в данной точке среды.
Из соотношения (17) следует, что на окончательный результат влияют не абсолютные значения концентрации клеток, а относительные, что позволяет не вводить в расчет ymax и записать следующее выражение для материальных констант перестраивающейся среды:
E;п) =M/(n)E;n)+(l-y(n) ))гап), М/(П) =^, п = 1, 2, ..., nmax, (18)
ф max
где ф(п) - относительные значения концентрации клеток в произвольной точке среды на п-м шаге по времени.
Таким образом, на последнем этапе алгоритма на каждом временном шаге происходит окончательное вычисление материальных констант регенерирующей ткани исходя из значений концентрации активных клеток ф(п), рассчитанных из соотношения (15), и значений соответствующих констант материала, полученных из конечноэлементного анализа напряженно-деформированного состояния рассматриваемой области ткани. Затем физико-механические свойства компьютерной модели обновляются, и численная схема повторяется для следующего временного шага. По результатам расчета на каждом временном шаге можно наблюдать за сходимостью процесса и стремлением изначально недифференцированной однородной ткани к преобразованию в тот или иной фенотип.
Одномерная модель структурной перестройки ткани
Для понимания приведенного механизма структурной перестройки костной ткани рассмотрим простейшую одномерную модель, описываемую уравнениями, следующими из общих соотношений. Будем считать, что изучаемая область репарации представляет собой небольшой участок вдоль оси x длиной 21, занимаемый сплошной пороупругой средой с постоянными коэффициентами. В этом случае динамические уравнения растяжения вдоль продольной оси x легко получить, если ввести допущения, что вектор перемещений, вектор внешних сил и тензор напряжений имеют ненулевые компоненты только вдоль оси x и зависят только от продольной координаты. Аналогичное требование зависимости только от продольной координаты
накладывается также на давление жидкой фазы и плотность объемных источников. Примем естественное в рассматриваемой постановке задачи условие отсутствия
специфических объемных сил, действующих на точки жидкой фазы (?/ = 0), и,
следовательно, отсутствия внутренних источников жидкости в порах ( к / = 0 ).
В результате преобразования системы (11) получаем два одномерных уравнения [3]:
- Ей" + (р - у ( ^ ) р г )2 и + (а - у ( ^ )) р ' = }Ух,
- К (я) р" + яф2Я-1 р + (а - у (я)) яй' = 0,
где й - изображение продольного перемещения точек стержня, й = йх (х, я); Е -модуль Юнга пороупругого материала в дренированном состоянии, Е = Е(хс1г); а -коэффициент эффективных напряжений Био (диагональная компонента тензора А вдоль оси х ), а = ах; у - коэффициент, характеризующий проницаемость среды и
взаимодействие фаз (диагональная компонента тензора 1" вдоль оси х ), у = ух (я); К -приведенная комплексная гидравлическая проницаемость пористой среды (диагональная компонента тензора К вдоль оси х ), К = Кх (я).
Выражения коэффициентов К и у следуют из общих выражений тензоров КС и 1" и имеют вид
К (я) =------Г—7Ф, ^ (я) = “/—Я ь /ф • (20)
П + яЬр т/ф П/Р/ + яь ТФ
Представляя переменные и правые части уравнений (19) в гармоническом виде и заменяя параметр Лапласа на выражение я = /ш, получим уравнения продольных вынужденных колебаний в частотной области:
- Ей"-(р-у (ш)р / )ш2 й + (а-у (ш))р ' = }гх,
I ^ /
- К (ш) р" + /шф2 Я-1 р + / (а-у (ш))ш й' = 0.
В системе (21) перемещение и давление, помеченные сверху крышечкой, представляют собой комплексные амплитуды соответствующих временно-зависимых
функций, изменяющихся по гармоническому закону. Коэффициенты К и у равны
соответственно
К (ш) =------Ь-----—, у (ш) = —-——----------—. (22)
П + /шЬр/ т/ф п Р/ +/шк т/ф
Система (21) допускает простые аналитические решения в виде разложения в ряд Фурье при определенных типах граничных условий, например:
й (0) = 0, р ' (0 ) = 0,
(23)
й ’ (I ) = 0, р (I ) = 0.
В случае моделирования структурной перестройки ткани, занимающей небольшую область длиной 21 вдоль продольной оси кости, будем рассматривать симметричную постановку задачи. Тогда физический смысл выражений (23) состоит в том, что в среднем сечении х = 0 продольное перемещение и поток жидкости в порах
равны нулю, а правое сечение х = I свободно от нагрузок и полностью проницаемо для поровой жидкости.
В этом случае, предполагая отсутствие внутренних источников, представим правую часть и решение в виде разложения в ряд по собственным функциям в соответствии с граничными условиями (23):
и(х) = е ип біп ^, р(х) = Е Рп
п=1,3,5... 21 п=1,3,5...
I
ппх ,
біп--------ах,
21
ппх
(24)
СОБ-
21 ’
где /п - коэффициенты разложения нагрузки в ряд Фурье.
Подставляя выражения (24) в уравнения (21) и принимая во внимание условие ортогональности тригонометрических функций, получим систему, состоящую из двух алгебраических уравнений, для расчета коэффициентов разложения ип и рп:
Е 02т) ^р^
ш
ип “(а“У(ш))—Рп = Ґп
т (а - у (ю))—ип
( /_._\2 ) 2 п-1
(25)
Рп = о.
Решение системы (25) может быть найдено любым известным алгебраическим методом. Например, пользуясь известным правилом Крамера решения систем линейных алгебраических уравнений, получаем следующие выражения для искомых переменных:
ип =
Рп =
(к (ш ) (пи/ 21 )2 + ішф2 Я -) /п А(ш)
іш (а-у (ш))(П/ 21 ))п
(26)
А(ю)
где знаменатель представляет собой определитель матрицы системы уравнений (25), являющийся характеристическим многочленом третьей степени относительно частоты ю с комплексными коэффициентами, также зависящими от ш :
Е(пп) -(-*>/)
(
Р-У Р/) ш
пп
Л
К\ ^ | + ішф2Я_1
іш(а-г)2 ^І = 0.
(27)
В результате подстановки (26) в (24) получаем комплексные частотные функции, характеризующие перемещения и давление поровой жидкости пороупругой среды, совершающей колебания под действием продольной гармонической силы:
<( х) = Е
(К (ш)(итс/ 21 )2 + ішф2 Я 1) /п ---------------^------------ —біп-
ппх
р( х) = Е
А(ш) 21
іш (а-у (ш))(п/ 21 ) ппх
(28)
А(ш)
-СОБ-
21
«=1,3,5
«=1,3,5
п=1,3,5
Для применения механорегулирующего критерия (2), определяющего, ткань какого фенотипа образуется в текущей точке среды, нужны эквивалентная сдвиговая деформация упругого скелета и скорость потока внутритканевой жидкости в порах. Поскольку рассматривается одномерная модель костной ткани, в которой сдвиг не учитывается, то в качестве деформации взята продольная составляющая (в = dw/ dx). Поток жидкости в пороупругой модели характеризуется переменной q = w, выражаемой через градиент давления по формуле (10). Следовательно, из выражений перемещения и давления (28) могут быть получены требуемые переменные:
Л ((®)(W21 )2 + )/21)) nnx
в( x) = > -------------——---- -----------cos-,
nil. А(ш) 21 (29)
V 4 ir( \ ^ Ца-У(®))(nV21 )2 fn . nnx q(x) = К (о) ^ — ----- ---- ---sin-
к=1,3,5..
A(o) 2l
Полученные формулы составляют основную расчетную часть алгоритма структурной перестройки ткани и определяют образование определенного вида ткани согласно механорегулирующему индексу (2), стимулируемого приложенной нагрузкой.
Образование клеток и их миграция в процессе сращения перелома еще недостаточно изучены. Этот процесс принимается случайным и ненаправленным, что соответствует диффузии частиц газа или жидкости, описываемой дифференциальным уравнением (13). В рассматриваемом одномерном случае соотношение (13) для расчета концентрации клеток примет вид
Af J VL.5v= 0, (30)
dxf dx) dt
где J - коэффициент диффузии (диагональная компонента тензора J вдоль оси x ), J = Jx .
Пусть в крайних сечениях зоны репарации x = ±l значение концентрации активных клеток-предшественников, мигрирующих внутрь области, поддерживается постоянным и равным ymax. Учитывая симметричную постановку задачи,
подразумевающую отсутствие потока вещества в центральном сечении, краевые условия могут быть записаны следующим образом:
У( °,t )= °, у(l, t)=^max. (31)
При условии, что коэффициент диффузии не зависит от координаты, уравнение (29) имеет аналитическое решение в виде суммы ymax и разложения по методу Фурье:
ад
^(x t) = Vmax +Z akWk (x)eXP (-J^kt), (32)
k=1
где ak - коэффициенты общего решения однородной задачи; wk, Xk - собственные
функции и собственные значения соответствующей однородной краевой задачи.
Уравнение (33) при однородных граничных условиях, следующих из соотношений (31), имеет собственные решения:
* к =
^(2к- 1)тсУ , ч (к- 1)nx
-------— , wk (x) = cos----------------------------------- —, к = 1, 2, ... . (33)
2l ) кУ ’ 2l
Для вычисления неопределенных коэффициентов ак используем начальное
условие, состоящее в том, что в исходный момент времени в области репарации отсутствуют активные клетки-предшественники. Тогда из (32) получаем
ад
-утах =Е а^к (х). (34)
к=1
Выражение (34) можно рассматривать как разложение функции утах в ряд по собственным функциям, имеющий смысл ряда Фурье для данного вида функций ~№к (33). Исходя из ортогональности собственных функций, могут быть найдены коэффициенты ак:
1 I I
ак = Й N [Х^тах К (Х), где |Ы = { ^к (Х)2 ^. (35)
\гк\\ о о
После вычисления коэффициентов (35) окончательное выражение распределения концентрации активных клеток примет вид
, ч 4 ■^(-1)Аг1 (к -1) их ( 3 (2к -1)2 п2/ ^
*(*. о=*_-у- 1^-1 со$-^г~ехр —4Р—
(36)
Соотношение (36) позволяет рассчитать концентрацию активных клеток в зоне репарации на каждом временном шаге, что необходимо для применения правила смешения (17) или (18), определяющего значения материальных констант ткани после завершения текущей итерации. Поскольку на окончательный результат влияют не абсолютные значения концентрации клеток, а относительные, то (36) удобно переписать в следующем виде:
\к-1 ( т1гм_ л\2 _2.Л
(37)
* (х,, )==, - 4 (-1^ С05 (2к - О** еХр (-3 (2к - 2) П <
К ’ * ..х п к=1 2к -1 2/ ^ 4/
Следовательно, в рассматриваемой схеме движения активных клеток-предшественников остается только один независимый параметр J. Его значение выбирается таким образом, чтобы за расчетный период, принятый равным 4 месяцам или приблизительно 120 дням, зона репарации была бы полностью заполнена активными клетками, способными к дифференциации [21, 23].
На рис. 2 показаны графики распределения концентрации клеток в области длиной 21 = 10 мм, рассчитанные по формуле (37) при различных значениях коэффициента диффузии и в различные моменты времени. Как следует из проведенного анализа (рис. 2, а), при J = 0,67 мм2/день за клинически достоверный период времени Т = 120 суток костная мозоль будет полностью заполнена активными дифференцирующимися клетками. Найденное значение коэффициента диффузии используется в дальнейшем расчете. Отметим также, что относительная концентрация клеток в рассматриваемой зоне структурной перестройки кости достаточно быстро стремится к единице (рис. 2, б). Так, уже на шестидесятый день в среднем сечении ф = 0,976, это означает, что концентрация клеток составляет почти 100% от своего максимального значения.
-4-2 0 2 4
х, мм
а
х, мм б
Рис. 2. Распределение относительной концентрации активных клеток у вдоль продольной оси х зоны репарации: а - в конечный момент времени ? = 120 суток (дней) при различных значениях коэффициента диффузии: 3 = 0,01 мм /день (••••),
3 = 0,1 мм2/день (-----), 3 = 0,67 мм2/день (—); б - в различные моменты времени
при фиксированном значении коэффициента диффузии 3 = 0,67 мм2/день:
I = 1 день (••••), I = 30 дней (-------------------), I = 60 дней (—)
Результаты численного анализа и обсуждение
Рассмотренная одномерная модель структурной перестройки была использована для тестирования общего алгоритма трансформации ткани при действии механического стимулирующего воздействия. Исследуемая область, моделирующая костную мозоль между двумя отломками трубчатой кости, представляла собой полый цилиндр, высота которого - 21, внешний диаметр - ё' и внутренний - ё ”. Поскольку в качестве основной цели численного анализа было понимание особенностей механизма репарации костной ткани, то размер костной мозоли был выбран достаточно произвольно, главным образом для удобства представления: 21 = 10 мм, ё' = 20 мм, ё" = 14 мм.
Численный анализ был проведен при значениях материальных констант, представленных в таблице. Плотность двухфазной среды, модули Юнга и сдвига в дренированном состоянии Е и О, коэффициенты Био эффективных напряжений а, гидростатическая константа Я получены расчетным путем исходя из типичных значений упругих модулей материала твердой и жидкой фаз и пористости материала [2, 21].
Соотношения (21) для параметров, характеризующих проницаемость биологических тканей, удобно переписать в виде, используемом в дальнейших расчетах:
к (ш) =т
к
+ ішКр г т/ф ’
/шр К
к=А,
(38)
где значения К взяты из [21], а характеристика искривленности поровых каналов Т принята равной 1,66, что часто используется в расчетах и соответствует порам неспецифической формы [16].
Амплитудные значения деформаций и скорости потока жидкости в порах, перемещений и давления жидкости в порах (29) используются для расчета механорегулирующего индекса (2), являющегося критерием структурной перестройки костной ткани. Графическая интерпретация уравнения (2) представлена на рис. 3. В проведенном численном анализе использованы экспериментальные параметры, взятые из [21, 23]. Эмпирические константы: а = 0,0375, Ь = 0,003 мм/с. Пороговые значения механорегулирующего индекса следующие: 1) М > 3 - образование
фиброзной ткани; 2) 1 < М < 3 - образование хрящевой ткани; 3) 0,267 < М < 1 -образование незрелой костной ткани с достаточно высокой пористостью, близкой к губчатому веществу кости; 4) 0,01 < М < 0,267 - образование зрелой костной ткани, приближающейся по своим механическим характеристикам к компактному веществу кости; 5) М < 0,01 - резорбция костной ткани, возникающая при недостаточном механическом стимулировании.
Эффективные модули биологических тканей
Тип ткани а р, кг/м3 а С а С Я, Па К, м4/Нс
Гранулированная 0,99 1,000 1021 1,36-105 0,57-105 2,29-109 1,0-10-14
Фиброзная 0,80 0,990 1100 1,15 -106 0,47-106 0,21 -109 1,0-10-14
Хрящевая 0,80 0,995 1120 5,82-106 2,35 -106 1,07-109 5,0-10-15
Незрелая кость 0,65 0,893 1182 3,73 -109 0,97-109 1,42-109 1,0-10-13
Зрелая кость 0,20 0,435 1416 1,31-1010 0,45-1010 0,39-109 3,7-10-13
Компактное вещество кости 0,10 0,367 1468 2,11 -1010 0,63-1010 0,19-109 1,0-10-17
д, мм/с
Рис. 3. Области значений переменных, стимулирующих образование различных типов тканей при механическом воздействии на зону структурной перестройки
кости
1
0,5
тах
0
20 40 60 80 100 120
і, дни
Рис. 4. График зависимости изменения амплитуды приложенной силы от времени
Несмотря на определенные ограничения, построенная математическая модель изменяющейся сплошной среды позволяет исследовать влияние характеристик прикладываемой нагрузки на структурную перестройку костной ткани. Известно, что если действующая на кость сила слишком велика, то процесс нестабилен и со временем образуется лишь фиброзная ткань. Поэтому рекомендуется проводить расчет с небольших значений нагрузки, постепенно увеличивая амплитуду силы.
В качестве распределенной нагрузки, действующей на ткань в зоне репарации, рассмотрим сумму медленно изменяющейся сжимающей силы, постоянной на каждом шаге по времени, и быстрой составляющей, изменяющейся по гармоническому закону с относительно высокой частотой ш. Пусть в течение времени 0 < / < ^ статическая
нагрузка на кость плавно увеличивается до своего максимального значения ¥та& = 500 Н
и сохраняет достигнутую величину до конца периода репарации / = 120 суток (рис. 4). Поскольку в формулы (19) входит распределенная нагрузка, то пусть для примера ее интенсивность не зависит от координаты, / = ^ (218), где £ - площадь поперечного
сечения цилиндрической области, моделирующей костную мозоль. Соотношения (28), (29) могут быть использованы при статической нагрузке /, если считать, что ш = 0. Однако заметим, что в этом случае поток жидкости в порах отсутствует: 7 = 0 .
Будем считать для определенности, что амплитуда гармонической нагрузки /Ух пропорциональна постоянной составляющей с некоторым заданным коэффициентом, принятым в расчете Р = 0,1. Тогда общее выражение приложенной силы на каждом временном шаге п можно представить в виде
Решение от действия суммарной силы (37) будет являться суперпозицией статического и колебательного решений согласно формулам (28), (29). Для деформаций упругого скелета двухфазного материала и потоков поровой жидкости, входящих в выражение биомеханического критерия (2), можно записать окончательные соотношения:
(39)
(х,ю) = |б£) (X)| + |ё(”) (X,и)|, где 8^ (X) = 8(") (Х,0),
Ям (х ю) = |<?( п) (х ю)|-
(40)
Полученные формулы (16), (18), (29), (37), (40) были реализованы
в компьютерной программе на языке высокого уровня Еог^ап. Результаты расчета структурной перестройки ткани при вариации некоторых ключевых параметров построенной модели представлены на рис. 5-8. В качестве основной переменной, характеризующей процесс восстановления костной ткани, взят продольный модуль упругости Е . Расчеты осуществлялись в произвольных точках области репарации при шаге по времени Дt = 1 день.
Графики изменения модуля упругости (см. рис. 5) в характерной точке х = 2,5 мм для типичной нагрузки при ^ = 60 дней (см. рис. 4) показывают, что существенное влияние на сходимость результатов к тому или иному значению модуля Юнга и, следовательно, на образование определенного вида ткани оказывает частота приложенной гармонической нагрузки. При этом зависимость не является монотонной, что позволяет предположить достаточно сложную зависимость механизма перестройки от частотных свойств биомеханической системы. Из рисунка видно, что на низких частотах 1-10 Гц при удовлетворительном восстановлении в начальный период в последующем наблюдаются резорбция ткани и переход в ткань с низким модулем упругости. Только статическая нагрузка, не способная генерировать вынужденные потоки внутритканевой жидкости в системе поровых каналов, как следует из рис. 4, приводит к образованию фиброзно-хрящевой ткани. Интересно заметить, что высокая частота, равная 30 Гц для данной числовой модели, также не способствует удовлетворительному остеогенезу и достаточно быстро ведет к формированию фиброзной ткани. Оптимальной в данном случае является нагрузка частотой 20 Гц, стимулирующая устойчивый процесс трансформации исходного недифференцированного вещества в зрелую костную ткань.
Е, кПа
1,2-107
4-106
8-10'
I6
0
30
60
90
120
^ дни
-*-10 Гц
—®—20 Гц -К-30 Гц
Рис. 5. Графики изменения модуля упругости перестраивающейся ткани
в характерной точке х = 2,5 мм при различных значениях частоты действующей
нагрузки и ^ = 60 дней
Е, кПа
Рис. 6. Графики изменения модуля упругости перестраивающейся ткани
в различных точках костной мозоли при частоте действующей нагрузки 20 Гц
и ^ = 60 дней
х, мм
Ь = 60 дней X ^ = 30 дней □ tl = 120 дней
Рис. 7. Результирующий профиль распределения модуля упругости вдоль продольной координаты при частоте действующей нагрузки 20 Гц и ^ = 60 дней
Е, кПа
г, дни
О г1 = 15 дней X г1 = 30 дней А ^ = 45 дней ^ = 60 дней [ = 120 дней
Рис. 8. Графики изменения модуля упругости перестраивающейся ткани
в характерной точке х = 2,5 мм при частоте действующей нагрузки 20 Гц и
значениях г1= 15-120 дней
Следующие графики (см. рис. 6) демонстрируют зависимость механизма репарации от координаты при фиксированной частоте 20 Гц, принятой оптимальной для данной числовой модели. Распределение деформаций упругого скелета твердой фазы и потоков жидкости в порах может иметь достаточно сложный характер даже для области простой формы, не говоря уже о трехмерных задачах, что приводит к значительным изменениям механорегулирующего индекса (2) и попаданию текущей точки (в, д) в ту или иную область управляющей диаграммы (см. рис. 3). Так, в среднем сечении костной мозоли х = 0 мм поток жидкости (29) равен нулю, что обусловлено симметричными граничными условиями (23), а это, как и в случае только статической нагрузки, приводит к образованию фиброзно-хрящевой ткани. Отрицательный эффект низких скоростей жидкости в порах остается заметным до сечений х = ±1 мм, что составляет 20% костной мозоли.
Результирующий профиль распределения модуля упругости вдоль продольной координаты показан на рис. 7. Кроме центральной области вблизи нуля небольшое снижение модуля упругости заметно на краю х = ±5 мм.
На этом же рисунке показаны значения модуля Юнга при вариации параметра г1, определяющего время, за которое приложенная статическая сила достигает своего максимума. Как видно из рис. 7, широкие пределы изменения г1 от 30 до 120 дней не
повлияли на конечные значения модуля упругости. Однако более детальное исследование выявило серьезные проблемы в средний период восстановления кости при значениях г1 = 15...45 дней (см. рис. 8). Из графиков следует, что раннее
нагружение костной мозоли силой, близкой к собственному весу пациента, в начальный период приводит к более быстрому восстановлению упругих свойств регенерата. Однако образующаяся ткань демонстрирует нестабильное поведение и при превышении силой некоторого предела, характерного для каждого случая нагружения,
в области перестройки возникают обратные процессы, приводящие к формированию ткани с характеристиками незрелой кости. Интересно отметить, что при любом из рассмотренных значений параметра г1 процесс регенерации все же возобновляется при
г = 70 дней, являющемся, видимо, характерным временем для рассмотренной комбинации параметров модели.
Заключение
В статье представлены обобщенная динамическая модель изменяющейся пороупругой сплошной среды и математический алгоритм, концептуально описывающий процесс структурной перестройки костной ткани под действием внешнего механического стимула периодического характера. Математическая модель дает возможность исследовать процессы восстановления поврежденных костных элементов опорно-двигательного аппарата человека при наличии динамической нагрузки и теоретически обосновать выбор оптимального периодического воздействия на поврежденные ткани с целью их скорейшего и устойчивого заживления.
Рассмотренная одномерная модель структурной перестройки костной ткани и ее компьютерная реализация были использованы для тестирования общего достаточно сложного алгоритма и оценки влияния отдельных физико-механических параметров модели на процесс регенерации костной мозоли. В частности, построенная модель позволила исследовать влияние частоты стимулирующей нагрузки на процесс перестройки ткани, что совершенно отсутствует в известных источниках, а также влияние раннего нагружения на восстановление упругих свойств костной мозоли.
Конечно, необходимо понимать, что поскольку аналитическое решение уравнений пороупругости справедливо только при постоянных коэффициентах, то модель структурной перестройки дает в основном качественный результат. Тем не менее полученные численные результаты представляются достаточно реалистичными и соответствующими известным медицинским исследованиям процессов регенерации костной ткани в зоне перелома. Разработанная одномерная модель также может быть использована для тестирования более сложных программ на основе метода конечных элементов, реализующих представленный алгоритм для пространственной задачи изменяющейся пороупругой среды.
Благодарности
Работа выполнена при поддержке гранта РФФИ № 12-01-00054-а.
Список литературы
1. Акулич Ю.В. Математическая модель процесса внутренней адаптационной перестройки спонгиозной и кортикальной костных тканей человека // Механика композиционных материалов и конструкций. - 2005. - Т. 11, № 2. - С. 157-168.
2. Арсеньев Д.Г., Зинковский А.В., Маслов Л.Б. Эффективные упругие характеристики анизотропной модели пористого биологического материала, насыщенного жидкостью // Научно-технические ведомости СПбГПУ. - 2008. - № 3 (59). - С. 230-236.
3. Маслов Л.Б. Математическое моделирование колебаний пороупругих систем: моногр. - Иваново: Изд-во Иванов. гос. энергет. ун-та, 2010. - 264 с.
4. Нигматулин Р.И. Основы механики гетерогенных сред. - М.: Наука, 1978. - 336 с.
5. Остеопороз / под ред. А.И. Воложина, В.С. Оганова. - М.: Практическая медицина, 2005. - 238 с.
6. Регирер С.А., Шадрина Н.Х. Движение крови и интерстициальной жидкости в костной ткани: обзор // Известия РАН: Механика жидкости и газов. - 1999. - № 5. - С. 4-28.
7. Регирер С.А., Штейн А.А. Методы механики сплошной среды в применении к задачам роста и развития биологических тканей // Современные проблемы биомеханики. - 1985. - Вып. 2. - С. 5-37.
8. Штейн А.А. О континуальных моделях растущего материала // Механика композитных материалов. - 1979. - № 6. - С. 1105-1110.
9. Biot M.A. General theory of three-dimensional consolidation // J. Appl. Phys. - 1941. - Vol. 12, No. 2. -P. 155-164.
10. Boccaccio A., Prendergast P.J., Pappalettere C., Kelly D.J. Tissue differentiation and bone regeneration in
an osteotomized mandible: a computational analysis of the latency period // Med. & Biolog. Eng. &
Comp. - 2008. - Vol. 46, No. 3. - P. 283-298.
11. Burr D.B., Martin R.B. Mechanisms of bone adaptation to the mechanical environment // Triangle. - 1992. -Vol. 31, No. 2/3. - P. 59-76.
12. Carter D.R. Mechanical loading history and skeletal biology // J. Biomech. - 1987. - Vol. 20. -P. 1095-1109.
13. Carter D.R., Beaupre G.S., Giori N.J., Helms J.A. Mechanobiology of skeletal regeneration // Clin. Orthop. - 1998. - Suppl. 355. - S. 41-55.
14. Carter D.R., Wong M. The role of mechanical loading histories in the development of diarthrodial joints //
J. Orthop. Res. - 1988. - Vol. 6. - P. 804-816.
15. Claes L.E., Heigele C.A. Magnitudes of local stress and strain along bony surfaces predict the course and type of fracture healing // J. Biomech. - 1999. - Vol. 32. - P. 255-266.
16. Coussy O. Poromechanics. - New York: Wiley, 2004. - 315 p.
17. Giori N.J., Ryd L., Carter D.R. Mechanical influence on tissue differentiation at bone-cement interfaces // J. Arthroplasty. - 1995. - Vol. 10. - P. 514-522.
18. Goodship A.E., Lawes T.J., Rubin C.T. Low-magnitude high-frequency mechanical signals accelerate and augment endochondral bone repair: Preliminary evidence of efficacy // J. Orthop. Res. - 2009. - Vol. 27. -№ 7. - P. 922-930.
19. Huiskes R., van Driel W.D., Prendergast P.J., Soballe K. A biomechanical model for periprosthetic fibrous-tissue differentiation // J. Mat. Sci.: Materials in Medicine. - 1997. - Vol. 8. - P. 785-788.
20. Isaksson H., Comas O., van Donkelaar C.C., Mediavilla J., Wilson W., Huiskes R., Ito K. Bone regeneration during distraction osteogenesis: Mechano-regulation by shear strain and fluid velocity // J. Biomech. -2007. - Vol. 40, No. 9. - P. 2002-2011.
21. Isaksson H., Wilson W., van Donkelaar C.C., Huiskes R., Ito K. Comparison of biophysical stimuli for mechano-regulation of tissue differentiation during fracture healing // J. Biomech. - 2006. - Vol. 39, No. 8. - P. 1507-1516.
22. Kenwright J., Gardner T. Mechanical influences on tibial fracture healing // Clin. Orthop. Rel. Res. -1998. - Vol. 355S, No. 10. - P. 179-190.
23. Lacroix D., Prendergast P.J. A mechano-regulation model for tissue differentiation during fracture healing: analysis of gap size and loading // J. Biomech. - 2002. - Vol. 35, No. 8. - P. 1163-1171.
24. Liu X., Niebur G.L. Bone ingrowth into a porous coated implant predicted by a mechano-regulatory tissue differentiation algorithm // Biomech. & Modeling in Mechanobiology. - 2008. - Vol. 7, No. 4. -P. 335-344.
25. Pauwels F. A new theory concerning the influence of mechanical stimuli on the differentiation of the supporting tissue // Biomechanics of the locomotor apparatus / eds. P. Maquet, R. Furlong. - Berlin: Springer-Verlag, 1980. - P. 375-407.
26. Prendergast P.J. Finite element models in tissue mechanics and orthopaedic implant design // Clin. Biomech. - 1997. - Vol. 12, No. 6. - P. 343-366.
27. Prendergast P.J., Byrne D.P., Kelly D.J., Lacroix D. Tissue differentiation steps towards bone formation // Europ. Cells and Mat. - 2007. - Vol. 14, No. 1. - P. 26.
28. Prendergast P.J., Huiskes R., Soballe K. Biophysical stimuli on cells during tissue differentiation at implant interfaces // J. Biomech. - 1997. - Vol. 30, No. 6. - P. 539-548.
29. Van der Meulen M., Huiskes R. Why mechanobiology? A survey article // J. Biomech. - 2002. - Vol. 35, No. 4. - P. 401-414.
30. Wolff J. Das Gesetz der Transformation der Knochen. - Berlin: Verlag von August Hirschwald, 1892. -300 s.
MATHEMATICAL MODEL OF THE BONE STRUCTURAL TRANSFORMATION
L.B. Maslov (Ivanovo, Russia)
The main research approaches for solving the basic problem of the living tissue biomechanics being in development of mathematical theory of the hard tissue structural transformation determined by cell differentiation and governed by external force are considered. The general dynamic model of the changing poroelastic continuous medium and the mathematical algorithm conceptually describing the bone tissue structural transformation process under impact of the external mechanical stimulus of periodic behavior are introduced in the paper. The model implies that forced interstitial fluid flows induced in the bone pore system due to osseous matrix deformation along with inherent elastic matrix strains are the key mechanical factors regulating bone tissue reparative regeneration. For numerical analysis of the developed mathematical model, the computational algorithms to develop the software for computer simulation of the bone regeneration are suggested. The mathematical model provides possibility of the investigation of the regeneration processes of the damaged bone elements of the human locomotion system with the availability of a dynamic load and the theoretical argumentation of the choice of the optimal periodic impact to the damaged tissues for the fastest and stable healing.
Key words: bone tissue, structural transformation, mechanical stimulus, poroelasticity, vibration, mathematical modelling.
Получено 16 января 2013
