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

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

а                          б

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

Пространство между передней и задней поверхностями дополнялось до сплошного тела при помощи части эллипсоида, вырезанной двумя этими поверхностями. Длины его полуосей соответствуют диаметрам роговицы в горизонтальной dh^ax и вертикальной плоскостях dV™ (рис. 2, в ). Окончательный вид 3D-модели роговицы приведен на рис. 2, г .

К исходным данным относились полученные на этой стадии величины: максимальный диаметр роговицы d _h ^max , толщины роговицы в апексе t_apex , глубина передней камеры h_c . Фиксировались величины минимальной толщины роговицы t _min , максимальной кривизны K _max передней и задней поверхностей роговицы и координаты соответствующих точек ( X p , y p , х_к , у к ) (табл.1). Остальные геометрические и топографические параметры роговицы (определяемые картами толщины, тангенциальной и сагиттальной кривизны) считались производными от координат и использовались для верификации расчетной геометрии построенной 3D-модели. Также задавалась величина внутриглазного давления p iop по данным тонометрии пациента (с использованием тонометра Corvis ST) .

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

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

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

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

V ( FS ) T + f v = 0 , (1) где V - оператор набла; T - символ транспонирования; f v - вектор объемных сил; S - второй тензор напряжений Пиолы - Кирхгофа; F - градиент деформации, определяемый как

F = I + V u , (2) где I - единичный тензор второго ранга, и - вектор перемещений.

а                            б                            в                            г

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

Таблица 1

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

Глаз	max dh , мм	max d v , мм	hc , мм	t apex , мм	t min , мм	xp , мм	y p , мм	К K max , дптр	xK , мм	yK , мм	p iop Па ,
Правый	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

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

S = S ext

₊ a w

5E

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

E = 1 ( F F - 1 ) .

Тензор напряжений Коши в данной формулировке которой конфигурации роговицы при решении обратной задачи

S • n = q b ,   Г еГ _ь ;

• 3)    действие на передней поверхности Г у роговицы поверхностно распределенной нагрузки интенсивностью q _f с параметрами, соответствующими параметрам воздушного импульса при бесконтактной тонометрии

S • n = q y ,   ГеГ у .           (7)

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

запишется в виде о = J-1FSFT ,

где J = det ( F ) - определитель градиента деформации F , характеризующий изменение объема, вызванное деформацией. Данные уравнения дополняются граничными условиями, заданными на соответствующих поверхностях Γ тела (рис. 4):

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

u = 0, ГеГ S ;

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

б

а

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

и = 0

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

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

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

W s = W iso + W vol + W fb + W is .           (8)

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

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

W iso = C 10 ( 11 - 3 ) + c 20 ( I 1 - 3 ) 2 ,           (9)

где c ₁₀ , c ₂₀ – эмпирические параметры модели; I ₁ – первый инвариант (след) тензора C , а именно

I . = tr ( C ) ,                          (10)

C = J 2 3 C ;

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

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

W_vo/ = - K ( J — 1 ) ²,              (12)

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

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

W fb = E k?- ‘ exP ^k 2 i ( I i* - 3 ) i = 4,6 ^{2 k} 2 i IL

—

- ^ ,

I i *= K i 1 1 + ( 1 — 3 K i ) I i ,             (14)

где k14 , k16 , k24 , k26 – параметры жесткости анизо- тропного материала; Ki e [0,1/3] - параметр структуры (ориентации) волокон, который может быть различным для каждого семейства фибрилл; I4 , I6 – четвертый и шестой псевдоинварианты модифицированного тензора деформаций C из (11), определяемые как [44; 45]

1 4 = § ■ с § , I 6 = £ ■ C Z , (15) где § и Z — единичные векторы, определяющие направление двух наборов армирующих коллагеновых волокон в исходной конфигурации.

Для материала роговицы с двумя характерными группами волокон в (13), (14) можно принять [52] k ₁₄ = k 16 = k 1 и k ₂₄ = k 26 = k 2 , при этом k_1i трактуются [54] как коэффициенты упругости волокон при простом растяжении и имеют размерность напряжений, k 2 i -как безразмерные параметры жесткости фибрилл при больших деформациях. Параметр к i в выражении (14) соответствует степени разориентации волокон вдоль рассматриваемых направлений: к _i = 0 - 100 % коллагеновых фибрилл ориентированы в заданном направлении; к _i = 1/3 - изотропное (неориентированное) распределение фибрилл [44].

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

m w ■ = Ут vis           а , а=1

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

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

О = 2

^а а с '

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

Q a + - Q a = в а S iso ■           (18)

^T a

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

6W-

S =   2---- ^lso ,                      (19)

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

На каждом из уровней сложности модели материала составляющие уравнения (8) вводятся в его структуру по мере необходимости описания особенностей механического поведения роговицы при сопоставлении расчетных и экспериментальных данных, а входящие в выражения (9), (12)-(14), (18) параметры (cw, c20 , К , k1, k2, кi, ва , та ) уточняются из серии вычислительных экспериментов. Предварительные значения указанных коэффициентов принимаются по литературным данным (например, [15; 44; 50; 52; 54; 57]) и приведены в табл. 2.

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

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

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

В этой зоне задается закон распределения эмпирических коэффициентов cw , c20 , устанавливающий их снижение от периферии к центру как по поверхности (по радиусу), так и по толщине роговицы при помощи следующей функции ук (в локальной системе координат Xk , yk , Zk , связанной с центром зоны кератокону-са - точкой Kо на рис. 5; ось Xk - по нормали к задней поверхности роговицы, ось yk - параллельна назальнотемпоральной плоскости, ось Zk - параллельна нижневерхней плоскости):

У k ( x k , y k , ^z k ) = ( 1 ^- у max ) ^exP^[9 k x k R k ⁺

.  2 Г 2  / 2/2 V (20)

где ym,,v - максимальное относительное снижение max жесткости в центре зоны кератоконуса; 9k , 9k , 9k - параметры, устанавливающие градиент изменения свойств по рассматриваемой области (вдоль соответствующей оси); ak , ak , ak - параметры, устанавливающие размер (длину полуоси) зоны кератоконуса вдоль каждой из осей.

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

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

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

C 10 , кПа	C 20 , кПа	К , МПа	k 1 , МПа	k 2	K i	P 1	T , с
140,9	80,13	5,5	0,04	200	0,3329	0,351	410

W K ⁽ x k , yk , z k ) = V k ⁽ x k , y k , z k ^{) +} [ ^C 10 ( 1 1 ^{- 3} ) ⁺ z 2            ,  ⁽²¹⁾

⁺ c 20 ( I 1 - 3 ) ]

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

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

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

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

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

Таким образом, увеличение модулей упругости, вызванное кросслинкингом роговичного коллагена, в любой точке роговицы (с координатами x_c , y_c , z_c , отсчитываемыми от центра зоны облучения, рис. 7, точка C ₀) является функцией радиуса и глубины, которую, по аналогии с функцией (20), примем в виде:

¥_r xx_r , y_r , z_r ) = (1 + Ф exp[9,, x 2/ a 2 + c c c c max c c c

2 . 2 / 2 , * ⁽² 2 ) + ⁹ cp P, + ⁹ c Z Ia2 ]

а

б

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

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

в, г – проекции на плоскость)

где Vmax - максимальное относительное увеличение параметров жесткости ( c10 , c20 ) в результате кросс-линкинга; θc , θc , θc – параметры, устанавливающие cx       y      cz градиент изменения свойств по рассматриваемой области (вдоль соответствующей оси); a , a , a – пара-cx     cy     cz метры, устанавливающие размеры зоны кросслинкинга вдоль каждой из осей (для области облучения круглой формы ac = ac , а также 9c = 9c ).

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

W C = V k ■ V _c - C 10 ( 7 1 - 3 ) + c 20 ( 7 1 - 3 ) 2 . (23)

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

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

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

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

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

qb = -рш - n, ГеГb , (24) при этом величина обратного давления pinv устанавли- вается равной ВГД, скорректированному на разницу площадей задней поверхности роговицы в деформированном и недеформированном состоянии

P inv = ^{- X} P iop , ⁽²⁵⁾ где p iop - величина ВГД, установленная по данным Corvis ST ; X - коэффициент, учитывающий изменение площади поверхности приложения прямого p_iop и обратного давлений p_inv в деформированном и недефор-мированном состоянии роговицы.

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

q b =^- P iop • ⁿ , ^{Г еГ} b . ⁽²⁶⁾

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

б

а

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

на графике ( б ) точками - профиль давления по данным Corvis ST

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

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

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

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

• q f = - P air • n , ^{Г еГ} f .

Пространственно-временная конфигурация давления Pair (x, У, z, t) определяется закономерностями взаимодействия воздушного потока, создаваемого импульсом пневмотонометра, с поверхностью роговицы. Распределение Pair (x, y, z, t) представляется в виде сочетания функции координат поверхности fs (y, z) (так как функция fs задается по передней поверхности Гf роговицы, координата x ее точек предопределена и может быть опущена) и функции времени ft (t) :

P air ( x , У , ^z ,¹ ) = P max ' f s ( У , ^z ) ' f t ( t ) , ^{x еГ} f . ⁽²⁸⁾

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

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

• ^f s ( У , ^z ) = ^exP ⁹ p ( У ² + ^{z 2} )/ ( d air 2 f . ⁽²⁹⁾

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

Временной профиль давления струи приборно фиксируется в ходе пневмотонометрического теста глаза пациента на базе Corvis ST и сохраняется в табличный файл (параметр Pressure Profile) для 140 точек исследуемого интервала времени (рис. 9, б). Полученные экспериментальные результаты используются для формиро- вания функции ft (t) в (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 модельных поверхностей, определяемых встроенными средствами Comsol Multiphysics , в величины кривизны данных поверхностей в меридиональной (тангенциальной) р tng (рис. 12, 13) и сагиттальной р sag

ρ sag

+ z ² ■

плоскостях.

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

Р tng = ( Р 1 ^- Р 2 ) ■

( У ■ n y + z ■ n z )

+ Р 2 ; (30)

ny nx

nz nx

,  (31)

где n x , П у , n z - компоненты вектора n нормали в данной точке поверхности.

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

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

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

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

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

Таблица 4

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

	№ точки
Параметр	1	2	3	4	5	6
Передняя поверхность
Нормальная роговица
Pentacam	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
Кератоконус
Pentacam	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
Кросслинкинг
Pentacam	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
Задняя поверхность
Нормальная роговица
Pentacam	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
Кератоконус
Pentacam	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
Кросслинкинг
Pentacam	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 ) и расчетных (модель) данных о параметрах деформации роговицы в характерных точках

	Роговица
Параметр	нормальная		после кросслинкинга
	DA, мм	PD, мм	DA, мм	PD, мм
Corvis	1,09	4,93	0,87	4,14
Модель	1,06	4,9	0,879	4,347
Отклонение, %	2,75	0,6	1,03	5

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

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

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

Основными контролируемыми параметрами являлись амплитуда деформации роговицы ( DA ), скорость перемещения в апексе ( CV ), расстояние между пиками ( PD ), изменение длины дуги ( dAL ), а также первая и вторая апланации при нагрузке и разгрузке давлением воздушной струи.

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

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

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

Заключение

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

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

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

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

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

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

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

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

Исследование      свойств      двухкомпонентной механической модели глазного яблока и возможности ее использования при практической оценке механических свойств глаза человека // Механика жидкости и газа. – 2014. – № 6. – С. 5–16.

(1–2). – P. 47–57.

properties for young corneal stroma after standard corneal cross-linking treatment with different  ultraviolet-A energies // Acta Biomaterialia. – 2020. – Vol. 113. – P. 438–451.

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