CC BY
96
DOI: 10.15593/RZhBiomeh/2016.3.03 УДК 616.314-08
к Российский
Журнал / Биомеханики
www.biomech.ru
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ КЛЕТОЧНЫХ ПРЕОБРАЗОВАНИЙ
ПРИ РЕГЕНЕРАЦИИ КОСТНОЙ ТКАНИ В УСЛОВИЯХ ИЗМЕНЯЮЩЕЙСЯ БИОХИМИЧЕСКОЙ СРЕДЫ С ВОЗМОЖНОЙ
МЕХАНОРЕГУЛЯЦИЕЙ
1 Кафедра травматологии, ортопедии и военно-полевой хирургии Ивановской государственной медицинской академии, Россия, 153012, Иваново, Шереметевский проспект, 8, e-mail: doc.kirpichev@yandex.ru
2 Кафедра высшей математики Ивановского государственного энергетического университета имени В.И. Ленина, Россия, 153003, Иваново, ул. Рабфаковская, 34, e-mail: dmitriyikorovin@list.ru, nikolay.tomin@gmail.com
3 Кафедра теоретической и прикладной механики Ивановского государственного энергетического университета имени В.И. Ленина, Россия, 153003, Иваново, ул. Рабфаковская, 34, e-mail: maslov@tipm.ispu.ru
Аннотация. Предлагается феноменологическая модель, описывающая возникающие при регенерации костной мозоли процессы диффузии, пролиферации, дифференциации и апоптоза основных типов клеток и тканевых структур под воздействием биохимических и механических факторов роста. Дополнительно к традиционно рассматриваемым четырем типам клеток (мезенхимальные стволовые клетки, фибробласты, хондроциты и остеобласты) и трем типам тканей (фиброзная, хрящевая и костная) в математическую модель включены кровеносные сосуды как четвертый тип ткани. В рамках принятой одномерной модели кость рассматривается как цилиндр с постоянным поперечным сечением, при этом явно выделяется зона повреждения. Все характеристики кости усредняются по поперечному сечению и изменяются только по длине кости. Основная система восьми дифференциальных уравнений в частных производных совместно с соответствующими граничными и начальными условиями описывает изменение концентраций указанных выше восьми типов клеток и тканей во времени для всех поперечных сечений кости. Вводится новый вид коэффициентов системы, определяющих миграцию клеток, взаимодействие и преобразование рассматриваемых типов клеток и тканей, в виде линейных комбинаций временно-зависимых функций, описывающих механические и биохимические факторы роста. Излагается план численного решения начально-краевой задачи для основной системы уравнений. Математическая модель дает возможность изучать процесс восстановления поврежденных костных элементов опорно-двигательного аппарата человека с целью нахождения оптимальных параметров стимулирования и медикаментозного воздействия на поврежденные ткани для их скорейшего и устойчивого заживления.
Ключевые слова: костная ткань, регенерация, клетка, диффузия, дифференциация, биомеханический фактор роста, морфогенетические белки кости, математическое моделирование.
© Кирпичев И.В., Коровин Д.И., Маслов Л.Б., Томин Н.Г., 2016
Кирпичев Иван Владимирович, к.м.н., доцент, завкафедрой травматологии, ортопедии и военно-полевой хирургии, Иваново
Коровин Дмитрий Игоревич, д.э.н., к.ф.-м.н., доцент, завкафедрой высшей математики, Иваново Маслов Леонид Борисович, д.ф.-м.н., доцент, завкафедрой теоретической и прикладной механики, Иваново
И.В. Кирпичев1, Д.И. Коровин2, Л.Б. Маслов3, Н.Г. Томин2
Введение
Теме моделирования процесса восстановления кости посвящены многочисленные работы. Особое внимание в них уделяется исследованию отношений между механической стимуляцией и преобразованием клеток и клеточных структур в процессе регенерации кости. Наиболее полноценной и адекватной для достижения целей управления механической стимуляцией процессов заживления костных тканей авторы считают модель, предложенную в [16]. Теоретическое исследование восстановления механических свойств ткани в зоне перелома с использованием осесимметричной конечно-элементной модели фрагмента кости было проведено с помощью уравнения диффузии, описывающего пространственно-временное распределение концентрации мезенхимальных стволовых клеток.
В работе [13] указывается на необходимость исследовать процессы миграции, пролиферации и дифференциации клеток различных типов. Предложенная в [13] модель была разработана для описания изменения концентраций фибробластов, хондроцитов и остеобластов, а также временного и пространственного распределения фиброзной (соединительной), хрящевой и костной тканей. Деятельность четырех типов клеток: мезенхимальных стволовых клеток, фибробластов, хондроцитов и остеобластов, взаимодействующих в зависимости от механической стимуляции, характеризовалась преобразованием друг в друга и в структуры тканей.
Кроме статьи [13], развивающей диффузионную модель [16] распространения активных протоклеток в область костной мозоли, возникающей после перелома, стоит упомянуть современные исследования, посвященные влиянию так называемых факторов роста на процессы биохимических преобразований мезенхимальных стволовых клеток в дифференцированные типы клеток [27, 28] при регенерации костной ткани. В данных работах подчеркивается их значительное стимулирующее влияние на пролиферацию, миграцию и дифференциацию клеток во время заживления кости. Можно отметить на основе собственных исследований [3, 19], что без учета факторов роста, активирующих процессы движения и преобразования клеток, представляется проблематичным построение адекватной модели, корректно описывающей восстановление костной ткани в объеме скаффолдов, применяемых в тканевой инженерии для замещения крупных дефектов длинных трубчатых костей человека. Нормальные значения скоростей процессов миграции, дифференциации и пролиферации клеток, используемые в модели регенерации костной мозоли [13, 16], без биохимического стимулирования могут оказаться недостаточными для заполнения клеточным материалом порового пространства имплантата значительного размера.
Математические модели, описывающие процессы преобразования клеток, в первую очередь дифференциацией мезенхимальных стволовых клеток, не ограничиваются только дифференциальными уравнениями типа приведенных в [13]. Активно развиваются как детерминистские цепные модели [21, 30], так и стохастические подходы [22] на их основе.
Уточнению модели [13, 16] были посвящены работы, учитывающие дополнительные ключевые аспекты заживления перелома кости, такие как направленная миграция клеток и развитие сети кровеносных сосудов [10, 23]. Действительно, существует тесная взаимосвязь между временем репарации кости и формированием кровеносных сосудов [9, 29]. Восстановление кровоснабжения признается важнейшим компонентом регенерации кости в работах [11, 24]. В [25] указывается, что рост новых капилляров кровеносных сосудов и формирование зрелой сосудистой сети имеют большое значение для дифференциации клеток в остеобласты и формирования костной ткани.
В настоящей работе будет рассмотрено обобщение предложенных в [10, 13] подходов к описанию процессов миграции клеток на основе уравнения диффузии в одномерном случае [3], клеточным преобразованиям и возникновению дифференцированных типов тканей с введением в рассмотрение моделей факторов роста.
Математическая модель преобразования клеток и тканей
Будем рассматривать диффузионную модель, учитывающую миграцию, пролиферацию, дифференциацию и апоптоз следующих восьми типов клеток и тканей, для которых введем сквозную нумерацию. Во-первых, четырех типов клеток: 1) мезенхимальных стволовых; 2) фибробластов; 3) хондроцитов; 4) остеобластов. Во-вторых, четырех типов тканей: 5) фиброзной (соединительной); 6) хрящевой; 7) костной и 8) кровеносных сосудов.
Перейдем к описанию предлагаемой математической модели восстановления кости, представляемой как динамический процесс установления нормальных концентраций (т.е. концентраций, присущих здоровой кости) для указанных выше восьми типов клеток и тканей, г = 1...8 .
Будем рассматривать кость как цилиндр длины I по оси х с постоянным поперечным сечением. Предполагается, что все характеристики кости усреднены по поперечному сечению и зависят только от положения сечения на отрезке [0,1] и
времени t > 0. Пусть начальный момент t = 0 соответствует повреждению кости, возникающему, например, в результате механической травмы, перелома кости или болезни.
Ограничимся рассмотрением указанных выше восьми основных типов г = 1...8 клеток и тканей (для краткости будем писать «клеток»), составляющих кость. Через С (х, t) обозначим объемную концентрацию клеток /-го типа в поперечном сечении х в
момент I. Пусть ci (х) - концентрация клеток /-го типа в неповрежденной кости в поперечном сечении х. Будем считать эту концентрацию известной при любом г = 1...8 на всем отрезке [0,1] и для определенности будем называть ее нормальной; в [13]
указаны максимальные концентрации различных типов клеток, и в качестве нормальной концентрации клеток того или иного типа в первом приближении можно взять половину от их максимальной концентрации.
Пусть [а, Ь] - зона повреждения кости по ее длине, [с, ё] - возможно, несколько более широкий участок, на котором явно выражено влияние зоны повреждения и может образоваться костная мозоль (0 < с < а < Ь < ё < I). Начальная
концентрация ^ (х) клеток /-го типа рассматривается как сумма соответствующей нормальной концентрации сг (х) и некоторой вызванной повреждением кости возмущающей функции ^ (х), где ^ (х) равна 0 вне отрезка [с, ё], носит хаотический
нерегулярный характер на [с, ё] и удовлетворяет неравенству | ^ (х)ёх < 0,
С
отражающему возможность потери костного материала во время повреждения кости.
Будем считать, что концы х = 0, х = I цилиндра, моделирующего кость, и его боковая поверхность на участке [с, ё] полупроницаемы, что свидетельствует о действии закона Ньютона для диффузии клеток каждого типа, т.е. поток внутрь объема
цилиндра пропорционален разности между концентрациями этого типа клеток в среде и на соответствующих участках поверхности цилиндра. В результате получаем неоднородные граничные условия третьего рода в сечениях х = 0, х = I, а для участка [с, ё] в правой части основного уравнения возникает дополнительное слагаемое.
Отметим, что в случае, когда коэффициент проницаемости стенок в сечениях х = 0, х = I обращается в 0 (случай непроницаемости кости в граничных сечениях), мы получаем однородные граничные условия второго рода.
Будем полагать, что на боковую поверхность кости на участке [с, ё] не
попадают клетки окружающих мягких тканей и не происходит прорастания соединительной ткани в пустоты, возможно имеющиеся в кости в результате механического повреждения, с последующим образованием рубцовой ткани, препятствующей нормальной регенерации кости. Для достижения этой цели при определенных обстоятельствах может быть поставлена остеоиндуктивная искусственная мембрана, отделяющая поврежденную поверхность кости на участке [с, ё] от окружающих мягких тканей, но пропускающая необходимые для регенерации
кости клетки. Кроме того, эта мембрана также защищает поврежденные участки кости от возможных локальных механических повреждений в процессе регенерации, а также она может при необходимости являться носителем морфогенетических белков и других факторов роста, способствующих росту клеток типов 1...4 и наносимых на внутреннюю поверхность мембраны. Важно, что этим процессом доставки морфогенетических белков можно управлять. Возможность этого также учитывается в рассматриваемой модели.
Выпишем основную систему восьми уравнений (1) относительно восьми неизвестных функций сг (х, 7) :
дс ^ d 2с С - ^ ( гЛ
— = D^—^r + a?rc dt ' дх2 11
1 - с,
v £rj
-Z
j*i
1 - ^ V
с
+
(1)
afCj 1 - f - a^c + kn ( ФП ( x, t)- сг) 1[cd]{ x)
j*i
v
г >0, 0<х<I, г = 1...8.
Здесь Д - коэффициент диффузии; арг (7) - функция, определяющая интенсивность пролиферации клеток; сг - нормальная концентрация клеток /-го типа в кости; а^ (г) - функция, определяющая интенсивность дифференциации /-го типа клеток в у'-й; а"р (г) - функция, определяющая интенсивность апоптоза клеток /-го типа; кл - коэффициент интенсивности проникновения клеток /-го типа из внешней среды (например из надкостницы) через боковую поверхность кости на участке [с, ё]; фг1 (х, 7) - концентрация клеток /-го типа вне кости вблизи от ее поверхности; ](х) - индикатор, принимающий значение 1, если х из [с, ё], и значение 0 в
противном случае. Величины арг, а^, а"13 являются неотрицательными функциями,
У
зависящими от времени I, а также от некоторых биохимических и механических факторов роста, в различной степени стимулирующих процессы пролиферации, дифференциации и апоптоза рассматриваемых типов клеток и количественно также зависящих от времени г. Умножение а^ на с в правой части уравнения (1) означает,
что скорость процесса дифференциации клеток типа у, отличного от г, в клетки типа г, пропорциональна имеющейся концентрации клеток типа у . При этом умножение на
с
относительный множитель 1 —обеспечивает уменьшение этой скорости при
£г
приближении концентрации сг клеток типа / к их нормальной концентрации сг. Знак
( ~ \
«плюс» перед слагаемым а^ с.
1 - ^ v CrJ
показывает, что такая дифференциация
способствует увеличению концентрации клеток типа г. Аналогично объясняется вид всех других слагаемых в правой части (1). В предпоследнем слагаемом правой части
с
(1), описывающем процесс апоптоза клеток типа г, множитель 1 —отсутствует,
с,
поскольку этот процесс обусловлен генетически и его скорость не зависит от близости концентрации сг к концентрации сг.
Подчеркнем, что при значениях индекса г = 5, 6, 7, 8 уравнения (1) определяют закон изменения плотностей тканей в результате возникновения и преобразования клеток соответствующих типов. При этом все коэффициенты диффузии должны быть равными нулю. Коэффициент, характеризующий апоптоз ткани, также обращается в нуль: ааР (7) = 0, в то время как два других коэффициента - арг (7) и а^ (7) -
понимаются как функции, определяющие интенсивность процессов роста и деградации массы (плотности) соответствующих костных структур.
Граничные условия третьего рода (2) для концентраций четырех типов клеток и начальные условия (3) для всех типов клеток и тканей примут следующий вид:
Ц ^ (0,7) + кг2 (с, (0,7) - ф,2 (7)) = 0, 7 > 0, г = 1...4, (2)
— I
Ц |с:(/,7) + к2 (с, (/,7)-ф(2 (7)) = 0, 7 > 0, г = 1...4,
С (х,0) = щ (х), 0 < х </, г = 1...8.
Отметим, что кг2 - это коэффициент проницаемости клеток /-го типа из внешней
среды внутрь кости через полупроницаемые граничные сечения х = 0 и х = /, а фг2 (7)
- концентрация клеток /-го типа во внешней среде вблизи этих граничных сечений. Стенки х = 0 и х = / становятся непроницаемыми при равенстве нулю коэффициента диффузии Ц = 0. Заданные функции щ (х) соответствуют характеру разрушения
тканей и определяют распределение концентраций клеток и плотности тканей в начальный момент времени.
Установить явный вид функций фг2 (7), фг1 (х, 7) можно только в хорошо
поставленном натурном эксперименте. При численном анализе для простоты и без потери общности можно задавать их равными константами, значения которых подбираются из виртуального эксперимента при многократном запуске модели для достижения наибольшего соответствия результатов расчетов реальному положению дел.
Очевидно, что система дифференциальных уравнений (1) содержит большое количество коэффициентов. Многие из этих параметров, в частности постоянные значения коэффициентов диффузии, скоростей процессов дифференциации, пролиферации и апоптоза клеток и соответствующих тканей, максимальные значения концентраций клеток и тканей, приведены в [13] на основе глубокого анализа
реферативных данных. Они могут быть использованы при решении системы (1) для моделирования регенерации ткани в области костной мозоли, возникающей при переломе кости вблизи поверхности имплантатов. Однако, как было отмечено во введении, в случае замещения дефектов кости пористыми скаффолдами значительных размеров постоянные значения интенсивностей арг, а^, аа вряд ли подходят и могут
приводить к неверным результатам, поскольку не учитывают биохимическую стимуляцию регенерации с помощью факторов роста. Конечно, и в других случаях лечения повреждений костей при использовании факторов роста вид коэффициентов системы (1) нуждается в уточнении. Поэтому следующие разделы статьи будут посвящены разработке коэффициентов системы как временно-зависимых функций, учитывающих биохимическую и механическую стимуляцию.
Описание биологических факторов роста
Регенерация костной ткани после повреждения представляет собой ряд последовательных стадий, включающих разрушение всех основных структур кости с формированием гематомы между отломками, воспалением и пролиферацией костной ткани, формированием соединительнотканной мозоли с последующей ее минерализацией и перестройкой [1, 27].
Каждая стадия сопровождается изменением концентрации биологически активных веществ и клеток в регенерате. Так, повреждение кости приводит к вазоконстрикции и снижению регионарного кровотока до 50 %. Данная физиологическая реакция направлена на снижение объема потери крови и уменьшение формирующейся гематомы. Сама гематома инициирует каскад биологических реакций, направленных на увеличение проницаемости сосудов и увеличение концентрации медиаторов воспаления и цинокинов, таких как ТКР-а, ГЬ-1, ГЬ-6, ГЬ-11, ГЬ-18, Ли§-1, Ли§-2. Это позволяет привлечь в зону регенерации полиморфно-ядерные лейкоциты, макрофаги, мезенхимальные стволовые клетки, начать прорастание новых капилляров, которые, в свою очередь, формируют сеть фибриновых, ретикулярных и коллагеновых волокон, стабилизирующих отломки.
Воспалительная реакция также инициирует активность остеокластов, которые удаляют некротизированные участки кости, подготавливая последнюю к заживлению. В последующем гематома видоизменяется в грануляционную ткань, являясь частью «мягкой» костной мозоли. Постепенно в регенерате повышается концентрация биологически активных веществ, направленных на образование хрящевой мозоли (ШЕ-р1, ШЕ-р2, ШЕ-р3, ОБЕ-5, ВМР-5, ВМР-6, ЮБ-1) путем пролиферации и дифференцировки мезенхимальных клеток. В последующем повышается концентрация биологических факторов роста (УЕОБ, ВМР-2, ЮТ-а), позволяющая сформировать незрелую грубую волокнистую костную ткань, которая в дальнейшем подвергается ремоделированию под влиянием механической среды.
Таким образом, процесс восстановления костной ткани после повреждения можно рассматривать как сложный взаимосвязанный каскадный биологический механизм, включающий остановку кровотечения, удаление омертвевших остатков костной ткани, привлечение в зону регенерата клеток, участвующих в восстановлении механической целостности и прочности поврежденной кости с последующей ее перестройкой. Можно сказать, что имеет место тонкая настройка структуры костного вещества в зависимости от механической нагрузки в каждой конкретной точке поврежденной кости. Это, в конечном итоге, позволяет вновь интегрировать поврежденный участок в общую структуру скелета. Маркерами начала и окончания той
или иной фазы можно считать изменение концентрации биологически активных веществ и клеток.
В табл. 1 приведены основные типы биохимических факторов роста, стимулирующих соответствующие стадии регенерации костной ткани.
Морфогенетические белки кости - это совокупность факторов роста кости, появляющихся только при объемных травмах и переломах. Это белковые молекулы-магниты TGF-P, притягивающие клетки кости в точки, нуждающиеся в восстановлении. Морфогенетические белки кости проходят три главных этапа последовательно с костным морфогенезом, а именно, хемотаксис, размножение и дифференциацию временной хрящевой основы, а также постоянную костную индукцию.
Внеклеточный белок TNF-а является многофункциональным провоспалительным цитокином, синтезирующимся в основном моноцитами и макрофагами. Влияет на липидный метаболизм, коагуляцию, устойчивость к инсулину, функционирование эндотелия. Трансформирующий ростовой фактор бета TGF-Ps -белок, который контролирует пролиферацию, клеточную дифференциацию и другие функции в большинстве клеток. Инсулиноподобный фактор роста IGF-1 - белок из семейства инсулиноподобных факторов роста, по структуре и функциям похожий на инсулин. В его функции входит эндокринная, аутокринная и паракринная регуляция процессов роста, развития и дифференциации клеток и тканей организма.
Таблица 1
Влияние факторов роста на биологические процессы
Фактор роста Стадия воспаления Хондрогенные процессы (образование хрящевой ткани) Образование и формирование костной ткани Кровеносные сосуды
BMP Ограничение воспалительного процесса (ЮТ-а) Морфогенез хряща, влияние на скелетные модели. Пролиферация и созревание хондроцитов, сильная индукция хондро-генной дифференциации мезенхимальных клеток Развитие миграции клеток, пролиферация и остеогенная дифференциация мезенхимальных клеток (BMP-2, TNF-a) [6, 8, 14, 17] -
TGF-Ps Миграция коллагена, остеонектина и остеопонтина, стимуляция недифференцированных мезенхимальных клеток [17, 20]
- Повышение пролиферации хондроци-тов. Стимуляция хондро-генной дифференциации мезенхимальных клеток Ингибирование пролиферации зрелых клеток и предшественников остеокластов [8] Стабилизация новых кровено^ ных сосудов
IGF-1 Посредник воспаления Хондрогенная дифференциация мезенхимальных стволовых клеток Производство костного матрикса [7]. Стимулирующие эффекты образования остеобластов [14] -
VEGF - - Обеспечение миграции и дифференциации остеобластов [12] Формирование транспорта крови [6, 11, 14], развитие сосудов
Эндотелиальный фактор роста сосудов VEGF - еще один из членов семейства структурно близких между собой белков. Он влияет на развитие новых кровеносных сосудов (ангиогенез) и выживание незрелых кровеносных сосудов (сосудистая поддержка), связываясь с двумя близкими по строению мембранными тирозинкиназными рецепторами и активируя их.
Математические подходы к моделированию факторов роста
Анализируя концентрацию факторов роста по данным, полученным в [15, 18, 26] и [27], можно выдвинуть следующую гипотезу. Поступление факторов роста в область повреждения кости подчиняется некоторому генетическому закону, который выделяет в зону повреждения количество фактора роста, необходимое для реализации процесса роста. При этом в нормальных условиях жизнедеятельности кости генерация факторов не ограничена в тех объемах, которые необходимы для процесса восстановления тканей, миграции, пролиферации и дифференциации клеток. Следовательно, можно предположить, что мы имеем простейшую модель распределения моментов поступления факторов роста в зону повреждения, известную в теории очередей как многоканальная система массового обслуживания со стационарным однородным потоком заявок, т.е. моментов поставки в зону «порции» соответствующего фактора роста [5]. В этом случае зависимость концентрации фактора роста от времени ФД^) необходимо рассматривать как функцию Эрланга вида
Ф (') = ^ • ^ • , (4)
где /., п., т . - параметры функции, отвечающей ] -му фактору роста, которые выбираются так, чтобы экспериментальные наблюдения были наиболее близки к значениям Ф. (7) .
На основе простейшей схемы оптимизации с использованием метода наименьших квадратов были получены аналитические выражения концентрации факторов роста, влияющих на развитие и подавление воспалительного процесса, процесса образования хрящевой ткани, процесса образования и формирования костной ткани. Рассчитанные параметры функций Эрланга представлены в табл. 2, а функциональные асимптотики соответствующих факторов роста - на рис. 1.
Таблица 2
Процесс образования и формирования костной ткани
№ п/п Процесс, аппроксимирующая функция Участие факторов роста Рассчитанные параметры модели Достоверность (коэффициент детерминации)
fj nj т j
1 Воспалительный, ф,(< ) ЮТ-а, IGF-1 16,663 2,3 2,875 0,685
2 Образование хряща, Ф 2 ( 7 ) ТОР-Рз, ВМР-5, ВМР-6, ЮР-1 0,00082 7 0,9333 0,845
3, 4 Окостенение (образование кровеносных сосудов), Фз (7), Ф4 (7) УБОР, ВМР-2, ТОТ-а 0,0028 3,5 0,2413 0,615
Концентрация факторов А
А
1 3
21
7 14
Время, день Воспалительный процесс Функциональная асимптотика ФКО
Концентрация факторов
1 3
7 14 21
Время, день Концентрация факторов TGF-Ps, ВМР-5, ВМР-6, IGF-1 Функциональная асимптотика Ф2(0 б
1
а
Концентрация факторов
Время, день Концентрация факторов VEGF, ВМР-2, ТОТ-а на стадии окостенения — — — Функциональная асимптотика Фз(0
в
Рис. 1. Модель динамики концентрации факторов роста, влияющих на развитие и подавление воспалительного процесса (а), на процесс образования хрящевой ткани (б), на процесс образования и формирования костной ткани (в)
Все функции Эрланга построены таким образом, чтобы их максимальные значения были равны единице. Тем самым нормированные на единицу функции Эрланга используются для задания временной зависимости, в то время как абсолютные значения представим с помощью масштабных коэффициентов (см. следующий параграф).
Процесс активного построения кровеносных сосудов, на который существенное влияние оказывает фактор VEGF, нами рассматривается как совпадающий с процессом образования и формирования костной ткани, что позволяет принять равенство
Ф4 ( 7) = Фз ( 7) .
Описание параметров модели
Рассмотрим функциональные зависимости между параметрами системы уравнений (1). Четыре первых уравнения системы описывают миграцию, пролиферацию, дифференциацию и апоптоз каждого из основных типов клеток г. Движение клеток в ткани моделируется с помощью диффузионного процесса. Предполагая, что изменения коэффициентов Д, кл и кг2 будут несущественны, положим их равными константам. В случаях г = 5, 6, 7, 8 имеем Д = 0, так как миграция сформированных тканей, естественно, отсутствует. В этих случаях кл = 0, кп = 0, так что первое и последнее слагаемые в правой части (1) и граничные условия (2) пропадают, а знание функций фг2 (7), фг1 (х, 7) не требуется.
Скорость изменения концентрации исследуемых клеток пропорциональна количеству фактора, инициирующего изменение этой концентрации. Кроме того, так как концентрация ограничена «нормальным» значением сг, процесс должен
стабилизироваться около сг. Следовательно, вклад влияния пролиферации на скорость
дс,.
изменения концентрации —можно оценить как
д7
ар (7) сг - . (5)
V Сг У
Аналогично для дифференциации одного типа клеток в другой можно получить упрощенные формулы. В частности, если г -й тип преобразуется в ] -й, то вклад в изменение скорости концентрации можно оценить как
- — U) с.
(' ) с
f
1 —
V ёи
(6)
Наоборот, если ] -й тип преобразуется в г -й, то вклад в скорость изменения концентрации примет вид
( ~ \
„Я (7 ) с
af (t) с 1 - Cl . (7)
V С,У
Рассмотрим схему, описывающую развитие костной ткани. На графе (рис. 2) выделены следующие состояния развития: 1) мезенхимальные стволовые клетки; 2) фибробласты; 3) хондроциты; 4) остеобласты; 5) соединительные ткани; 6) хрящевые ткани; 7) костные ткани и 8) кровеносные сосуды.
Стрелки, представляющие собой ориентированные ребра, соединяющие разные вершины, означают, что возможны дифференциация одного типа клеток в другой или преобразование (структурирование) клеток в ткань. Наличие петлеобразных стрелок, выходящих и входящих в одну вершину, указывает на процесс пролиферации или развития структуры.
В рассматриваемой математической модели каждому ребру приписана функция, определяющая интенсивность соответствующего процесса во времени. Так как присутствует большое количество нулевых функций (что означает отсутствие процесса преобразования из одной вершины в другую), то удобно указать список этих функций с помощью матрицы сопряженности
j = 1 j = 2 j = 3 j = 4 j = 5 j = 6 j = 7 j = 8
= 1 ai (t) а12 (t) ai .(t) а14 (t) 0 0 0 ai С)
= 2 0 a2 (t) 0 0 a25 (t) 0 a27 (t) 0
= 3 0 0 a3 (t) 0 0 a36 (t) 0 0
= 4 М (t ) = 0 0 0 a4 (t) 0 0 a47 (t) 0
= 5 0 0 0 0 a5 (t) 0 0 0
= 6 0 0 0 0 0 a6 (t) a67 (t) 0
= 7 0 0 0 0 0 0 a7 (t) 0
= 8 0 0 0 0 0 0 a87 (t) a8 (')
Элемент матрицы М({) в /-й строке и у-м столбце есть интенсивность реализации соответствующего преобразования и коэффициент в модели, описываемой уравнениями (1). Нули в матрице определяют равенство нулю соответствующего слагаемого в модели.
Опишем теоретические подходы к заданию ненулевых элементов матрицы сопряженности. Коэффициент а[г (7) соответствует интенсивности пролиферации
мезенхимальных стволовых клеток. Будем считать, что интенсивность пролиферации мезенхимальных стволовых клеток постоянна, поэтому имеет смысл принять этот коэффициент равным константе К, значение которой можно определить лабораторным путем.
Элемент (7) соответствует интенсивности дифференциации
мезенхимальных стволовых клеток в фибробласты. Процесс этой дифференциации плохо изучен, но можно предположить, что поскольку фибробласты активно участвуют в ремоделировании костной ткани путем синтеза компонентов межклеточного вещества, то этот процесс должен активизироваться еще на стадии воспаления. Поэтому полагаем
< (7) = Л2Ф1 (7) + (7), (8)
где Д2 - масштабный коэффициент; 512 (7) - функция, отвечающая за механическую стимуляцию процесса.
Отметим, что Л12 подбирается из сравнения теоретического прогнозирования с реальными данными для более точного описания процесса. Величина 512 (7) может
рассматриваться как управляющая переменная, в частности отвечающая за механическую стимуляцию процесса. Здесь и далее подобные переменные предполагаются кусочно-постоянными на временном отрезке функциями.
Коэффициент а^* (7) соответствует дифференциации мезенхимальных
стволовых клеток в хондроциты. Известно, что мезенхимальные стволовые клетки соединяются в хондрогенные островки, клетки которых размножаются и дифференцируются в хондробласты. Хондробласты затем дифференцируются в хондроциты. Существенно на процесс влияют факторы ТОР-Рз и 1ОР-1. Поэтому полагаем
< (7) = Аз (в 1зФ1 (^) + ЙФ (7)) + 5,3 (7) , (9)
где Аи, Р^, Р^ - масштабные коэффициенты; 513 (7) - функция, отвечающая за механическую стимуляцию процесса.
Коэффициент а^ (7) отвечает за процесс дифференциации мезенхимальных
стволовых клеток в остеобласты. В ходе прямого остеогенеза мезенхимальные стволовые клетки при образовании остеогенных островков непосредственно дифференцируются в остеобласты. Активное влияние на эти процессы оказывают все рассматриваемые нами факторы: ТОТ-а, 1ОР-1, ТОР^, ВМР-5, ВМР-6, УБОР, ВМР-2. Поэтому удобным представляется следующее выражение коэффициента:
< (^) = (&Ф (7) + Р?4Ф2 (7) + РЗ4Ф3 (7)) + 5,4 (7) , (10)
где А14, Р^, Р^, Р^ - масштабные коэффициенты; 514 (7) - функция, отвечающая за
механическую стимуляцию процесса.
Аналогично получим выражения коэффициентов, моделирующих действие фибробластов по созданию межклеточного вещества для организации костной ткани:
арг (7) = 42Ф1 (7), а* (7) = 45Ф2 (7), а* (7) = А7ф 2 ('). (11)
Коэффициенты, описывающие интерстициальный рост хряща, можно задать следующим образом:
ар (7) = А33Ф2 (7), а% (7) = 46Ф2 (7) . (12)
С другой стороны, аппозиционный рост хряща можно задать с помощью следующего выражения:
а? (7) = АФ (7). (13)
Образование остеобластов и дальнейший остеогенез опишем, используя соотношения
аРг (7) = Л4 (Р44Ф2 (7) + в44Фз (7)), а% (7) = Л47Ф3 (7), (14)
а учитывая существенное влияние факторов TGF-Ps, ВМР-5, ВМР-6, IGF-1 на развитие соединительных тканей, получим
а? (7) = 45Ф2 (7) . (15)
Результат непрямого остеогенеза предлагается описать следующим образом:
а6Т (7 ) = 47Ф3 (7). (16)
Функциональный коэффициент, который определяет активность таких биологических процессов, как минерализация костной мозоли, образование первичного губчатого вещества (ретикулофиброзной ткани), формирование надкостницы и наращивание остеоцитов посредством разрушения ретикулофиброзной костной ткани (создание пластинчатой ткани), зададим в виде
а7РГ (7) = А77Ф3 (7)-А77екрФ1 (7) , (17)
где Л7н7екр - подбираемый коэффициент, определяющий интенсивность разрушения костной ткани при некрозе на стадии воспаления.
Для последних коэффициентов матрицы сопряженности, характеризующих образование сети кровеносных сосудов, примем следующие функциональные зависимости:
а? (7) = АФ4 (7), а8Г (7) = (7) . (18)
Наконец будем предполагать, что интенсивность апоптоза существенно не зависит от процессов регенерации, а следовательно, можно установить для всех типов клеток и тканей простые выражения
ааР (7) = р, (19)
где р - константы.
План численного решения основной краевой задачи
Для значений г = 1...4 уравнения (1) относительно неизвестных функций сг,
г = 1...4 являются дифференциальными уравнениями в частных производных второго порядка вида «реакция-диффузия», которые нередко возникают в химической кинетике. Система (1) первых четырех уравнений (г = 1...4) в случае известных
функций с. (у = 5...8) при граничных условиях (2) и начальных условиях (3) может
быть численно решена методом прямых, или методом Роте [2]. При г = 5... 8 в системе (1), как было отмечено выше, коэффициенты Ц и кл должны быть равными нулю, что приводит к дифференциальным уравнениям в частных производных первого порядка.
Тогда система (1) последних четырех уравнений (г = 5...8) в случае известных функций с ., у = 1...4 при начальных условиях (3) может быть численно решена относительно неизвестных функций сг, г = 5...8 тем или иным сеточным методом [4].
Отсюда естественным образом возникает следующий план численного решения дифференциальной краевой задачи (1)-(3) с восьмью неизвестными функциями сг,
г = 1...8 (конкретизацию, обоснование и численную реализацию этого плана предполагается изложить в последующих работах авторов). Согласно (3) нам известны значения всех восьми неизвестных функций сг, г = 1...8 при 7 = 0 и всех х из [0,1]. Поэтому, применяя явную разностную схему (с шагом т по переменной 7 и с шагом И
по переменной х, где И = — и т - натуральное число) к системе последних четырех
т
уравнений (1), находим приближенные значения функций сг, г = 5...8 при 7 = т и всех х = кИ, где к = 0...т. Используя интерполяцию функций, получаем приближенные аналитические выражения для функций сг, г = 5...8 при 7 = т и всех х из [0,1]. Зная
эти аналитические выражения и применяя метод Роте [2] к первым четырем уравнениям (1), с учетом (2) и (3) находим приближенные аналитические выражения для функций с, г = 1. 4 при 7 = т и всех х из [0,1]. Теперь нам известны
приближенные аналитические выражения для всех восьми неизвестных функций с , г = 1...8 при 7 = т и всех х из [0,1]. Это позволяет нам, повторяя описанный только что вычислительный процесс с заменой 7 = 0 на 7 = т, получить приближенные аналитические выражения для функций сг, г = 1...8 при 7 = 2т и всех х из [0,1], затем
на их основе аналогично - при 7 = 3т, потом последовательно при 7 = 4т, 7 = 5т и т.д.
Предложенная в настоящей статье математическая модель описывает основные биологические процессы, происходящие при восстановлении или перестройке костной ткани. Полученные данные, с учетом использования в модели факторов роста, позволят смоделировать ряд патологических процессов, происходящих в костях при врожденных, эндокринных и других системных заболеваниях. Это позволяет применять ее не только в научных целях при описании патогенеза заболеваний, но и при моделировании состояния кости при использовании фармакологических средств и изделий медицинского назначения. Таким образом, после всесторонней верификации применение предложенной модели будет актуальным не только в травматологии и ортопедии, но и в эндокринологии, ревматологии, терапии, онкологии, фармакологии.
Заключение
В статье представлена одномерная модель процесса восстановления кости в условиях механической и биохимической стимуляции. В нее заложено большое количество исследуемых в работе разнообразных связей и зависимостей между основными типами клеток и тканей, составляющих кость и влияющих на процесс ее восстановления. Рассматриваются также вопросы, связанные с численной реализацией модели.
Ее реализация и применение требуют проведения тестовых вычислительных и натурных экспериментов, что выходит за рамки настоящей статьи. Представленные результаты служат основой для дальнейших теоретических исследований в достаточно новой интенсивно развивающейся области на стыке вычислительной биомеханики и тканевой инженерии, посвященной теоретическому изучению влияния биохимического
и механического воздействия на биологические процессы ремоделирования и репаративной регенерации ткани после переломов, эндопротезирования, применения перспективных пористых имплантатов.
Благодарности
Работа выполнена при поддержке гранта Российского фонда фундаментальных исследований № 15-29-04825.
Список литературы
1. Аврунин А.С. Остеоцитарное ремоделирование костной ткани: история вопроса, морфологические маркеры // Морфология. - 2011. - № 1. - С. 86-94.
2. Мартинсон Л.К., Малов Ю.И. Дифференциальные уравнения математической физики. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2002. - 368 с.
3. Маслов Л.Б. Математическая модель структурной перестройки костной ткани // Российский журнал биомеханики. - 2013. - Т. 17, № 2 (60). - С. 39-63.
4. Самарский А.А. Введение в теорию разностных схем. - М.: Наука-Физмат, 1971. - 553 с.
5. Хинчин А.Я. Работы по математической теории массового обслуживания. - М.: Физматгиз, 1963. -263 с.
6. Allori A.C., Sailon A.M., Warren S.M. Biological basis of bone formation, remodeling, and repair - part I: biochemical signaling molecules // Tissue Eng. Part B: Reviews. - 2008. - Rev. 14, № 3. - P. 259.
7. Arnsdorf E.J., Tummala P., Kwon R.Y., Jacobs C.R. Mechanically induced osteogenic differentiation - the role of RhoA, ROCKII and cytoskeletal dynamics // J. Cell Sci. - 2009. - Vol. 122, № 4. - P. 546-553.
8. Bose S., Tarafder S. Calcium phosphate ceramic systems in growth factor and drug delivery for bone tissue engineering: a review // Acta Biomater. - 2011. - Vol. 8, № 4. - P. 1401-1421.
9. Buschmann J., Welti M., Hemmi S., Neuenschwander P., Baltes C., Giovanoli P., Rudin M., Calcagni M. Three-dimensional co-cultures of osteoblasts and endothelial cells in DegraPolfoam: histological and high-field magnetic resonance imaging analyses of pre-engineered capillary networks in bone grafts // Tissue Eng. Part A. - 2011. - Rev. 17, № 3, 4. - P. 291-299.
10. Geris L., Gerisch A., Sloten J.V., Weiner R.D., Oosterwyck H.V. Angiogenesis in bone fracture healing: a bioregulatory model // J. Theor. Biol. - 2008. - Vol. 251. - P. 137-158.
11. Hankenson K.D., Dishowitz M., Gray C., Schenker M. Angiogenesis in bone regeneration // Injury. -2011. - Vol. 42, № 6. - P. 556-561.
12. Huang Y.C., Kaigler D., Rice K.G., Krebsbach P.H., Mooney D.J. Combined angiogenic and osteogenic factor delivery enhances bone marrow stromal cell-driven bone regeneration // J. Bone Miner. Res. -2005. - Vol. 20, № 5. - P. 848-857.
13. Isaksson H., van Donkelaar C.C., Huiskes R., Ito K. A mechano-regulatory bone-healing model incorporating cell-phenotype specific activity // Journal of Theoretical Biology. - 2008. - Vol. 252. -P. 230-246.
14. Kempen D.H., Creemers L.B., Alblas J., Lu L., Verbout A.J., Yaszemski M.J., Dhert WJ. Growth factor interactions in bone regeneration // Tissue Eng., Part B. - 2010. - Rev. 16, № 6. - P. 551-566.
15. Kon E., Delcogliano M., Filardo G., Pressato D., Busacca M., Grigolo B., Desando G., Marcacci M. A novel nano-composite multi-layered biomaterial for treatment of osteochondral lesions: technique note and an early stability pilot clinical trial // Injury. - 2010. - Vol. 41, № 7. - P. 693-701.
16. Lacroix D., Prendergast P.J., Li G., Marsh D. Biomechanical model to simulate tissue differentiation and bone regeneration: application to б^сШте healing // Med. Biol. Eng. Comput. - 2002. - Vol. 40, № 1. -P. 14-21.
17. Lieberman J.R., Daluiski A., Einhorn T.A. The role of growth factors in the repair of bone. Biology and clinical applications // J. Bone Joint Surg. Am. - 2002. - Vol. 84-A, № 6. - P. 1032-1044.
18. Liu C., Han Z., Czernuszka J.T. Gradient collagen/nanohydroxyapatite composite scaffold: development and characterization // Acta Biomater. - 2009. - Vol. 5. - P. 661-669.
19. Maslov L.B. Mathematical modeling of the callus mechanical properties restoration // J. Appl. Math. Mech. - 2015. - Vol. 79, № 2. - P. 195-206.
20. Mundy G.R., Chen D., Zhao M., Dallas S., Xu C., Harris S. Growth regulatory factors and bone // Rev. Endocr. Metab. Disord. - 2001. - Vol. 2, № 1. - P. 105-115.
21. Nakata Y., Getto P., Marciniak-Czochra A., Alarcon T. Stability analysis of multi-compartment models for cell production systems // J. Biol. Dyn. - 2012. - Vol. 6. - P. 2-18.
22. Pazdziorek P.R. Mathematical model of stem cell differentiation and tissue regeneration with stochastic noise // Bull. Math. Biol. - 2014. - Vol. 76, № 7. - P. 1642-1669.
23. Pivonka P., Dunstan C.R. Role of mathematical modeling in bone fracture healing // Bone Key Reports 1, Article number: 221. - 2012. - № 11. - P. 1-10.
24. Rhinelander F.W. Tibial blood supply in relation to fracture healing // Clin. Orthop. Relat. Res. - 1974. -Vol. 105. - P. 34-81.
25. Risau W. Mechanisms of angiogenesis // Nature. - 1997. - Vol. 386. - P. 671-674.
26. Rodrigues M.T., Gomes M.E., Reis R.L. Current strategies for osteochondral regeneration: from stem cells to pre-clinical approaches // Curr. Opin. Biotechnol. - 2011. - Vol. 22. - P. 726-733.
27. Santo V.E., Gomes M.E., Mano J.F., Rui L.R. Controlled release strategies for bone, cartilage, and osteochondral engineering. Part I: Recapitulation of native tissue healing and variables for the design of delivery systems // Tissue Engineering. Part B: Reviews. - 2013. - Vol. 19, № 4. - P. 308-326.
28. Santo V.E., Gomes M.E., Mano J.F., Rui L.R. Controlled release strategies for bone, cartilage, and osteochondral engineering. Part II: Recapitulation of native tissue healing and variables for the design of delivery systems // Tissue Engineering. Part B: Reviews. - 2013. - Vol. 19, № 4. - P. 327-352.
29. Santos M.I., Unger R.E., Sousa R.A., Reis R.L., Kirkpatrick C.J. Crosstalk between osteoblasts and endothelialcells co-cultured on a polycaprolactone-starch scaffold and the in vitro development of vascularization // Biomaterials. - 2009. - Vol. 30, № 26. - P. 4407-4415.
30. Stiehl T., Marciniak-Czochra A. Characterization of stem cells using mathematical models of multistage cell lineages // Math. Comput. Model. - 2011. - Vol. 53, № 7-8. - P. 1505-1517.
MATHEMATICAL MODEL OF CELL TRANSFORMATIONS AT BONE TISSUE REGENERATION UNDER ALTERATING BIOCHEMICAL MEDIUM WITH POSSIBLE MECHANOREGULATION
I.V. Kirpichev, D.I. Korovin, L.B. Maslov, N.G. Tomin (Ivanovo, Russia)
A phenomenological model is proposed which describes diffusion, proliferation, differentiation, and apoptosis processes appearing during regeneration of callus for the basic cell types and tissue structures under the impact of biochemical and mechanical growth promoting factors. Additionally to the traditionally considered four cell types (mesenchymal stem cells, fibroblasts, cartilage cells, and osteoblasts) and three tissue types (fibrous, cartilaginous and bone), blood vessels as the fourth tissue type are included to the mathematical model. In the context of the accepted one-dimensional model, a bone is considered as a cylinder with a constant cross-section, in doing so the damage area obviously is emphasized. All the bone characteristics are averaged by the cross-section and are changed along the bone length only. The basic set of eight differential equations in the partial derivatives along with the corresponding boundary and initial conditions describes changes in concentrations of the described above eight cell types and tissues through time for all cross-sections of a bone. A new system of coefficient types defining cell migration, interaction and transformation of the considered cell types and tissues in the form of linear combinations of the time-sensitive functions describing mechanical and biochemical growth promoting factors is intriduced. A numerical calculation plan for solving the initial boundary value problem for the basic set of equations is stated. The mathematical model gives the opportunity to study the process of regeneration of the damaged bone elements of the human locomotor system with the purpose of consequent identifying optimal parameters of simulation and medicated impact on the damaged tissues for their fastest and sustained healing.
Key words: bone tissue, regeneration, cell, diffusion, differentiation, biomechanical growth factor, bone morphogenetic proteins, mathematical modelling.
Получено 7 июля 2016
