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

Солодкова Е.Г.; Малюгин Б.Э.; Захаров И.Н.; Багмутов В.П.; Фокин В.П.; Балалин С.В.; Лобанов Е.В.; Лэ В.Х.

Солодкова Е.Г., Малюгин Б.Э., Захаров И.Н., Багмутов В.П., Фокин В.П., Балалин С.В., Лобанов Е.В., Лэ В.Х. Разработка комплекса математических моделей биомеханических параметров роговицы с диагностированным кератоконусом до и после лечения кросслинкингом роговичного коллагена. Российский журнал биомеханики, 2022, № 3, С. 10-28. DOI: 10.15593/RZhBiomeh/2022.3. 01

РОССИИСКИИ ЖУРНАЛ БИОМЕХАНИКИ № 3, 2022

RUSSIAN JOURNAL OF BIOMECHANICS

https://ered.pstu.ru/index.php/rib

Научная статья

DOI: 10.15593/RZhBiomeh/2022.3.01 УДК 531/534: [57+61]

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

Е.Г. Солодкова1, Б.Э. Малюгин2, И.Н. Захаров3, В.П. Багмутов3, В.П. Фокин1, С.В. Балалин1, Е.В.Лобанов1, В.Х. Лэ3

1 Волгоградский филиал Национального медицинского исследовательского центра, межотраслевой научно-технический комплекс «Микрохирургия глаза» им. акад. С.Н. Федорова Минздрава России, Россия

2 Национальный медицинский исследовательский центр, межотраслевой научно-технический комплекс «Микрохирургия глаза» им. акад. С.Н. Федорова Министерства здравоохранения Российской Федерации, Москва, Россия

3 Волгоградский государственный технический университет, Волгоград, Россия

О СТАТЬЕ

АННОТАЦИЯ

Получена: 22 января 2022 Одобрена: 05 сентября 2022 Принята к публикации: 09 сентября 2022

Ключевые слова:

математическое моделирование, роговица, кератоконус, кросслинкинг, метод конечных элементов, РеЫаоат, Оомв.

Целью работы является построение трехмерной геометрической и механической конечно-элементной модели роговицы и численного анализа ее напряженно-деформированного состояния и механического поведения в условиях действия внутриглазного давления и внешнего давления воздушной струи при диагностике пациентов на разных стадиях керато-конуса, а также после лечения с помощью кросслинкинга роговичного коллагена. Выстраивается геометрическая модель роговицы в виде пространственного сегмента выпуклой тонкостенной оболочки с переменной толщиной стенки и произвольной формой передней и задней поверхностей, задаваемых путем интерполяции экспериментальных данных томографического исследования роговицы конкретного пациента с помощью кератотопографа РеЫаоат ДХ1. Материал оболочки, моделирующей роговицу, считается неоднородным, изотропным и нелинейно-упругим. Его коэффициенты жесткости устанавливаются на основе известных экспериментальных данных и уточняются из сопоставления результатов численного моделирования параметров деформации роговицы при воздействии воздушного им-

© Захаров Игорь Николаевич - д.т.н., заведующий кафедрой, e-mail: sopromat@vstu.ru :0000-0001-7177-7245 © Багмутов Вячеслав Петрович - д.т.н., профессор кафедры, e-mail: sopromat@vstu.ru :0000-0003-3648-8450 © Фокин Виктор Петрович - д.м.н., директор филиала, e-mail: fokin@isee.ru :0000-0002-2513-9709 © Балалин Сергей Викторович - д.м.н., заведующий научным отделом, e-mail: s.v.balalin@amail.com ID: 0000-0002-5250-3692

Эта статья доступна в соответствии с условиями лицензии Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0)

This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0)

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

©ПНИПУ

Введение

Математическое моделирование является эффективным инструментом для изучения кератотопогра-фии, биомеханики и физиологии роговицы (и глаза в целом) в норме и при различной патологии. При этом современный уровень разработки таких моделей, подразумевающий по сути создание цифрового «двойника» исследуемого объекта, позволяет не только существенно дополнить сложные, дорогостоящие экспериментально-клинические исследования, но и расширить возможности изучения ряда физиологических процессов, недоступных для наблюдения in vivo (а часто и in vitro) с использованием существующих технических средств, а также применять компьютерные модели в качестве виртуального тренажера в учебном процессе или при отработке различных стратегий лечения заболевания [5; 12].

В настоящее время известны работы как в области построения полных биомеханических моделей глаза [5; 12], так и моделей его отдельных анатомо-оптических структур (хрусталика и цинновых связок, склеры, роговицы, глазодвигательных мышц и др.), предназначенных для решения широкого круга офтальмологических задач [1; 5; 12; 35]. Для их реализации используются различные методы - аналитические [1; 2], статистические [6; 13; 34], численные (метод «стрельбы» [4; 7-9; 14; 42], метод конечных разностей [4], метод конечных элементов [41; 42; 45; 48; 49]), при этом наибольшее распространение в последнее время получает метод конечных элементов (МКЭ), в связи с его гибкостью при описании особенностей геометрии и свойств элементов глаза при минимальном наборе исходных допущений [42].

Компьютерные модели роговицы глаза сегодня с успехом применяются при изучении, диагностике и прогнозировании результатов лечения при миопии, гиперметропии, астигматизме, кератоконусе [3 ; 41; 45; 48; 49]. Областями приложения такого рода моделей являются:

• задачи описания топографических и томографических показателей роговицы во взаимосвязи с особенностями распределения характеристик ее жесткости в норме и патологии [10; 17; 23], в том числе при

моделировании процедур лечения заболеваний роговицы (например, при имплантации интрастромальных колец [10; 32], кератопластике [22] или кросслинкинге роговичного коллагена [50; 51]);

• задачи моделирования механического поведения и напряженно-деформированного состояния роговицы при контактной [6-9; 13; 14] и бесконтактной тонометрии [16; 30; 37; 40; 43; 46; 53; 54; 57] в зависимости от внутриглазного давления [6; 8; 13], особенностей строения (неоднородность, анизотропия) [7; 9; 14; 16; 54; 56; 57] и физических моделей деформирования материала роговицы (нелинейно-упругий, гиперупругий, вязкогиперупругий материал и др.) [22; 30], включая анализ приложений полученных результатов при диагностике и лечении [9; 22]. При построении комплекса моделей немаловажную роль играют современные экспериментальные исследования геометрии, структуры и свойств материала роговицы (и других структур глаза) как на основе передового клинического оборудования, так и уникальных методик испытания, позволяющих восполнить ограниченный объем данных об особенностях напряженно-деформированного состояния, механического поведения и разрушения трудно исследуемого материала роговицы (в том числе у здоровых лиц и пациентов с различной патологией роговицы до и после лечения) [19; 31; 33; 47; 55].

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

К особенностям разрабатываемого подхода, отличающих его от известных работ [30; 43; 50; 56], рассматривающих идеализированные геометрические модели роговицы или усредненные для группы пациентов ее свойства, можно отнести использование пациенто-индивидуализированных методик определе-

ния геометрических и биомеханических параметров роговицы. Они устанавливаются непосредственно из клинических данных конкретного пациента, либо, при отсутствии или невозможности прямых замеров, - по результатам серии вычислительных экспериментов, воспроизводящих соответствующие диагностические процедуры.

Исходные данные. Геометрическая модель роговицы

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

Построение 3^-модели роговицы (с использованием компьютерной системы конечно-элементного моделирования Comsol Multiphysics) заключалось в линейной интерполяции экспериментально найденного массива точек (7000-10 000 точек) ограничивающими поверхностями (передней и задней), «натянутыми» на указанный массив высот (координат), как на каркас (рис. 2). При таком подходе исключается промежуточная операция аппроксимации экспериментальных данных аналитическими зависимостями (например, часто для этих целей применяются полиномы Цернике [10; 50], уравнения эллипсоида [43] или сферы [3; 35]), что позволяет «индивидуализировать» исследуемую модель роговицы, наиболее полно описать особенности ее геометрии, характерные для возможных случаев нормы или патологии, произвести «тонкую» настройку параметров математической модели на контрольных группах пациентов. При этом за пределами области роговицы, доступной для измерения при помощи Pentacam, расчетная поверхность задавалась уравнением некото-

а б

Рис. 1. Поля координат передней (а) и задней (б) поверхностей роговицы здорового пациента, построенные по данным керато-топографии Pentacam AXL

рой трехмерной аналитической функции, автоматически подбираемой при помощи программы TableCurve 3D, параметры которой устанавливались из аппроксимации экспериментальных данных.

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

плоскостях а Г (рис. 2, в). Окончательный вид 3D-модели роговицы приведен на рис. 2, г.

К исходным данным относились полученные на этой стадии величины: максимальный диаметр роговицы атах, толщины роговицы в апексе tapex, глубина передней камеры Нс. Фиксировались величины минимальной толщины роговицы , максимальной кривизны Ктах передней и задней поверхностей роговицы и координаты соответствующих точек (Хр , Ур, хк , ук )

(табл.1). Остальные геометрические и топографические параметры роговицы (определяемые картами толщины, тангенциальной и сагиттальной кривизны) считались производными от координат и использовались для верификации расчетной геометрии построенной 3D-модели. Также задавалась величина внутриглазного давления р^ар по данным тонометрии пациента (с использованием тонометра Corvis ST).

Значения указанных параметров применительно к роговицам левого и правого глаза конкретного (анонимного) пациента приведены в табл. 1. Соответствующая твердотельная модель роговицы с указанием ее основных размеров показана на рис. 3 (при этом выделены области поверхности, построенные по данным Pentacam, и области, аппроксимированные аналитическим уравнением).

Механическая модель роговицы. Основные соотношения

Основные уравнения механики деформирования роговицы (в формулировках и обозначениях, принятых в системе конечно-элементного моделирования Comsol Multiphysics) записываются в следующем виде [18; 28].

Уравнения квазистатического равновесия:

V( )Т + Яу = 0, (1)

где V - оператор набла; T - символ транспонирования; Яу - вектор объемных сил; S - второй тензор напряжений Пиолы - Кирхгофа; F - градиент деформации, определяемый как

Д = I + V, (2)

где I - единичный тензор второго ранга, и - вектор перемещений.

а б в г

Рис. 2. Интерполяция экспериментального массива точек для передней (а) и задней (б) поверхности роговицы, формирование сплошного тела (в), конечная форма 3.0-модели роговицы (г)

Таблица 1

Исходные экспериментальные данные для построения модели

Глаз ,max dh jmax dv hc, 'apex 'min xp Ур Kmax xK , Ук , piop

мм мм мм мм мм мм мм дптр мм мм Па

Правый 12,3 11,76 3,09 0,521 0,512 -0,68 -0,96 46,1 -0,68 -1,84 2053,2

Левый 12,32 11,66 3,1 0,519 0,507 0,604 0,82 48,1 -0,68 -2,19 3199,7

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

д Ж

S — So.

ÖE

(3)

где Бех{ - тензор внешних (дополнительных) напряжений; - тензор деформации Грина - Лагранжа

* = 1 (FTF - /).

Тензор напряжений Коши в данной формулировке запишется в виде

о = J~1Е8ЕТ , (4)

где J = det (F) - определитель градиента деформации

¥, характеризующий изменение объема, вызванное деформацией. Данные уравнения дополняются граничными условиями, заданными на соответствующих поверхностях Г тела (рис. 4):

1) отсутствие перемещений на поверхностях роговицы, примыкающей к склере Г £

и = 0, ГеГ^ ; (5)

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

которой конфигурации роговицы при решении обратной задачи

£• п = чъ, ГеГь; 3) действие на передней поверхности Г^ роговицы

поверхностно распределенной нагрузки интенсивностью Ч[ с параметрами, соответствующими параметрам воздушного импульса при бесконтактной тонометрии

£• п = ч/, Г еГу. (7)

Вид функций для чъ и Ч[ в условиях (6), (7) будет уточнен далее.

г.

аппроксимация функцией

а б

Рис. 3. Общий вид и основные параметры геометрической 3Б-модели роговицы: а - вид спереди; б - вид сбоку

Рис. 4. Расчетная схема модели роговицы: а - начальное состояние; б - деформированное состояние

Модель материала роговицы

Физические соотношения для роговицы устанавливаются в соответствии с принятой в данном исследовании обобщенной моделью вязкогиперупругого анизотропного материала, в рамках которой функция плотности энергии деформации Щ может быть записана в виде [52]

Щ = + Кы + ^ ъ + ^ . (8)

Составляющие функции Щ в выражении (8) интерпретируются следующим образом.

Величина представляет собой изотропную

часть удельной энергии упругой деформации при постоянном объеме без учета анизотропии деформации вдоль характерных направлений коллагеновых волокон роговицы. Для ее нахождения может использоваться, например, эмпирическая модель Муни - Ривлина (Mooney - Rivlin) для почти несжимаемого гиперупругого материала [50]

Що = С10 (II - 3) + с20 (/1 - 3)2, (9)

где Сю, С20 - эмпирические параметры модели; 11 -первый инвариант (след) тензора С, а именно

д = 1г (С), (10)

С = J"2/3 С; (11)

где C - правый тензор деформаций Коши - Грина.

Составляющая Щуо1 соответствует энергии объемной упругой деформации и может рассматриваться как «штрафной» член для обеспечения требования несжимаемости [45]

1 ? = ^ к( J-1)2,

(12)

где К - модуль объемной упругости.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Слагаемое задает функцию плотности энергии

деформации коллагеновых волокон и описывает анизотропное поведение материала роговицы, включая влияние структуры и расположения характерных групп фибрилл [44; 54]

ЩАЬ = X

к

=л е, 2к->

ехр

къ (/;- з)

(13)

где к14, к16, к2

1*=К 11 +(1 - 3кг) I ,, (14)

к24 , к26 - параметры жесткости анизотропного материала; к■ е [0,13] - параметр структуры (ориентации) волокон, который может быть различным для каждого семейства фибрилл; 14, 16 - четвертый и

шестой псевдоинварианты модифицированного тензора деформаций С из (11), определяемые как [44; 45]

14 = \ • С \, 76 = £• С С, (15)

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

Для материала роговицы с двумя характерными группами волокон в (13), (14) можно принять [52] &!4 = ^ = и ¿24 = ¿26 = ¿2, при этом ки трактуются [54] как коэффициенты упругости волокон при простом растяжении и имеют размерность напряжений, ^г -как безразмерные параметры жесткости фибрилл при больших деформациях. Параметр к 1 в выражении (14) соответствует степени разориентации волокон вдоль рассматриваемых направлений: кг = 0 - 100 % колла-геновых фибрилл ориентированы в заданном направлении; кг = 1/3 - изотропное (неориентированное) распределение фибрилл [44].

Согласно [28; 29; 57], вязкостные свойства материала роговицы задаются в (8) «диссипативной» составляющей

ж = У у

угх .У а а=1

(16)

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

В таком случае вторые тензоры напряжений Пиолы - Кирхгофа дополняются при помощи вспомогательных компонент неравновесных напряжений Qа [18; 28; 29; 57]

дУ

Qа = 2-

д С

(17)

Динамика изменения тензора вспомогательных напряжений Qa в каждом структурном элементе модели может быть описана уравнением [18; 28; 29; 57]

1

В уравнении (18)

6а 4 6а ~ Ра^гзо ■ Та

д Ж = 2 ■

д С

(18)

(19)

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

На каждом из уровней сложности модели материала составляющие уравнения (8) вводятся в его структуру по мере необходимости описания особенностей механического поведения роговицы при сопоставлении расчетных и экспериментальных данных, а входящие в выражения (9), (12)-(14), (18) параметры (с10, с20 , К ,

¿1, ¿2, к г, Ра, та) уточняются из серии вычислительных экспериментов. Предварительные значения указанных коэффициентов принимаются по литературным данным (например, [15; 44; 50; 52; 54; 57]) и приведены в табл. 2.

Модель кератоконуса

В первом приближении нормальная роговица рассматривалась как изотропный гиперупругий материал, закон деформирования которого задан соотношениями (8)-(12). Модель анизотропного материала, учитывающая распределение коллагеновых фибрилл в ткани роговицы (ортогональное распределение в центральной части, периферическое - около лимба), как и модель вязкоупругости, при необходимости вводится на следующих уровнях сложности для уточнения особенностей механического поведения материала и будет рассмотрена дополнительно. При этом принятое в данной работе начальное допущение о пространственной изотропии представляется оправданным для роговицы с кератоконусом, где предпочтительная ориентация кол-лагеновых фибрилл и типичное распределение их плотности нарушены и имеют более случайный характер, особенно в области конуса [39; 50].

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

В этой зоне задается закон распределения эмпирических коэффициентов Сю, С20 , устанавливающий их снижение от периферии к центру как по поверхности (по радиусу), так и по толщине роговицы при помощи следующей функции Уk (в локальной системе координат хк , ук , , связанной с центром зоны кератоконуса - точкой Ко на рис. 5; ось хк - по нормали к задней поверхности роговицы, ось ук - параллельна назально-темпоральной плоскости, ось - параллельна нижневерхней плоскости):

Vk (хк, Ук, 2к ) = (1 - Ушах ) ехР[0^ хк /4Х +

+ekУk2/+ 4 4 «2 ]

(20)

У ?

где Ушах - максимальное относительное снижение жесткости в центре зоны кератоконуса; 0k , 0k , 0k -

параметры, устанавливающие градиент изменения свойств по рассматриваемой области (вдоль соответствующей оси); ak , % , % - параметры, устанавли-

вающие размер (длину полуоси) зоны кератоконуса вдоль каждой из осей.

С учетом выражения (20) в области пониженной жесткости выражение (9) для упругого потенциала можно записать так [50]:

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

плоскость)

Таблица 2

Значения параметров механической модели материала роговицы

cio, кПа c2o , кПа К , МПа k1, МПа ¿2 кг Pi Ti , с

140,9 80,13 5,5 0,04 200 0,3329 0,351 410

щКо (хк.Ук. 2к) = У к (хк >Ук > 2к) + [с10 (11 - 3) +

\2 , (21)

+ С20 (/1 - 3 ) ]

Отметим, что количество таких областей и их расположение может быть произвольным, за счет чего появляется возможность описания достаточно сложных конфигураций зоны кератоконуса (например, на рис. 5 использовано наложение двух областей с различной степенью снижения свойств). Центр зоны кератоконуса (точка K0) устанавливается по результатам клинического обследования пациента с помощью Шаймпфлюг-анализатора переднего отрезка глазного яблока Pentacam AXL.

Модель кросслинкинга

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

Современные источники УФ-излучения для лечения кератоконуса обеспечивают различные закономерности распределения интенсивности луча, например [27; 50]: равномерное, со снижением интенсивности вдоль радиуса пятна облучения около 10 % (рис. 6, а; кривая 1); с

повышенной интенсивностью на периферии (рис. 6, а; кривая 2). С другой стороны, развиваются методики кросслинкинга, предполагающие плавное снижение интенсивности УФ-излучения от максимума в центре до нуля на границе зоны облучения (рис. 6, а; кривая 3) для предотвращения возможного резкого скачка свойств на границе облученного и необлученного материала. Снижение интенсивности УФ-излучения происходит также и по глубине роговицы по закону, близкому к экспоненциальной зависимости Бира - Ламберта [50].

Распределение концентрации рибофлавина в объеме роговицы определяется диффузионным уравнением 2-го закона Фика и в первом приближении также может быть задано некоторой экспоненциальной зависимостью, что согласуется с известными расчетными и экспериментальными данными (рис. 6, б) [38].

Таким образом, увеличение модулей упругости, вызванное кросслинкингом роговичного коллагена, в любой точке роговицы (с координатами xc, yc, zc, отсчитываемыми от центра зоны облучения, рис. 7, точка СО) является функцией радиуса и глубины, которую, по аналогии с функцией (20), примем в виде:

Yc (*c,Ус> Zc ) = (1 + Ymax )exp[9c x2c Ia2c +

2/2 * 2 / 2 ' (22) + Qcry2c /a\ + 0cz2/a2 ]

а б

Рис. 6. Экспериментальные данные [27; 38; 50] о характерных распределениях интенсивности ф0 ультрафиолетового излучения (а) и диффузии % 0 рибофлавина при различном времени выдержки по диаметру ё0 зоны воздействия (б)

(все параметры нормализованы по максимальным значениям)

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

(б, г) поверхностям (а, б - трехмерный вид; в, г - проекции на плоскость)

где ^шах - максимальное относительное увеличение параметров жесткости ( С10 , С20 ) в результате кросс-линкинга; 0С , 9С , 0С - параметры, устанавливающие

градиент изменения свойств по рассматриваемой области (вдоль соответствующей оси); аС , аС , аС - пара-

Сх у 1

метры, устанавливающие размеры зоны кросслинкинга вдоль каждой из осей (для области облучения круглой формы аС = аС , а также 0С = 0С ).

Полная функция, устанавливающая величину упругой деформации при постоянном объеме в модели Муни - Ривлина материала роговицы (9) с учетом деградации

свойств в зоне кератоконуса (20) и повышения параметров жесткости материала после операции кросслинкинга роговичного коллагена, (22) примет вид:

с10 (Ä - з) + С20 (Ä - з)2

Wuo = W

(23)

На рис. 7 иллюстрируется возможное расчетное распределение параметра жесткости роговицы с керато-конусом (заданной на рис. 5) после повышения модулей упругости операцией кросслинкинга в локальной области с центром в точке Со, координаты ( yc , zc ) которой назначаются в зависимости от выбранной методики кросслинкинга (например, соответствуют центру кера-токонуса).

и = О

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

и под действием ВГД

Модель начальной геометрии роговицы

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

Известны два основных подхода к решению таких задач [20]. В первом, который используется в данной работе, к задней поверхности роговицы прикладывается давление p¡nv, обратное ВГД по направлению (знаку) [24; 25; 36; 43]. В таком случае компоненты поверхностной нагрузки qb в граничном условии (6) для задней поверхности роговицы будут записаны в виде (рис. 8, а):

qb =-Pmv ■ n ГеГь, (24)

при этом величина обратного давления p nv устанавли-

вается равной ВГД, скорректированному на разницу площадей задней поверхности роговицы в деформированном и недеформированном состоянии

Ргш = -ЪРор , (25)

где р0р - величина ВГД, установленная по данным

Corvis ST; X - коэффициент, учитывающий изменение площади поверхности приложения прямого р ор и обратного давлений рш в деформированном и недефор-мированном состоянии роговицы.

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

qь =-рЮр • п ГеГь. (26)

Отметим, что в качестве одного из ограничений этого подхода авторы [20; 26] указывают возможность проявления бифуркации («прощелкивания») тонкостенной гибкой оболочки (которую представляет собой роговица) при приложении обратного давления согласно (6). Однако в данной работе подобные эффекты не наблюдались.

МС

а б

Рис. 9. Картины распределения давления воздушной струи при бесконтактной тонометрии по диаметру ё зоны воздействия (а) и изменения максимального давления по времени t импульса (б). На графике (а) сплошными линиями показаны различные моменты времени по данным [16; 37; 40], точками - расчет по формуле (29); на графике (б) точками - профиль давления по данным Свп>18 5Т

Другой подход, используемый, например, в [20; 46], основан на расчете деформаций роговицы с помощью МКЭ и постепенном итерационном уточнении перемещений узлов сетки конечных элементов к конфигурации без напряжений. К ограничениям такого подхода можно отнести его сравнительную трудоемкость и необходимость определения и сохранения на каждой итерации трехмерных полей перемещений и координат точек роговицы с определением степени отклонения текущего состояния от экспериментально найденной формы роговицы.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Модель внешних нагрузок при пневмотесте

Для идентификации параметров материала в соотношениях (8)-(12), (20)-(23) разрабатывается модель бесконтактной тонометрии, воспроизводящая процесс клинико-экспериментального исследования биомеханических параметров роговицы пациентов при помощи пневматического тонометра Свп>18 5Т.

При моделировании на передней поверхности роговицы задавалось граничное условие (7), устанавливающее действие поверхностных нагрузок ч/, соответствующих давлению рагг дозированного коллимиро-ванного воздушного импульса

Чг =-Рагг • П ГеГГ . Пространственно-временная конфигурация давления раГ (х, у, 7, t) определяется закономерностями взаимодействия воздушного потока, создаваемого импульсом пневмотонометра, с поверхностью роговицы. Распределение раГ (х, у, 7, t) представляется в виде сочетания

функции координат поверхности (у, 7) (так как

функция fs задается по передней поверхности Гf роговицы, координата x ее точек предопределена и может быть опущена) и функции времени ft (t):

Par (^ У z=t) = Pmax " fs (У z)" ft (t) , x еГ f • (28)

Здесь pmax - максимальное давление в центре воздушного импульса, величина которого, так же, как вид и параметры функций fs (y, z) и ft (t), устанавливались по экспериментальным данным (замеры Corvis ST) и результатам работ [16; 30; 37; 40; 43; 46; 53; 54; 57], посвященных моделированию динамики взаимодействия воздушной струи и роговицы глаза (рис. 9, а).

В частности, по аналогии с картинами расчетного распределения давления воздушной струи по поверхности роговицы, полученными авторами [16; 37] по результатам компьютерного моделирования бесконтактной тонометрии (рис. 9, а), функция пространственного распределения давления fs (y, z) задается по круговой области (диаметр которой определяется диаметром dair сопла прибора и коэффициентом «сосредоточенности» 9 p ; принято dair = 2,4 мм [37] и 0p = -0,75) на передней поверхности роговицы в виде

fs (y, z) = exp 0p (y2 + z2)/(dair /2)2

(29)

График функции fs (y, z) , построенный в соответствии с (29), показан на рис. 9, а.

Временной профиль давления струи приборно фиксируется в ходе пневмотонометрического теста глаза пациента на базе Corvis ST и сохраняется в табличный файл (параметр Pressure Profile) для 140 точек исследуемого интервала времени (рис. 9, б). Полученные экспериментальные результаты используются для формиро-

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

На рис. 10 показаны картины распределения нагрузок, действующих на передней и задней поверхностях роговицы в различные моменты времени в ходе бесконтактной тонометрии воздушным импульсом: на задней поверхности - внутриглазное давление согласно (6), (26); на передней поверхности - переменное давление струи согласно (7) с учетом (27)-(29).

Расчетные кератотопографические карты роговицы

На основе разработанной модели (1)-(29) с использованием компьютерной системы конечно-элементного анализа Comsol Multiphysics производилось моделирова-

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

Результаты сопоставительного анализа иллюстрируются на примере пациента, один глаз которого не имел значимых кератэктатических отклонений (считался здоровым), на другом - был диагностирован прогрессирующий кератоконус I степени. Для этого глаза пациенту было назначено хирургическое лечение с выполнением операции кросслинкинга роговичного коллагена по модифицированной методике, разработанной авторами [11].

На рис. 11 показаны карты толщин роговицы указанного пациента, полученные по результатам измерений при помощи Pentacam AXL и расчетным путем на основе разработанной модели. Приводятся экспериментальные данные и результаты математического моделирования для нормальной роговицы (рис. 11, а, б), роговицы с кератоконусом (рис. 11, в, г) и той же роговицы с кератокон

б в г д

Рис. 10. Картины расчетного распределения давления воздушной струи по передней поверхности и внутриглазного давления по задней поверхности роговицы в различные моменты времени от начала импульса (а - 6 мкс, б - 9 мкс, в - 12 мкс, г - 16 мкс, д - 19 мкс, е - 22 мкс)

Таблица 3

Сопоставление экспериментальных (Pentacam) и расчетных (модель) данных о толщине роговицы (мкм) в характерных точках

а

е

№ точки

Параметр 1 2 3 4 5 6

Нормальная роговица

Pentacam 509 520 517 546 551 657

Модель 508,7 517,5 515,9 533,5 535,7 638,8

Отклонение, % 0,06 0,48 0,21 2,28 2,77 2,77

Кератоконус

Pentacam 507 519 519 558 537 671

Модель 508,7 518,7 518,7 554,1 537,1 656,5

Отклонение, % 0,34 0,06 0,06 0,7 0,02 2,16

Кросслинкинг

Pentacam 502 514 514 548 531 663

Модель 501,02 512,7 513,5 546 527,6 651,2

Отклонение, % 0,2 0,25 0,1 0,36 0,64 1,78

усом через две недели после операции кросслинкинга (рис. 11, д, е).

Для всех рассмотренных состояний роговицы (нормальная; с кератоконусом; с кератоконусом после кросс-линкинга) достигнуто удовлетворительное качественное (см. рис. 11) и количественное (табл. 3) соответствие опытных и расчетных данных - максимальные взаимные отклонения результатов не превышают 2-3 %.

Для этого производился пересчет главных кривизн Р1, Р2 модельных поверхностей, определяемых встроенными средствами Сот^о1 Multiphysics, в величины кривизны данных поверхностей в меридиональной (тангенциальной) р^ (рис. 12, 13) и сагиттальной р

плоскостях.

Соответствующие формулы пересчета записывались так [21]:

{ \ 2 (у • Пу + г • п2 )

Р^ =(Р1 - Р2 )•

)

+р2; (30)

у •

Р sag

2 , 2 у + г •

22 у + г +

Л 2

у •

где Пх

Пу, пг - компоненты вектора п нормали в дан-

ной точке поверхности.

В качестве примера на рис. 12, 13 приведены экспериментальные (по данным РвМасат ЛХЬ) и расчетные (по результатам моделирования) карты тангенциального радиуса кривизны передней (см. рис. 12) и задней (см. рис. 13) поверхностей роговицы в различных состояниях (для пациента, указанного ранее): здоровой роговицы (см. рис. 12, 13; а, б), роговицы с кератоконусом (см. рис. 12, 13; в, г) и роговицы с кератоконусом после операции кросслинкинга (см. рис.12, 13; д, е). Количественное сопоставление соответствующих данных приведено в табл. 4 (для тех же характерных точек, что использовались в табл. 3).

п

Рис. 11. Экспериментальные (а, в, д) и расчетные (б, г, е) карты распределения толщины роговицы: а, б - нормальной; в, г - с кератоконусом; д, е - после кросслинкинга

Рис. 12. Экспериментальные (а, в, д) и расчетные (б, г, е) карты распределения тангенциального радиуса кривизны передней поверхности роговицы: а, б - нормальной; в, г - с кератоконусом; д, е - после кросслинкинга

Рис. 13. Экспериментальные (а, в, д) и расчетные (б, г, е) карты распределения тангенциального радиуса кривизны задней поверхности роговицы: а, б - нормальной; в, г - с кератоконусом; д, е - после кросслинкинга

Рис. 14. Экспериментальные и расчетные картины деформаций роговицы при действии воздушного импульса (а - нормальная роговица, максимальный прогиб; б - роговица после кросслинкинга, максимальный прогиб)

Таблица 4

Сопоставление экспериментальных (РвМасаш) и расчетных (модель) данных о радиусах тангенциальной кривизны роговицы (мм) в характерных точках

№ точки

Параметр 1 2 3 4 5 6

Передняя поверхность

Нормальная роговица

РвШасат 7,73 7,65 7,73 7,36 7,93 9,72

Модель 7,52 7,53 7,69 7,36 7,82 8,33

Отклонение, % 2,72 1,57 0,51 0 1,39 14,3

Кератоконус

РвШасат 7,72 7,56 7,6 7,09 8,03 9,09

Модель 8,06 7,69 7,83 7,29 8,62 9,59

Отклонение, % 4,4 1,72 3,03 2,82 7,34 5,5

Кросслинкинг

РвШасат 7,68 7,5 7,53 7,03 8,05 9,06

Модель 7,53 7,26 7,47 7,27 8,57 9,53

Отклонение, % 1,95 3,2 0,8 3,41 6,46 5,19

Задняя поверхность

Нормальная роговица

РвШасат 6,13 6,57 6,78 5,73 6,41 8,96

Модель 5,95 7,02 7,21 5,52 6,07 11,79

Отклонение, % 2,94 6,85 6,34 3,66 6,08 31,6

Кератоконус

РвШасат 6,21 6,62 6,45 5,91 6,68 8,3

Модель 6,26 6,31 6,59 5,94 7,09 8,12

Отклонение, % 0,81 4,68 2,17 0,51 6,13 2,17

Кросслинкинг

РвШасат 6,11 6,61 6,39 5,64 6,73 8,7

Модель 5,95 6,27 6,41 5,81 7,29 9,53

Отклонение, % 2,62 5,14 0,31 3,01 8,32 9,54

Таблица 5

Сопоставление экспериментальных (Со/ги) и расчетных (модель) данных о параметрах деформации роговицы в характерных точках

Роговица

Параметр нормальная после кросслинкинга

БЛ, мм РБ, мм БЛ, мм РБ, мм

Corvis 1,09 4,93 0,87 4,14

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Модель 1,06 4,9 0,879 4,347

Отклонение, % 2,75 0,6 1,03 5

Приведенные данные указывают на достигнутое качественное соответствие результатов математического моделирования основных топографических параметров роговицы их экспериментальным распределениям: наибольшие отклонения расчетных величин наблюдаются в периферийных точках - от 7-10 до 15-30 % в некоторых случаях, тогда как в центральной зоне роговицы погрешности не превышают 4-7 %.

Результаты моделирования бесконтактной тонометрии

При моделировании пневмотонометрического теста сопоставлялись расчетные профили роговицы в различные моменты процесса ее деформирования импульсом воздушной струи (на основе (27)-(29)) с соответствующими динамическими картинами, полученными (для сечения роговицы в назальной плоскости) при помощи Шаймпфлюг-камеры бесконтактного тонометра Corvis 8Г (рис.14).

Основными контролируемыми параметрами являлись амплитуда деформации роговицы (ДЛ), скорость перемещения в апексе (СУ), расстояние между пиками (РД), изменение длины дуги (*МЕ), а также первая и вторая апланации при нагрузке и разгрузке давлением воздушной струи.

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

Численные результаты такого сопоставления в момент времени, соответствующий максимальной деформации роговицы (17,56 мс - для здоровой роговицы; 17,79 мс - для роговицы после кросслинкинга), приведены в табл. 5. Для указанного момента наибольшие отклонения составляют не более 3-5 %. Однако процесс деформации роговицы под действием воздушной струи характеризуется достаточно высокой степенью неустойчивости, возникновением колебаний и отклонением траектории апекса от прямолинейной, что в отдельные моменты может сопровождаться несовпадением расчетных и экспериментальных точек на 25-30 %. В целом же принятая стратегия моделирования показала свою адекватность и работоспособность.

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

Заключение

По результатам проделанной работы можно сделать следующие заключения:

1. На базе системы конечно-элементного анализа Сот^о1 Multiphysics разработан комплекс математических моделей роговицы, позволяющий производить численный анализ ее геометрических и томографических параметров, биомеханических свойств и напряженно-деформированного состояния в условиях действия внутриглазного давления и внешних воздействий при диагностике и лечении. Рассмотрены основные этапы построения такого комплекса, его базовые модули, определяющие соотношения и результаты сопоставительного анализа.

2. Разработана геометрическая модель роговицы в виде пространственного сегмента выпуклой тонкостенной оболочки с переменной толщиной стенки и произвольной формой передней и задней поверхностей. Учет реальной геометрии роговицы конкретного пациента производился путем прямой (поточечной) интерполяции экспериментальных данных (карт высот), полученных с помощью кератотопографа РвШасат ЛХЬ. Для отыскания недеформированной поверхности роговицы применялась методика решения обратной задачи с приложением давления, противоположного внутриглазного давления и обнулением полученных напряжений (начальное ненапряженное состояние).

3. В первом приближении использована физическая модель гиперупругого материала в форме Муни - Ривли-на, идентификация параметров которой производилась на основе сопотавительного анализа результатов численных экспериментов о деформировании исследуемой модели роговицы и опытных данных по бесконтактной тонометрии глаза пациента при помощи пневмотонометра Corvis 5Г.

4. Разработаны модели кератоконуса и кросслинкин-га роговичного коллагена, представляющие собой функциональные зависимости, устанавливающие закон рас-

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

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

6. Сопоставительный анализ результатов математического моделирования и полученных экспериментальных данных показал удовлетворительное их согласование. Относительные отклонения при моделировании большинства топографических и биомеханических параметров роговицы в различных состояниях не превышали 5-10 %.

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

1. Бауэр С.М. Модели теории оболочек и пластин в офтальмологии: автореф. дис. д-ра физ. мат. наук. -СПб., 2002. - 22 с.

2. Бауэр С.М., Замураев Л.А., Котляр К.Е. Модель трансверсально-изотропного сферического слоя для расчета изменения внутриглазного давления при интрасклеральных инъекциях // Российский журнал биомеханики. - 2006. - Т. 10, № 2. - С. 43-49.

3. Бауэр С.М., Венатовская Л.А., Франус Д.В., Федотова Л.А. Оценка изменения напряженно-деформированного состояния глаза и показателей внутриглазного давления после рефракционной коррекции гиперметропии // Российский журнал биомеханики. - 2015. - Т. 19, № 2. - С. 136-143.

4. Ермаков А.М. Напряженно-деформированное состояние склеры и роговицы как ортотропных неоднородных сопряженных сферических оболочек // Российский журнал биомеханики. - 2009. - Т. 13, № 1.

- С. 49-60.

5. Иомдина, Е.Н., Полоз М.В. Биомеханическая модель глазного яблока человека: описание и верификация // Математическая биология и биоинформатика. - 2014.

- Т. 9, № 1. - С. 286-295.

6. Качанов А.Б. , Бауэр С.М., Воронкова Е.Б. и др. Статистическая оценка влияния некоторых параметров глазного яблока на тонометрическое внутриглазное давление // Российский журнал биомеханики. - 2018. - Т. 22, № 4. - С. 527-536.

7. Любимов Г.А., Моисеева И.Н., Штейн А.А. Исследование свойств двухкомпонентной механической модели глазного яблока и возможности ее использования при практической оценке механических свойств глаза человека // Механика жидкости и газа. - 2014. - № 6. - С. 5-16.

8. Моисеева И.Н., Штейн А.А. Анализ зависимости «давление-объем» для глазного яблока, нагруженного плоским штампом, на основе двухсегментной упругой модели // Механика жидкости и газа. - 2011. - № 5. -С. 3-15.

9. Моисеева И.Н., Штейн А.А. Математическое моделирование эластотонометрии по Маклакову в случае искусственно созданной неоднородности роговицы // Российский журнал биомеханики. - 2019.

- Т. 23, № 1. - С. 8-21.

10. Никитин И.С., Журавлев А.Б., Ирошников Н.Г. и др. Механико-математическая модель интрастромальной коррекции формы роговицы глаза при кератоконусе // Российский журнал биомеханики. - 2017. - Т. 21, № 4.

- С. 404-417.

11. Пат. 2760482 Российская Федерация, МПК A61F 9/007, A61F 9/008. Способ лечения прогрессирующего кератоконуса / Е.Г.Солодкова, В.П. Фокин, Е.В. Лобанов; заявитель и патентообладатель ФГАУ «МНТК "Микрохирургия глаза" им. акад. С. Н. Федорова» Минздрава России. - № 2021104215/14 ; заявл. 19.02.2021; опубл. 25.11.2021, Бюл. № 33. - 7 с.

12. Полоз М.В. Биомеханическая модель глазного яблока человека: автореф. дис. физ.-мат. наук. - Саратов, 2011. - 20 с.

13. Слесорайтите Е. Статистический и численный анализ влияния толщины роговицы на показатели внутриглазного давления // Российский журнал биомеханики. - 2006. - Т. 10, № 2. - С. 58-63.

14. Штейн А.А., Моисеева И.Н., Любимов Г.А. Математическая модель роговицы глаза с учетом экспоненциальной нелинейности ее упругих свойств при условии геометрической малости деформаций // Российский журнал биомеханики. - 2019. - Т. 23, № 3.

- С. 375-390.

15. Abyaneh M.H., Wildman R.D., Ashcroft I.A., Ruiz P.D. A hybrid approach to determining cornea mechanical properties in vivo using a combination of nano-indentation and inverse finite element analysis // Journal of the Mechanical Behavior of Biomedical Materials. - 2013. -Vol. 27. - P. 239-248.

16. Ariza-Gracia M. A., Zurita J.F., Pinero D.P. et al. Coupled Biomechanical Response of the Cornea Assessed by Non-Contact Tonometry. A Simulation Study // PLoS ONE. -2015. - No. 10 (3) : e0121486. - DOI: 10.1371/journal.pone.0121486.

17. Carvalho L.A., Prado M., Cunha R.H. et al. Keratoconus prediction using a finite element model of the cornea with local biomechanical properties // Arquivos Brasileiros de Oftalmologia. - 2009. - No. 72(2). - P. 139-145.

18. COMSOL Documentation. Structural Mechanics Module

- Large Strain Viscoelasticity. [Электронный ресурс] -URL: https://doc.comsol.com/5.6/docserver/#!/com. comsol.help.sme/sme_ug_theory.06.27.html (дата обращения 04.07.2022).

19. Dias J., Diakonis V.F., Kankariya V.P. et al. Anterior and posterior corneal stroma elasticity after corneal collagen crosslinking treatment // Experimental Eye Research. -2013. - Vol. 116. - P. 438-451.

20. Elsheikh A., Whitford C., Hamarashid R. Stress free configuration of the human eye // Medical Engineering & Physics. - 2013. - Vol. 35. - P. 211-216.

21. Estrada-Molina A., Campos-García M., Díaz-Uribe R. Sagittal and meridional radii of curvature for a surface with symmetry of revolution by using a null-screen testing method // Applied Optics. - 2013. - Vol. 52, No. 4. - P. 625-634.

22. Fraldi M., Cutolo A., Esposito L., Guarracino F. The role of viscoelasticity and stress gradients on the outcome of conductive keratoplasty // Biomechanics and Modeling in Mechanobiology. - 2011. - Vol. 10 (3). - P. 397-412.

23. Gefen A., Shalom R., Elad D., Mandel Y. Biomechanical analysis of the keratoconic cornea // Journal of the Mechanical Behavior of Biomedical Materials. - 2009. -Vol. 2 (3). - P. 224-236.

24. Govindjee S, Mihalic P.A. Computational methods for inverse finite elastostatics // Computer Methods in Applied Mechanics and Engineering. - 1996. - Vol.136 (1-2). - P. 47-57.

25. Govindjee S., Mihalic P.A. Computational methods for inverse deformations in quasi-incompressible finite elasticity // International Journal for Numerical Methods in Engineering. - 1998. - Vol. 43 (5). - P. 821-841.

26. Grytz R., Downs J.C. A forward incremental prestressing method with application to inverse parameter estimations and eye-specific simulations of posterior scleral shells // Computer Methods in Biomechanics and Biomedical Engineering. - 2013. -Vol. 16, No. 7. - P. 768-780.

27. Herber R., Kunert K.S., Veliká V. et al. Influence of the beam profile crosslinking setting on changes in corneal topography and tomography in progressive keratoconus: Preliminary results // Journal of Cataract & Refractive Surgery. - 2018. - Vol. 44, No. 6. - P.718-724.

28. Holzapfel G.A. Nonlinear solid mechanics : A continuum approach for engineering. - Chichester [etc.]: Wiley, 2000. - 455 p.

29. Holzapfel G.A., Gasser T.C. A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications // Computer methods in applied mechanics and engineering. - 2001. -Vol. 190. - P. 4379-4403.

30. Jannesari M., Mosaddegh P., Kadkhodaei M. et al. Numerical and clinical investigation on the material model of the cornea in Corvis tonometry tests: differentiation between hyperelasticity and viscoelasticity // Mechanics of Time-Dependent Materials. - 2019. - Vol. 23. - P. 373-384.

31. Kazaili A., Geraghty B., Akhtar R. Microscale assessment of corneal viscoelastic properties under physiological pressures // Journal of the Mechanical Behavior of Biomedical Materials. - 2019. - Vol. 98 (1). - P. 31-38.

32. Lago M.A., Ruperez M.J., Monserrat C. et al. Patient-specific simulation of the intrastromal ring segment implantation in corneas with keratoconus // Journal of the Mechanical Behavior of Biomedical Materials. - 2015. -Vol. 51. - P. 260-268.

33. Liu T., Shen M., Li H. et al. Changes and quantitative characterization of hyper-viscoelastic biomechanical properties for young corneal stroma after standard corneal cross-linking treatment with different ultraviolet-A energies // Acta Biomaterialia. - 2020. - Vol. 113. - P. 438-451.

34. Liu X., Wang L., Ji J. et al. A Mechanical Model of the Cornea Considering the Crimping Morphology of Collagen Fibrils // Investigative Ophthalmology & Visual Science. - 2014. - Vol. 55, No. 4. - P. 2739-2746.

35. Ljubimova D. Biomechanics of the Human Eye and Intraocular Pressure Measurements: Doctoral Thesis in Mechanics. - Stockholm: Royal Institute of Technology, 2009. - 200 p.

36. Lu J., Zhou X., Raghavan M.L. Inverse elastostatic stress analysis in pre-deformed biological structures: Demonstration using abdominal aortic aneurysms // Journal of Biomechanics. - 2007. - Vol. 40, No. 3. -P. 693-696.

37. Maklad A., Eliasy A., Chen K.-J. et al. Fluid-Structure Interaction Based Algorithms for IOP and Corneal Material Behavior // Frontiers in Bioengineering and Biotechnology. - 2020. - No. 8. - P. 970.

38. McQuaid R.M. Diffusion of oxygen and riboflavin during corneal cross-Linking (CXL): Doctoral Thesis. - Dublin: University College, 2017. - 119 p.

39. Meek K.M., Tuft S.J., Huang Y. et al. Changes in Collagen Orientation and Distribution in Keratoconus Corneas // Investigative Ophthalmology & Visual Science. - 2005. - Vol. 46, No. 6. - P. 1948-1956.

40. Montanino A., Angelillo M., Pandolfi A. A 3D fluid-solid interaction model of the air puff test in the human cornea // Journal of the Mechanical Behavior of Biomedical Materials. - 2019. - No. 94. - P. 22-31.

41. Mousavi S.J., Nassiri N., Masoumi N. et al. Finite Element Analysis of Blunt Foreign Body Impact on the Cornea After PRK and LASIK // Journal of Refractive Surgery. -2012. - Vol. 28, No. 1. - P. 59-64.

42. Nejad T.M., Foster C., Gongal D. Finite element modelling of cornea mechanics: a review // Arquivos Brasileiros de Oftalmologia. - 2014. - No. 77 (1). - P. 6065.

43. Nguyen B.A., Roberts C.J., Reilly M.A. Biomechanical Impact of the Sclera on Corneal Deformation Response to an Air-Puff: A Finite-Element Study // Frontiers in Bioengineering and Biotechnology. - 2019. - No. 6 - P. 210.

44. Pandolfi A., Holzapfel G.A. Three-Dimensional Modeling and Computational Analysis of the Human Cornea Considering Distributed Collagen Fibril Orientations // Journal of Biomechanical Engineering. - 2008. - Vol. 130. - P. 12.

45. Pandolfi A., Fotia G., Manganiello F. Finite element simulations of laser refractive corneal surgery // Engineering with Computers. - 2009. - No.25:15. - DOI: 10.1007/s00366-008-0102-5.

46. Pandolfi, A. Cornea modelling // Eye and Vision. - 2020. - No. 7 - P. 2.

47. Qiao X., Chen D., Huo H. et al. Full-field strain mapping for characterization of structure-related variation in corneal biomechanical properties using digital image correlation (DIC) technology // Medicine in Novel Technology and Devices. - 2021. - No. 11. - P. 100086.

48. Roy A.S., Dupps W.J. Effects of altered corneal stiffness on native and postoperative LASIK corneal biomechanical behavior: a whole-eye finite element analysis // Journal of Refractive Surgery. - 2009. - Vol. 25, No. 10. - P. 875887.

49. Roy A.S., Dupps W.J. Patient-specific modeling of corneal refractive surgery outcomes and inverse estimation of elastic property changes // Journal of Biomechanical Engineering. - 2011. - No. 113(1). - P. 10.

50. Roy A.S., Dupps W.J. Patient-specific computational modeling of keratoconus progression and differential responses to collagen cross-linking // Investigative Ophthalmology & Visual Science. - 2011. - Vol. 52, No. 12. - P. 9174-9187.

51. Roy A.S., Rocha K.M., Randleman J.B. et al. Inverse computational analysis of in vivo corneal elastic modulus change after collagen crosslinking for keratoconus // Experimental Eye Research. - 2013. - Vol. 113. - P. 92104.

52. Roy A.S., Dupps W.J., Roberts C.J. Comparison of biomechanical effects of small-incision lenticule extraction and laser in situ keratomileusis: Finite-element analysis // Journal of Cataract & Refractive Surgery. -2014. - Vol.40, No.6. - P. 971-980.

53. Roy A.S., Kurian M., Matalia H., Shetty R. Air-puff associated quantification of non-linear biomechanical properties of the human cornea in vivo // Journal of the Mechanical Behavior of Biomedical Materials. - 2015. -Vol. 48. - P. 173-182.

54. Simonini I., Angelillo M., Pandolfi A. Theoretical and numerical analysis of the corneal air puff test // Journal of the Mechanics and Physics of Solids. - 2016. - Vol. 93. -P.118-134.

55. Vellara H.R., Patel D.V. Biomechanical properties of the keratoconic cornea: a review // Clinical and Experimental Optometry. - 2015. - No. 100 - P. 103375.

56. Wang S., Chester S.A. Multi- physics modeling and finite element formulation of corneal UV cross- linking // Biomechanics and Modeling in Mechanobiology. - 2021. - Vol. 20. - P. 1561-1578.

57. Whitford C., Movchan N.V., Studer H., Elsheikh A. A viscoelastic anisotropic hyperelastic constitutive model of the human cornea // Journal of the Mechanical Behavior of Biomedical Materials. - 2018. - Vol. 17. - P. 19-29.

Финансирование. Исследование выполнено за счет гранта Российского научного фонда (проект № 22-21-20085). Конфликт интересов. Авторы заявляют об отсутствии конфликта интересов.

DEVELOPMENT OF A COMPLEX OF MATHEMATICAL MODELS OF THE CORNEA BIOMECHANICAL PARAMETRS WITH DIAGNOSED KERATO-CONUS BEFORE AND AFTER TREATMENT WITH CORNEAL COLLAGEN CROSSLINKING

E.G. Solodkova1, B.E. Malyugin2, I.N. Zakharov3, V.P. Bagmutov3, V.P. Fokin1, S.V. Balalin1, E.V. Lobanov1, V.H. Le3

1 Volgograd's National Medical Research Center «Eye Microsurgery» named after S.N. Fedorov, Volgograd Russia

2 National Medical Research Center «Eye Microsurgery» named after S.N. Fedorov, Moscow, Russia

3 Volgograd State Technical University, Volgograd, Russia

ABSTRACT

This study is aimed at creation of 3-D geometrical and mechanical finite-element models of the cornea and numerical analysis of its stress-strain state and mechanical behavior under conditions of intraocular pressure and external air jet pressure when diagnosing patients at different stages of keratoconus as well as after treatment with corneal collagen cross-linking.

There is created a geometrical model of the cornea in a form of spatial segment of thin-shell convex structure where shell thickness is variable and form of anterior and posterior surfaces is determined arbitrarily by extrapolation of experimental data of corneal topographic study of a specific patient with keratotopograph Pentacam AXL. Material, used for emulation of the cornea, is regarded as heterogeneous, isotropic and nonlinearly elastic. Its coefficients of rigidity are assigned basing on the known experimental data and are adjusted by correlation of parameters on corneal deformation under airflow impulse resulting from numerical modeling and their true estimation by noncontact tonometer Corvis ST. To describe airflow impact occurring in course of tonometry a distributed impulse strain, with parameters corresponding to the settings of airflow impulse of the applied equipment, is effectuated to a part of anterior surface of the emulated shell.

If ectatic process is present adjacent to the corneal region in question, local segments with diminishing coefficients of rigidity (from periphery to the center) are applied. Dimensions, form, place, quantity of such regions and character of rigidity decrease inside them are derived from solution of inverse problems aimed at minimization of difference between the results of emulation and experimental measurement of keratotopometric and deformational parameters in specific points of the cornea. To simulate corneal collagen crosslinking procedure, an additional circular zone with enhanced rigidity coefficients is provided. Modification of rigidity characteristics by depth and radius of this zone, is set in accordance with experimental data on photosensitive agent distribution in course of its diffusion within the cornea and flux density of UV emission.

©PNRPU

ARTICLE INFO

Received: 22 January 2022

Approved: 05 September 2022

Accepted for publication: 09 September 2022

Key words:

mathematical modeling, cornea, keratoconus, cross-linking, finite element method, Pen-tacam, Corvis

Аннотация научной статьи по медицинским технологиям, автор научной работы — Солодкова Е.Г., Малюгин Б.Э., Захаров И.Н., Багмутов В.П., Фокин В.П.

Похожие темы научных работ по медицинским технологиям , автор научной работы — Солодкова Е.Г., Малюгин Б.Э., Захаров И.Н., Багмутов В.П., Фокин В.П.

DEVELOPMENT OF A COMPLEX OF MATHEMATICAL MODELS OF THE CORNEA BIOMECHANICAL PARAMETRS WITH DIAGNOSED KERATOCONUS BEFORE AND AFTER TREATMENT WITH CORNEAL COLLAGEN CROSSLINKING