CC BY
141
Солодкова Е.Г., Малюгин Б.Э., Захаров И.Н., Лэ В.Х., Фокин В.П., Балалин С.В., Лобанов Е.В. Идентификация параметров модели роговицы с кератоконусом при численно-экспериментальном исследовании ее геометрии и механического поведения. Российский журнал биомеханики, 2023, Т. 27, № 3, С. 67-80. БО!: 10.15593/Я7ЬВютесЬ/2023.3.05
V. РОССИИСКИИ ЖУРНАЛ БИОМЕХАНИКИ
пермскии №3 2023
П 0/1И тех RUSSIAN JOURNAL OF BIOMECHANICS
https ://ered.pstu. ru/index.php/rjb
Научная статья
DOI: 10.15593/RZhBiomech/2023.3.05 УДК: 531/534: [57+61]
ИДЕНТИФИКАЦИЯ ПАРАМЕТРОВ МОДЕЛИ РОГОВИЦЫ С КЕРАТОКОНУСОМ ПРИ ЧИСЛЕННО-ЭКСПЕРИМЕНТАЛЬНОМ ИССЛЕДОВАНИИ ЕЕ ГЕОМЕТРИИ И МЕХАНИЧЕСКОГО ПОВЕДЕНИЯ
Е.Г. Солодкова1, Б.Э. Малюгин2, И.Н. Захаров3, В.Х. Лэ3, В.П. Фокин1, С.В. Балалин1, Е.В. Лобанов1
1 Волгоградский филиал Национального медицинского исследовательского центра, межотраслевой научнотехнический комплекс «Микрохирургия глаза» им. акад. С.Н. Федорова Минздрава России, Россия
2 Национальный медицинский исследовательский центр, межотраслевой научно-технический комплекс «Микрохирургия глаза» им. акад. С.Н. Федорова Министерства здравоохранения Российской Федерации, Москва, Россия
3 Волгоградский государственный технический университет, Волгоград, Россия
О СТАТЬЕ
АННОТАЦИЯ
Получена: 30 марта 2023 Одобрена: 11 сентября 2023 Принята к публикации: 12 сентября 2023
Ключевые слова:
математическое моделирование, роговица, кератоконус, кросслинкинг, метод конечных элементов, Pentacam,
Целью работы является построение расчетной методики идентификации параметров модели механического поведения роговицы с использованием трех клинических биомаркеров в качестве входных данных: амплитуда деформаций роговицы в реальном времени, измеренная в ходе бесконтактной пневмотонометрии, внутриглазное давление и геометрические особенности роговицы пациента, полученные при его томографическом обследовании.
Разработанный инструмент прогнозирования свойств роговицы в области кератоконуса и за ее пределами опирается на набор данных, сформированный по результатам конечно-элементного моделирования бесконтактного тонометрического теста, и процедуру многопараметрической оптимизации (минимизации) отклонения расчетных результатов от измеренных данных о амплитуде деформации роговицы, полученных на пневмотонометре Сомб ЭТ. Моделирование основано на результатах обследования 78 пациентов реальной клинической базы данных с учетом специфических для каждого пациента геометрии роговицы, размеров и положения зоны кератоконуса, внутриглазного давления и деформационных параметров материала при пневмотесте.
©ПНИПУ
© Солодкова Елена Геннадиевна- к.м.н., заместитель директора по научной работе, e-mail: solo23el@mail.ru ¡D: 0000-0002-7786-5665
© Малюгин Борис Эдуардович - д.м.н., заместитель генерального директора по научной работе, e-mail: boris.maluain@gmail.com 0000-0001-5666-3493
© Захаров Игорь Николаевич - д.т.н., заведующий кафедрой, e-mail: zaxap@mail.ru ¡D 0000-0001-7177-7245 © Лэ Ван Хоанг - аспирант кафедры: hoanale.vol@gmail.com 0000-0002-1536-3061 © Фокин Виктор Петрович - д.м.н., директор филиала, e-mail: fokin@isee.ru 0000-0002-2513-9709 © Балалин Сергей Викторович - д.м.н., заведующий научным отделом, e-mail: s.v.balalin@gmail.com ¡D 0000-0002-5250-3692
© Лобанов Евгений Валерьевич - инженер, e-mail: lobanoff95@vandex.ru 0000-0001-9112-3230
Эта статья доступна в соответствии с условиями лицензии 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)
Введение
Знание механических свойств глазных тканей, в частности роговицы, необходимо для корректной оценки параметров, характеризующих физическое состояние глаза (включая определение внутриглазного давления), для уточнения диагностических методов, а также при планировании операций [2]. Между тем роговица обладает сильно индивидуальными, весьма сложными и при этом недостаточно изученными механическими свойствами [2].
На разных уровнях исследования роговица может рассматриваться как гиперупругий, существенно нелинейный материал, наделяемый характеристиками неоднородности, анизотропии по различным направлениям, вязкоупругости и т.д. [8, 18, 21, 25 35, 37, 38, 46,]. Кроме того, роговица неоднородна геометрически - ее толщина и кривизна существенно меняются от точки к точке и от апекса к периферии [2]. Все это делает задачу определения механических свойств роговицы крайне нетривиальной даже в случае нормального ее состояния, наличие же отклонений и паталогических изменений еще более усложняет картину.
Важно знать и понимать ограничения и преимущества различных методик, используемых для получения экспериментальных и теоретических (расчетных) данных, присущие им источники возможных ошибок для определения степени надежности полученных результатов. Особое значение результаты натурных (клинических) экспериментов приобретают при построении математических моделей роговицы (и глаза в целом), которые могут быть чрезвычайно чувствительны к исходным параметрам и предположениям [26].
В работе [26] приводится обзор исследований по механике роговицы, которые делятся на две группы: исследования in vivo и исследования in vitro на различных образцах (как интактных, так и разрезанных на полоски, диски и т.д.).
Показано, что многие методы оценки биомеханических параметров роговицы in vivo (например, голографическая интерферометрия [19], контактная тонометрия при динамическом индентировании [12], ультразвуковая микроскопия [14]) являются технически очень сложными, часто - трудно переносимы пациентами (иногда - на пределе безопасного применения), поэтому большинство из них не реализованы в клинике. В настоящее время наибольшее распространение получили
диагностические методики на основе бесконтактного пневмотонометрического воздействия на роговицу (ORA, Corvis и др.), в том числе, с визуализацией и
измерением ее деформаций при помощи Шаймпфлюг-метода [26, 35, 46].
Исследования in vitro в основном базируются на испытаниях полосок роговицы на растяжение [13, 24, 30, 43, 47, 48, 50], либо испытаниях роговичного диска [10, 49] или кадаверного глаза в целом [20] под действием нагнетаемого давления.
При этом в подобных исследованиях может наблюдаться достаточно широкий разброс между свойствами материалов (отличия достигают нескольких порядков). Это может быть связано с различиями в методиках подготовки и испытания образцов; различиями в уровнях напряжений или деформаций, при которых определяются свойства, что особенно актуально учитывая нелинейный характер деформаций; различиями в скорости деформации при испытании [26]. Кроме того, следует учитывать существенные индивидуальные различия в свойствах ткани для разных пациентов, глаз (тем более при использовании донорского материала) и т.д.
Отдельную группу составляют исследования свойств роговицы in silico, то есть на основе компьютерного моделирования ее структуры и механического поведения. Основу таких работ составляют решения обратных задач по определению параметров материала роговицы, обеспечивающих согласование экспериментальных данных и расчетных результатов о геометрии и деформациях роговицы, полученных с помощью, пациент-ориентированных моделей на основе метода конечных элементов (МКЭ). Количество подобных работ сравнительно невелико, при этом расчетное определение числовых значений параметров, характеризующих свойства роговицы, можно найти лишь в некоторых из них [3-5, 7, 11, 16, 17, 22, 23, 32, 36, 40, 41].
В работах [5, 16, 17, 22, 23, 36, 40,] для идентификации параметров роговицы, представляющих собой коэффициенты в принятых законах деформирования материала, используются
экспериментальные результаты о его деформациях в ходе бесконтактного пневмотонометрического теста (при помощи Corvis ST). В частности, в [5, 16, 23, 40] роговица рассматривается как гиперупругий материал: в форме модели Муни-Ривлина с одним [16] или двумя [5, 40] параметрами или в форме модели Огдена с 4 параметрами [23]. При этом в [16, 40] рассматривается анизотропия свойств роговицы за счет дополнительных параметров, учитывающих влияние плотности распределения коллагеновых волокон. Указанные параметры устанавливаются на базе оптимизационных процедур (например, с использованием алгоритма Левенберга-Марквардта [40]), минимизирующих разницу между расчетными и экспериментальными данными о амплитуде деформаций роговицы в ходе
действия воздушной струи пневмотонометра. Подобным же образом устанавливаются характеристики вязкоупругого материала в соответствующих моделях роговицы - например, для модели Прони [22] или Бюргерса [17]. Модель линейно-вязкоупругого материала используются в работе [36] для идентификации параметров здоровой роговицы и роговицы с кератоконусом. При этом для больной роговицы устанавливаются пониженные значения свойств, равномерно распределенных по ее объему.
Отметим, что в большинстве из приведенных работ рассматриваются упрощенные одномерные [17], двумерные [5, 22, 36], либо осесимметричные (сфера [23], эллипсоид [16]) геометрические модели роговицы. Полноценная трехмерная геометрия передней и задней поверхностей роговицы, построенная по данным Шаймпфлюг топографии (Pentacam) используется нечасто [40].
Другая группа расчетных методов оценки механических свойств роговицы базируется на математическом моделировании экспериментов in vitro по деформированию роговицы или глаза (на специальных установках) под действием нагнетаемого внутриглазного давления [4, 7], при индентировании [32] или другими способами. В таких случаях часто используются достаточно сложные физические модели анизотропного (по характерным направлениям распределения фибрилл) гиперупругого материала роговицы (модели Муни-Ривлина [7], Гассера-Холцапфеля-Огдена [4], Нео-Гука [32]) и трехмерные геометрические модели ее формы.
В работе Roy [41] параметры гиперупругого материала (для модели Муни-Ривлина) роговицы с кератоконусом, а также степень их повышения после операции кросслинкинга устанавливаются из итерационной процедуры согласования расчетной кривизны передней поверхности с экспериментальными данными, полученными при помощи кератотопографа Pentacam AXL. При этом в области кератоконуса (заданной при моделировании МКЭ в виде круговой области радиусом 3,5 мм с центром в точке максимальной кривизны передней поверхности) с использованием опытных данных работ [3, 30] устанавливается падение упругих модулей на 50 % (для образцов роговиц после пересадки). По результатам оптимизации расчетное приращение модулей упругости для рассмотренных моделей роговиц после кросслинкинга составило от 75 до 150 % [41].
Как видим, проблема определения механических характеристик роговицы в приложениях к диагностике и лечению ее заболеваний (в частности, кератоконуса) достаточно остро стоит перед современной офтальмологией. При этом требуют развития как тонкие экспериментальные и клинические методы
локального изучения свойств и структуры роговицы, так и расчетные методики восстановления и идентификации параметров ткани на основе решения обратных задач при компьютерном моделировании процессов деформирования роговицы при различных воздействиях.
Постановка задачи. Основные соотношения и граничные условия
Этапы построения и базовые соотношения пациент-ориентированной трехмерной конечно-элементной модели биомеханических свойств и поведения роговицы подробно рассмотрены авторами в предыдущей статье [1]. Поэтому в этом разделе ограничимся кратким изложением наиболее характерных особенностей данной модели.
Первый уровень индивидуализации модели заключается в построении геометрического 3^-прототипа роговицы (с использованием компьютерной системы конечно-элементного моделирования (СКМ) на основе результатов топографических, томографических исследований конкретного пациента с помощью кератотопографа Pentacam AXL. Производится линейная интерполяция экспериментально найденного массива точек (700010000 точек) ограничивающими поверхностями (передней и задней), «натянутыми» на указанный массив высот (координат), как на каркас, с последующей конечно-элементной дискретизацией полученной расчетной области. Исходными данными к задаче являются: максимальные диаметры роговицы в назально-темпоральной dт™ и нижне-верхней dт™
плоскостях, толщина роговицы в апексе tapex, глубина передней камеры hc. Другие геометрические и топографические параметры роговицы (определяемые экспериментальными картами толщины,
тангенциальной и сагитальной кривизны) используются для верификации расчетной геометрии построенной 3D модели.
Материал роговицы считается гиперупругим, его закон деформирования задается соотношением функции плотности энергии упругой деформации Ws и второго тензора напряжений Пиолы-Кирхгофа (здесь и далее используются обозначения, принятые в СКМ Сот^'о1 Multiphysics) [38, 39]:
5 = S„, +-
ж
(1)
где 5 - тензор внешних (дополнительных) напряжений; Е - тензор деформации Грина-Лагранжа. В общем случае функция Ws может учитывать
соотношения, соответствующие различным механизмам деформирования материала, включая, например [15, 34], составляющие для изотропной части удельной
Рис. 1. Общий вид геометрической модели роговицы, схема нагружения и закрепления граничных поверхностей: а - поверхность примыкания Гх; б - задняя поверхность Г6; в - передняя поверхность Г^
энергии упругой деформации при постоянном объеме Ж0 , для энергии объемной упругой деформации , для энергии деформации коллагеновых волокон Ж^ (при описании анизотропных свойств роговицы), для «диссипативной» составляющей Ж^ (при описании вязкостных свойств роговицы):
Ж = Ж^ + ЖуЫ + ЖГЛ + Жь, (2)
Общая постановка рассматриваемых задач будет задаваться системой уравнений равновесия и уравнений совместности деформаций в виде, приведенном в [1], физическими соотношениями (1), (2), а также граничными условиями, заданными на соответствующих поверхностях Г роговицы (рис. 1):
а) отсутствие перемещений на поверхностях Г , примыкающей к склере
и = 0, ГеГх ; (3)
б) действие на задней поверхности Гь роговицы
поверхностно распределенной нагрузки
интенсивностью д4 (вдоль составляющих вектора нормали п)
Б • п = дь, ГеГ4; (4)
в) действие на передней поверхности Г^ роговицы
поверхностно распределенной нагрузки
интенсивностью ^
S ■ n = I
■■Ч/, ГеГ /• (5)
Вид функций для ч и в условиях (4), (5)
устанавливается в зависимости от решаемой задачи. Например, ч может соответствовать внутриглазному давлению (ВГД) пациента при решении прямой задачи или некоторой компенсирующей нагрузке при решении обратной задачи, тогда как ^ используется в задачах о
действии на переднюю поверхность роговицы воздушного импульса при бесконтактной тонометрии.
В частности, согласно [31], для определения исходной, недеформированной конфигурации роговицы
к ее задней поверхности прикладывается фиктивное давление pim,, обратное по знаку ВГД (p ) и равное
Р,™ = —^Рир . Здесь X - коэффициент, учитывающий
изменение площади задней поверхности роговицы в деформированном и недеформированном состояниях (принимали X = 0,96 -1). Найденная из численного расчета геометрия роговицы экспортируется в качестве исходной в новую модель, при этом поля напряжений и деформаций обнуляются.
На полученных границах задаются краевые условия, соответствующие прямой задаче определения текущего НДС материала - задняя поверхность роговицы нагружается ВГД pfo . Поверхностные силы
q в условии (4) в таком случае будут задаваться в
виде:
Чь = -Piop ■ n ГеГь . (6)
Полученная из решения прямой задачи расчетная геометрия роговицы сопоставляется с
топографическими и томографическими картами пациента (карты высот, толщин, кривизны передней и задней поверхностей). При этом используются различные приборы (например, кератотопографы Pentacam на (рис. 2, а, б) или Sirius на (рис. 2, в, г) и методики обработки первичных данных (с использованием или без использования аппроксимации полиномами Цернике (рис. 2, е-з).
Для моделирования действия струи воздуха при бесконтактной пневмотонометрии (Corvis ST) на передней поверхности роговицы, согласно условию (5), задается нагрузка q , соответствующая давлению
воздушного импульса Р
qf = - Pai
ГеГ
f ■
(7)
Распределение давления p^ (x, y, z, t) по передней поверхности роговицы представляется в виде сочетания функции координат f (y, z) и времени f (t):
Parr ^ y, z, t) = Pmx • f (У, z)-f (t), x eTf . (8)
д е ж з
Рис. 2. Экспериментальные (а—г) и расчетные (д-з) карты распределения толщины (а, д), тангенциальной (б, в, е, ж) и сагиттальной (г, з) кривизны передней поверхности роговицы с кератоконусом: а, б - данные Pentacam; в, г - данные Sirius; д, ж - без использования аппроксимации; е, з - после аппроксимации Цернике
Здесь р_ - максимальное давление в центре воздушного импульса.
Величина ртах, вид и временной профиль £(£)
давления струи задается путем интерполяции экспериментальных данных, фиксируемых в ходе тонометрического теста на базе Corvis БТ для 140 точек исследуемого интервала времени.
Для описания функции пространственного распределения давления £ (у, г) используется
аппроксимация расчетных результатов работ [5, 16, 17, 22, 28, 40] в виде [1]
£ (у,г) = ехр[бр (у2 + г2)/К,/2}2 ], (9)
где йа1г - диаметр сопла прибора (= 2,4 мм [42]); 0 - коэффициент формы импульса по результатам аппроксимации (в первом приближении 0 = -0,75).
Модель материала роговицы
Выбор модели гиперупругости материала роговицы
Механические характеристики роговицы и их изменения на различных стадиях заболевания остаются неизвестными для данного пациента и практически не поддаются надежному экспериментальному
определению (на современном уровне развития технических средств диагностики). Использование усредненных характеристик, назначаемых, например,
по литературным данным, оказывается малоэффективным с точки зрения пациент-ориентированного подхода, а в ряде случаев -невозможным, в виду существенных индивидуальных различий свойств роговицы и течения заболеваний у различных пациентов. Возникает необходимость решения обратных задач отыскания неизвестных параметров модели материала роговицы конкретного пациента.
Закон деформирования материала задается выражением (2). Проведенные тестовые расчеты показали, что учет анизотропии свойств роговицы ^ь
и ее вязкостных характеристик Wv¿s приводит к некоторому уточнению получаемых расчетных результатов (в пределах 5-10 %) относительно экспериментальных данных, но сопровождается многократным увеличением времени счета, критичным для разрабатываемых в статье процедур. Поэтому ограничимся первыми двумя слагаемыми (, ) в выражении (2), вид которых устанавливается следующим образом.
Составляющая
соответствует энергии
объемной упругой деформации [33, 34]
1 9
Wvol =- к(7-1) , (10)
где J = det (- определитель градиента деформации
¥, характеризующий изменение объема, вызванное деформацией; К - модуль объемной упругости (в да-
О 5 10 15 20 25 г, мс Рис. 3. Сопоставление расчетных и экспериментальных данных по амплитуде деформаций роговицы при действии воздушного импульса: а - для некоторых моделей гиперупругости (1 - Муни-Ривлина с 2 параметрами; 2 - Муни-Ривлина с 5 параметрами; 3 - Огдена; 4 - Йо); б - для различных алгоритмов оптимизации параметров материала (I -ММА; II- 1РОЯТ; III - БОБУдА; IV- СОБУ1А; V- ЫеМег-Меа^; в-е - картины деформаций в назально-темпоральном сечении роговицы
нной работе принимается постоянным - по данным [34] К = 5,5 МПа).
Слагаемое W¡so задается в зависимости от принятой модели гиперупругого деформирования материала роговицы, например, Муни-Ривлина, Огдена, Сторакерса, Йо и т.д. [4, 5, 7, 16, 32, 40, 41].
Для выбора конкретной модели проводится серия сравнительных вычислительных экспериментов о деформациях роговицы при пневмотонометрическом тесте. Задаются распределения нагрузок: на задней поверхности - внутриглазное давление, согласно (4), (6); на передней - переменное давление струи, согласно (5), (7) с учетом (8), (9). Находятся закономерности деформации роговицы во время действия воздушного импульса, в том числе, строятся кривые изменения амплитуды в апексе (рис. 3). Из сопоставления полученных расчетных результатов и клинических данных пациентов (более 30) установлено, что наименьшие отклонения (для рассмотренных моделей гиперупругости - Муни-Ривлина с 2 и 5 параметрами, Огдена, Сторакерса, Йо) соответствуют модели Йо (рис. 3, г).
В рамках модели Йо функция плотности энергии деформации IV. в соотношениях (2) может быть
записана в виде [6, 39]
W = с,
= с (Л - 3)+с (I - 3) + Сз (I - 3)
(11)
где с, с2, с - эмпирические параметры модели; I - первый инвариант (след) тензора С
I, = ^ (С), С = J-2/3 С , (12)
С - правый тензор деформаций Коши-Грина.
Модель материала в зоне кератоконуса
Для описания патомеханики кератоконуса в объеме роговицы вводится локальная область Qfc пониженной жесткости с центром в заданной точке K0 (рис. 4). Считается, что точке Ko соответствует вершина кератоконуса и максимальное снижение абсолютных величин коэффициентов «жесткости» с, с2, с3 в модели гиперупругого материала (11). По мере удаления от центра Ko кератоконуса к периферии роговицы задается плавное повышение параметров с,
с2, с , например, при помощи функции yfc :
Vfc = Vrnx • exp[e^x2/al +Qkry2kjalr +Qk2z2k/al ] (13) где x, У, z* - локальная система координат, связанная с центром зоны кератоконуса - точкой K0; Ушах - максимальное относительное снижение
жесткости в центре зоны кератоконуса; e , e , e -параметры, устанавливающие градиент изменения
а б в
Рис. 4. Зона пониженной жесткости в области кератоконуса в объеме геометрической модели: а - передняя поверхность; б - задняя поверхность; в - в сечениии по толщине свойств по рассматриваемой области вдоль
соответствующей оси; ак
ак - параметры,
устанавливающие размер зоны кератоконуса вдоль каждой из осей.
Таким образом, в случае роговицы с кератоконусом выражение (11) для упругого потенциала с учетом функции (13) можно записать так [38]:
Жко = (1 - Чс) с (Л - 3) + с (I, - 3)2 + С3 (I, - 3)3
. (14)
Приняты следующие значения коэффициентов в (13): е^ = е^ = е^ = -0,693; ^ = а^ = аК = Як, что
соответствует пятидесятипроцентному снижению значений функции на одинаковом расстоянии Кк от точки К0.
Идентификация параметров модели
Для определения неизвестных биомеханических свойств материала решается серия обратных задач по отысканию параметров модели в случаях здоровой роговицы (11) и роговицы с кератоконусом (14).
Используется процедура многопараметрической оптимизации (встроенный модуль СКМ Сош^'о1 Multiphysics), состоящая в отыскании минимума целевой функции Ъ (с, с2, с3, ^, К ), в неявном виде
содержащей параметры материала роговицы в интактных областях (с, с2, с3) и в зоне кератоконуса (V™, К). Данная функция формируется в виде суммы квадратов отклонений расчетных величин амплитуды деформаций ВЛс от измеренных на пневмотонометре Corvis БТ действительных деформаций БЛР роговицы пациента, полученных для каждой из п=140 контрольных точек временного интервала действия воздушного импульса:
Ъ( с!
^ с2 , ^ V,
ттах^ Ч
п 2
К) = Ё{°ЛС -Щ") ^тт. (15)
Выбор конкретного алгоритма минимизации целевой функции (например, Монте Карло, Левенберга-Марквардта, Нелдера-Мида, координатного спуска, внутренней точки и т.д.; (см. рис.3, б) производится на основе серии вычислительных экспериментов по наименьшему отклонению расчетных и опытных данных, а также времени счета. По указанным критериям в качестве основного для здоровых роговиц был выбран алгоритм Левенберга-Марквардта (ЬеуепЬе^-МагдиагЛ), для роговиц с кератоконусом -метод Нелдера-Мида (ИеМег-Меа^ [6].
С учетом вышеизложенного процедура идентификации параметров механической модели здоровой роговицы (11) выглядит следующим образом.
1. На основе данных кератотопографии переднего отрезка глаза пациента (при помощи РеШасаш ЛХЬ) строится индивидуализированная геометрическая 3Б-модель роговицы, задаются ее геометрические параметры.
2. В первом приближении механические характеристики материала назначаются по результатам серии предварительных расчетов, как средние для здоровых роговиц (в нашем случае: с = 0,225 МПа; с2 = -28,3 МПа; с = 2460 МПа).
3. Решается обратная задача, по результатам которой устанавливается форма и размеры роговицы в недеформированном состоянии.
4. Для найденной геометрии недеформированного состояния роговица нагружается ВГД в соответствии с (6) - устанавливается естественная конфигурация роговицы под действием ВГД и соответствующее ей напряженно-деформированное состояние (НДС).
5. При необходимости расчет повторяется с п.3 до сходимости с заданной погрешностью отклонения апекса от действительного положения.
6. После того, как естественная конфигурация и начальное НДС роговицы установлены с необходимой точностью, ставится задача (1)-(5) о деформациях роговицы под действием давления воздушного импульса при тонометрии с учетом (8), (9). Данная задача решается, как обратная, то есть производится поиск параметров с , с , с , соответствующих минимуму целевой функции (15), полученной для данного пациента при обследовании на пневмотонометре Corvis БТ. Для минимизации целевой функции Ъ при вариации трех параметров с, с, с на
данном этапе используется алгоритм Левенберга-Марквардта. При этом для здоровой роговицы параметры и К считаются равными нулю и из
анализа исключаются, то есть Ъ(с, с2, с3).
7. Для найденных коэффициентов с, с2, с3 форма роговицы может быть дополнительно уточнена по
а
г=1
а б в
Рис. 5. Расчетные картины передней поверхности здоровой роговицы под действием ВГД: а- перемещения; б - интенсивность напряжений; в - интенсивность деформаций
Рис. 6. Схема определения размеров зоны кератоконуса и соотношения отрезков назально-темпоральной оси (NT) роговицы, приходящихся на области пониженной (KA, AC) и нормальной (NK, CT) жесткости
процедуре 3-5 с повторным решением задачи оптимизации до достижения заданной точности.
На рис. 5 показаны некоторые результаты расчетов - картины распределения перемещений, эквивалентных напряжений по Мизесу и интенсивности деформаций по задней поверхности здоровой роговицы под действием ВГД.
Для роговицы с кератоконусом процедура идентификации свойств материала усложняется за счет увеличения количества свободных параметров целевой функции (15) с трех до пяти Z(q,с2,с3,vmax,Rk). В
таком случае решение задачи ищется в два этапа при следующем порядке операций.
1. На первом этапе описанная выше процедура используется для поиска «усредненных» коэффициентов упругости с*, с*2, с*ъ эктатической
роговицы с условно равномерным изменением жесткости по ее объему.
2. На втором этапе вводится локальная область кератоконуса эффективным радиусом R с
пониженными кратно коэффициентами с, с2, с3 . На удалении от этой зоны роговица считается неизменной.
Граница зоны кератоконуса ( R ) предварительно
устанавливается при помощи полученной с использованием Pentacam пахиметрической карты роговицы, как внешний контур желтого цвета по интуитивной шкале Белина (рис. 6).
Устанавливается соотношение между отрезками, приходящимися на области пониженной (KA, AC) и нормальной (NK, CT) жесткости (рис. 6), в сечении NT роговицы назально-темпоральной плоскостью.
где Рис. 7. Сопоставление экспериментальных результатов с расчетными на роговице с кератоконусом под действием ВГД: а, г - карта толщин; б, д - передная кривизна, в, е - задняя элевация
В зависимости от величины параметров с*, с\, с3 и принятой в «нулевом» приближении величины Ушах =^тах, вычисляются ориентировочные значения коэффициентов упругости с\, с\, с30 роговицы за пределами зоны кератоконуса:
" NK CT' +(i-vL )• " KA AC'
_ "NX _|__ AT _ __|__ _ NA AT _
Начальное значение коэффициента
максимального снижения параметров упругости в зоне кератоконуса задается в зависимости от стадии заболевания следующим образом: I стадия -у0 = 0,25; II - у0 = 0,5; III - у0 = 0,75.
Tmax ' > Tmax ' > т max '
Начальный эффективный радиус R° зоны кератоконуса принимается, как наибольший из отрезков KA или AC: R = max (KA, AC) (рис. 6)
3. Используя найденные значения q0, c\, c30, у0ах и R° в качестве начальных, решается обратная задача (1)-(9), (15) многопараметрической оптимизации (алгоритм Нелдера-Мида) для роговицы с кератоконусом. В качестве центра зоны с пониженными свойствами в разных сериях расчетов выбирались точки: минимальной пахиметрии, максимальной задней
элевации, максимальной кривизны передней поверхности. Поиск пяти неизвестных свободных параметров с, с2, с3, у^, производится путем
минимизации целевой функции ,с2,с3,у^,Як) из
сопоставления расчетной амплитуды деформаций с данным согу1&' БТ для конкретного пациента (рис. 3).
Расчетные результаты и их анализ
На рис. 7, а-в показаны экспериментальные карты толщины, тангенциальной кривизны передней поверхности и задней элевации роговицы, полученные в ходе дооперационной диагностики пациента при помощи анализатора Рв^асат. Соответствующие расчетные карты толщины, кривизны и картины распределения интенсивности деформаций под действием ВГД приводятся на (рис. 7, г-е).
Сопоставление экспериментальных и расчетных картин показывает, что локальной области роговицы с наименьшей толщиной (желтая область по шкале Белина; (рис. 6; рис. 7, а, г) соответствуют наибольшие значения деформаций (50-100 % от максимума; на (рис. 7, е) - темнокрасная область), тогда как остальная
д е ж з
Рис. 8. Карты кривизны задней поверхности а-в и интенсивности деформаций г-е роговицы под действием ВГД для различных положений центра зоны снижения упругих характеристик: в точке минимальной пахиметрии ртт а, г, максимальной задней элевации етах б, д и максимальной кривизны задней поверхности Ктах в, е
часть роговицы мало деформирована (0-50 % диапазона значений).
Точка максимальной интенсивности деформаций достаточно случайным образом может локализовываться вблизи любой из характерных точек роговицы - минимальной пахиметрии (pmin), наибольшей задней элевации (emax), максимальной кривизны (Kmax) передней поверхности. По-видимому, ни одна из этих точек не может однозначно ассоциироваться с центром области пониженных свойств роговицы с кератоконусом. Формирование этой зоны определяется как уменьшением толщины роговицы (геометрическое разупрочнение), так и снижением ее упругих и прочностных свойств (физическое разупрочнение), при этом вершина конуса будет располагаться на промежутке между этими двумя точками. В этой связи существующие гипотезы о том или ином расположении центра наиболее ослабленной области - в зоне минимальной пахиметрии [27, 44], максимальной задней элевации [9, 42], максимальной кривизны передней [29, 45] или задней [27] поверхностей требуют дальнейшего исследования.
На (рис. 8) приведены экспериментальные и расчетные карты тангенциальной кривизны (а-г), пахиметрии (д) и интенсивности деформаций (е-з) для трех возможных случаев локализации зоны пониженной жесткости - в точке минимальной пахиметрии, максимальной задней элевации и максимальной кривизны передней поверхности. Показано (рис. 8), что изменение положения области пониженных свойств принципиально не изменяет ни основные расчетные топографические параметры (взаимное отклонение
величин - не более 5 %), ни распределение параметров НДС материала (отклонение до 15 %).
В результате вычислительных экспериментов для роговиц с кератоконусом разных стадий, а также здоровых роговиц (общее количество пациентов всех групп - 78 чел.) были установлены значения упругих коэффициентов с , с , с для моделей (11), (14) гиперупругого поведения материала вне зоны кератоконуса (рис. 9, а). Показаны средние значения каждого из коэффициентов (толстая линия), границы первого и третьего квартиля (прямоугольная область) и среднеквадратическое отклонение (тонкая линия с засечкой).
На (рис. 9, б) приведены расчетные данные о характеристиках зоны пониженной жесткости в зависимости от стадии заболевания - максимальное снижение коэффициентов упругости в центре зоны Утах, ее эффективный радиус Rk, а также отношение
деформаций в роговице MMSI
MMSI =Smax/ S,
(17)
Как видим, за пределами зоны кератоконуса средние расчетные значения параметров упругости с , с , с для различных стадий заболевания мало отличаются от характеристик здоровой роговицы, укладываясь в диапазон стандартного отклонения
(рис. 9). Для роговиц с кератоконусом наблюдается разброс значений коэффициентов модели, что связано с отсутствием четких границ области пониженной жесткости и нестабильной динамикой деформаций экта-
максимальной е v и минимальной s • интенсивности
max ^Qin
с2х(-107) а х 109 V™. Л,
а б
Рис. 9. Зависимости: а - упругих коэффициентов С1, с2, с3 для модели гиперупругого поведения материала роговицы; б- характеристик зоны пониженной жесткости - коэффициента снижения жесткости утах, эффективного радиуса Як, отношения максимальной и минимальной интенсивности деформаций ММБ1
в роговице от стадии заболевания
тической роговицы в ходе пневмотеста.
С развитием заболевания возрастает коэффициент ^nax падения свойств в зоне кератоконуса (I стадия -
V m
¥ m
% ; II стадия - ¥ CL ~ 11 ; Ш стадия -
! 21 % ), так же как и соотношение максимальной и минимальной интенсивности деформаций MMSI (с 2,2 до 4,05). При этом радиус Rk данной зоны уменьшается - процессы максимальной деградации свойств и роста деформаций локализуются вокруг вершины конуса. Повышенный разброс данных может объясняться как индивидуальными особенностями протекания заболевания (разбросом свойств) у различных пациентов, так и более высокой погрешностью экспериментальных данных при обследовании больных на поздних стадиях кератоконуса.
Заключение
По результатам проделанной работы можно сделать следующие заключения:
1. Обзор публикаций показал, что проблема идентификации свойств роговицы в различных состояниях и с разной степенью паталогических изменений является весьма актуальной для современной офтальмологии. Приходится
констатировать отсутствие однозначных клинических и экспериментальных данных о величине и распределении по объему упругих и прочностных свойств роговицы в нормальном состоянии и, тем более, при наличии, например, эктатического процесса. Возможности современного испытательного оборудования остаются ограниченными, как в плане построения пространственных карт распределения свойств, так и повторяемости результатов на широкой базе пациентов, разных приборах и методиках. В этой
связи большое значение приобретают методы расчетного определения характеристик роговицы на базе пациент-ориентированных математических моделей и вычислительных экспериментов, воспроизводящих распространенные натурные методики испытаний - контактной и бесконтактной тонометрии, растяжения образцов, действия давления и т.д. Перспективными представляются модели развития эктатических процессов, позволяющие описать динамику изменения структуры и свойств ткани по мере развития заболевания. 2. Разработанная на базе компьютерной системы конечно-элементного моделирования Сот'о1 МиШрИус трехмерная модель роговицы включает несколько уровней персонализации, позволяющих с высокой степенью достоверности получать расчетные результаты, соответствующие клиническим данным реальных пациентов, и, как следствие, рассматривать ее в качестве пациент-ориентированного вычислительного инструмента исследования специфики заболеваний и методик их лечения.
Указанные уровни персонализации включают прямое использование приборных клинических данных пациента, импортируемых в модель из рабочих файлов кератотопографа Рв^асат ЛХЬ и бесконтактного пневмотонометра Согу!' БТ. Так, геометрическая 3^-модель роговицы строится путем интерполяции координат точек передней и задней поверхности, полученных из Рв^асат ЛХЬ, а параметры ее нагружения - внутриглазное давление и давление в центре воздушной струи в каждый момент времени в ходе пневмотеста задаются по данным СОГУ!' БТ.
3. Ряд параметров моделей роговицы остаются неизвестными (например, геометрия ненагруженного состояния, параметры механического поведения материала и их распределение по объему роговицы в
норме и паталогии, распределение давления воздушной струи по поверхности роговицы в различные моменты времени и т.д.), поскольку их прямое экспериментальное определение либо невозможно, либо крайне затруднительно.
Для отыскания таких параметров формулируется и решается ряд обратных задач, в которых в качестве исходных данных используется экспериментальная информация о текущем состоянии роговицы и ее поведении под нагрузками, а искомые неизвестные параметры модели устанавливаются по результатам серии вычислительных (виртуальных) экспериментов на модели и их соответствия опытным данным.
В частности, исходная геометрия роговицы восстанавливается из решения задачи о действии на заднюю поверхность роговицы давления, обратного внутриглазному. При этом величина поправочного коэффициента, учитывающего изменение геометрии нагружаемой поверхности, уточняется в итерационном процессе.
Идентификация параметров модели гиперупругого поведения материала роговицы связана с решением достаточно сложной задачи, так как требует отыскания нескольких параметров (в нашем случае, трех - для здоровой роговицы и пяти - для роговицы с кератоконусом), которые не поддаются прямому экспериментальному определению. В этом случае используется процедура многопараметрической
Список литературы
1. Солодкова Е.Г., Малюгин Б.Э., Захаров И.Н., Багмутов В.П., Фокин В.П., Балалин С.В., Лобанов Е.В., Лэ В.Х. Разработка комплекса математических моделей биомеханических параметров роговицы с диагностированным кератоконусом до и после лечения кросслинкингом роговичного коллагена // Российский журнал биомеханики. - 2022. - № 3. - С. 10-28.
2. Штейн А.А. Математическая модель роговицы глаза с учетом экспоненциальной нелинейности ее упругих свойств при условии геометрической малости деформаций / А.А. Штейн, И.Н. Моисеева, Г.А. Любимов // Российский журнал биомеханики. - 2019. -Т. 23, № 3. - С. 375-390.
3. Andreassen T.T., Simonsen A.H., Oxlund H. Biomechanical properties of keratoconus and normal corneas // Experimental Eye Research. - 1980. - Vol. 31, No. 4. -P. 435-441.
4. Ariza-Gracia M.Á., Zurita J., Piñero D.P., Calvo B., Rodríguez-Matas J. F. Automatized Patient-Specific Methodology for Numerical Determination of Biomechanical Corneal Response // Annals of Biomedical Engineering. - 2015. - Vol. 44, No. 5. - P. 1753-1772.
5. Bekesi N., Dorronsoro C. de la Hoz A., Marcos S. Material Properties from Air Puff Corneal Deformation by Numerical Simulations on Model Corneas // PLoS ONE. - 2016. -Vol. 11, No. 10. D0I:10.1371/journal.pone.0165669.
оптимизации с поиском минимума целевой функции метода наименьших квадратов на основе измеренных данных о амплитуде деформации роговицы пациента, полученных на пневмотонометре Corvis ST. Выбор алгоритма минимизации производится на основе серии предварительных вычислительных
экспериментов по оптимизации параметров моделей здоровых роговиц и роговиц с кератоконусом по наименьшему отклонению расчетных и опытных данных.
4. На основе разработанной методики расчетного восстановления упругих свойств модели материала роговицы были определены значения соответствующих параметров для роговиц в нормальном состоянии и с кератоконусом различных стадий на примере 78 пациентов. Установлены закономерности их изменения в зависимости от степени развития заболевания и различных вариантов описания и положения зоны пониженных свойств. На основе сопоставления расчетных картин распределения геометрических параметров роговицы и интенсивности напряжений и деформаций по ее объему под действием ВГД с картами пахиметрии, элевации и кривизны, полученными при обследовании конкретного пациента, проанализированы основные возможные механизмы локализации области кератоконуса.
6. COMSOL Documentation. Structural Mechanics Module -Large Strain Viscoelasticity. URL: https://d0c.c0ms0l.c0m/5.6/d0cserver/#i/c0m.c0ms0l.help.s me/sme ug theory.06.27.html.
7. Cornaggia A., Boschetti F., Mazzotta C., Pandolfi A. Numerical investigation on epi-off crosslinking effects on porcine corneas // Mechanics of Soft Materials. - 2020. -Vol. 2, No. 5.
8. Dias, J. Anterior and posterior corneal stroma elasticity after corneal collagen crosslinking treatment / J. Dias, V.F. Diakonis, V.P. Kankariya, S.H. Yoo, N.M. Ziebarth // Experimental Eye Research. - 2013. - Vol. 116. -P. 438-451.
9. Eliasy, A., Abass, A., Lopes, B. T., Vinciguerra, R., Zhang, H., Vinciguerra, P. et al. Characterization of cone size and centre in keratoconic corneas // Journal of The Royal Society Interface. - 2020. - Vol. 17, No. 169. DOI: 10.1098/rsif.2020.0271.
10. Elsheikh, A., Anderson, K. Comparative study of corneal strip extensometry and inflation tests // Journal of the Royal Society Interface. - 2005. - Vol. 2, No. 3. - P. 177-185.
11. Gizzi A., De Bellis M. L., Vasta M., Pandolfi A. Diffusion-based degeneration of the collagen reinforcement in the pathologic human cornea // Journal of Engineering Mathematics. - 2021. - Vol. 127, No. 3.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
Grabner, G., Eilmsteiner, R., Steindl, C. et al. Dynamic 27. corneal imaging // Journal of Cataract and Refractive Surgery. - 2005. - Vol. 31, No. 1. - P. 163-174. Hoeltzel, D., Altman, P., Buzard, K., Choe, K. Strip extensiometry for comparison of the mechanical response of bovine, rabbit, and human corneas // Journal Of Biomechanical Engineering. - 1992. - Vol. 114. P. 202-215. Hollman, K. W., Emelianov, S. Y., Neiss, J. H. et al. Strain 28. imaging of corneal tissue with an ultrasound elasticity microscope // Cornea. - 2002. - Vol. 21. - P. 68-73. Holzapfel G.A. Nonlinear solid mechanics : A continuum approach for engineering. - Chichester [etc.]: Wiley, 2000. -455 p. 29.
Huang L., Shen M., Liu T. et al. Inverse solution of corneal material parameters based on non-contact tonometry: A comparative study of different constitutive models // Journal of Biomechanics. - 2020. - No. 112. Jannesari M., Mosaddegh P., Kadkhodaei M. et al. Numerical and clinical investigation on the material model 30. of the cornea in Corvis tonometry tests: differentiation between hyperelasticity and viscoelasticity // Mechanics of Time-Dependent Materials. - 2019. - Vol. 23. - P. 373-384. 31. DOI: 10.1007/s11043-018-9390-3.
Jannesari, M. Numerical and clinical investigation on the material model of the cornea in Corvis tonometry tests: differentiation between hyperelasticity and viscoelasticity / M. Jannesari,-P. Mosaddegh, M. Kadkhodaei,-H. Kasprzak, 32. M. J. Behrouz // Mechanics of Time-Dependent Materials. -2019. - Vol. 23. - P. 373-384.
Jaycock, P.D., Lobo, L., Ibrahim, J., Tyrer, J., Marshall, J. Interferometric technique to measure biomechanical changes 33. in the cornea induced by refractive surgery // Journal of Cataract and Refractive Surgery. - 2005. - Vol. 31(1). - P. 175-184.
Jue, B., Maurice, T. The mechanical properties of the rabbit 34. and human cornea // Journal of Biomechanics. - 1986. -Vol. 19, No. 10. - P. 847-853.
Kazaili, A. Microscale assessment of corneal viscoelastic properties under physiological pressures / A. Kazaili, B. Geraghty, R. Akhtar // Journal of the Mechanical Behavior 35. of Biomedical Materials. - 2019. - Vol. 98 (1). - P. 31-38. Kling S., Bekesi N., Dorronsoro C. et al. Corneal Viscoelastic Properties from Finite-Element Analysis of In Vivo Air-Puff Deformation // PLoS ONE. - 2014. -Vol. 9, No. 8. DOI: 10.1371/journal.pone.0104904. Lago M.A., Rupérez M.J., Martínez-Martínez F. A new 36. methodology for the in vivo estimation of the elastic constants that characterize the patient-specific biomechanical behavior of the human cornea // Journal of Biomechanics. -2015. - No. 48. - P. 38-43.
Liu T., Shen M., Li H. et al. Changes and quantitative 37. 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. 38.
Liu, T. Changes and quantitative characterization of hyper-viscoelastic biomechanical properties for young corneal stroma after standard corneal cross-linking treatment with different ultraviolet-A energies / T. Liu, M. Shen, H. Li et al. 39. // Acta Biomaterialia. - 2020. - Vol. 113. - P. 438-451. Ljubimova, D. Biomechanics of the Human Eye and Intraocular Pressure Measurements // Doctoral Thesis in Mechanics. - Stockholm: Royal Institute of Technology, -2009. - 200 p.
Mahmoud, A.M., Nuñez, M.X., Blanco, C., Koch, D.D., Wang, L., Weikert, M.P., Roberts, C.J. Expanding the Cone Location and Magnitude Index to Include Corneal Thickness and Posterior Surface Information for the Detection of Keratoconus // American Journal of Ophthalmology. -2013. - Vol. 156, No. 6. - P. 1102-1111. DOI: 10.1016/j.ajo.2013.07.018.
Maklad, A. Fluid-Structure Interaction Based Algorithms for IOP and Corneal Material Behavior / O. Maklad, A. Eliasy, K.-J. Chen, J. Wang, A. Abass, B.T. Lopes, V. Theofilis, A. Elsheikh // Frontiers in Bioengineering and Biotechnology. -2020. - No. 8. DOI: 10.3389/fbioe.2020.00970. Mazzotta, C., Moramarco, A., Traversi, C., Baiocchi, S., Iovieno, A., Fontana, L. Accelerated Corneal Collagen Cross-Linking Using Topography-Guided UV-A Energy Emission: Preliminary Clinical and Morphological Outcomes // Journal of Ophthalmology. - 2016. - Vol. 2016. DOI: 10.1155/2016/2031031
Nash, I. Comparison of mechanical properties of keratoconus and normal corneas // Experimental Eye Research. - 1982. - Vol. 35. - P. 413-423. 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. DOI: 10.3389/fbioe.2018.00210.
Nguyen T.D., Boyce B.L. An inverse finite element method for determining the anisotropic properties of the cornea // Biomechanics and Modeling in Mechanobiology. - 2011. -No. 10. - P. 323-337.
Pandolfi A., Fotia G., Manganiello F. Finite element simulations of laser refractive corneal surgery // Engineering with Computers. - 2009. - No. 25. DOI: 10.1007/s00366-008-0102-5.
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. - 12 p.
Qiao, X. Full-field strain mapping for characterization of structure-related variation in corneal biomechanical properties using digital image correlation (DIC) technology / X. Qiao, D. Chen, H. Huo et al. // Medicine in Novel Technology and Devices. - 2021. - No. 11 DOI: 10.1016/j.medntd.2021.100086. Rahmati S.M. Razaghi R., Karimi A. Biomechanics of the keratoconic cornea: Theory, segmentation, pressure distribution, and coupled FE-optimization algorithm // Journal of the Mechanical Behavior of Biomedical Materials. - 2021. - Vol. 113.
Roy A.S. Air-puff associated quantification of non-linear biomechanical properties of the human cornea in vivo / A.S. Roy, M. Kurian, H. Matalia, R. Shetty // Journal of the Mechanical Behavior of Biomedical Materials. - 2015. -Vol. 48. - P. 173-182.
Roy A.S., Dupps W.J. Jr. Patient-specific modeling of corneal refractive surgery outcomes and inverse estimation of elastic property changes // J Biomech Eng. - 2011, Vol. 113, No. 1.
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.
40. Roy A.S., Kurian M., Matalia H., Shetty R. Air-puff 45. 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.
41. Roy, A.S. Inverse computational analysis of in vivo corneal 46. elastic modulus change after collagen crosslinking for keratoconus / A.S. Roy, K.M. Rocha, J.B. Randleman, R.D. Stulting, W.J. Dupps Jr. // Experimental Eye Research. -
2013. - Vol. 113. - P. 92-104. 47.
42. Seiler, T. G., Fischinger, I., Koller, T., Zapp, D., Frueh, B. E., Seiler, T. Customized Corneal Cross-linking: One-Year Results // American Journal of Ophthalmology. - 2016. - 48. Vol. 166. - P. 14-21. DOI: 10.1016/j.ajo.2016.02.029.
43. Seiler, T., Matallana, M., Sendler, S., Bende, T. Does Bowman's layer determine the biomechanical properties of the cornea? // Refractive and Corneal Surgery. - 1992. -
Vol. 8, No. 2. - P. 139-142. 49.
44. Shetty, R., Matalia, H., Srivatsa, P., Ghosh, A., Dupps, W. J., Sinha Roy, A. A Novel Zernike Application to Differentiate Between Three-dimensional Corneal Thickness of Normal Corneas and Corneas With Keratoconus // American Journal 50. of Ophthalmology. - 2015. - Vol. 160, No. 3. - P. 453-462.
DO: 10.1016/j.ajo.2015.06.001.
Shetty, R., Pahuja, N., Roshan, T., Deshmukh, R., Francis, M., Ghosh, A., Sinha Roy, A. Customized Corneal Cross-linking Using Different UVA Beam Profiles // Journal of Refractive Surgery. - 2017. - Vol. 33, No. 10. - P. 676-682. DOI: 10.3928/1081597x-20170621-07. Vellara, H.R. Biomechanical properties of the keratoconic cornea: a review / H.R. Vellara, D.V. Patel // Clinical and Experimental Optometry. - 2015. - No. 100. DOI: 10.1016/j.jmbbm.2019.103375.
Wollensak, G., Spoerl, E. Collagen crosslinking of human and porcine sclera // Journal of Cataract and Refractive Surgery. - 2004. - Vol. 30, No. 3. - P. 689-695. Wollensak, G., Spoerl, E., Seiler, T. Stress-strain measurements of human and porcine corneas after riboflavin-ultraviolet-A-induced cross-linking // Journal of Cataract and Refractive Surgery. - 2003. - Vol. 29, No. 9. -P. 1780-1785.
Woo, S.L., Kobayashi, A.S., Schlegel, W.A., Lawrence, C. Nonlinear material properties of intact cornea and s clera // Experimental Eye Research. - 1972. - Vol. 14, No. 1. - P. 29-39.
Zeng Y., Yang J., Huang K., Lee Z., Lee X.A Comparison of biomechanical properties between human and porcine cornea // Journal of Biomechanics. - 2001. - Vol. 34, No. 4. - P. 533-537.
Финансирование. Исследование выполнено за счет гранта Российского научного фонда (проект № 22-21-20085). Конфликт интересов. Авторы заявляют об отсутствии конфликта интересов.
IDENTIFICATION OF THE MODEL PARAMETERS IN CORNEA
WITH KERATOCONUS IN TERMS OF THE NUMERICAL-EXPERIMENTAL
STUDY OF CORNEAL GEOMETRY AND MECHANICAL BEHAVIOR
E.G. Solodkova1, B.E. Maliugin2, I.N. Zaharov3, V.H. Le3, V.P. Fokin1, S.V. Balalin1, E.V. Lobanov1
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
The aim of this work is to derive a computational technique for model parameter estimation of corneal mechanical behavior basing upon the following input clinical biomarkers: the real-time amplitude of corneal deformations, measured in course of non-contact pneumotonometry, intraocular pressure and geometric features of the patient's cornea, revealed during his tomographic examination.
The developed technique of corneal properties prediction within the area of keratoconus and beyond it, is based upon the data resulting from finite element modeling of a non-contact tonometric test, and upon the procedure of multi-parameter optimization (minimization) of the deviation of the calculated results from the amplitude of corneal deformation measured by means of pneumotonometer Corvis ST. The simulation is based on the results of a survey of 78 patients of a real clinical database, taking into account the geometry of the cornea, the size and position of the keratoconus zone, the intraocular pressure, and the deformation parameters of the material during the pneumotest, which are specific for each patient.
©PNRPU
ARTICLE INFO
Received: 30 March 2023 Approved: 11 September 2023 Accepted for publication: 12 September 2023
Key words:
mathematical modeling, cornea, keratoconus, cross-linking, finite element method, Pentacam, Corvis.
