Анализ чувствительности моделей биомеханики глаза человека

Журавлева Дарья Ильинична

Анализ чувствительности моделей биомеханики глаза человека

Биомеханика глаза человека на сегодняшний день является быстро развивающейся областью. За последние несколько десятилетий множество исследований показало, что методы биомеханики предоставляют клинически надёжные результаты, объясняющие механизмы многих офтальмологических заболеваний. Интравитриальные инъекции являются эффективным способом лечения (или частью лечения) различных заболеваний глаза, например, отёков и воспалений внутренних оболочек глаза, диабетического макулярного отёка и возрастной макулярной дегенерации. При интравитреальных инъекциях непосредственно в стекловидное тело вводится лекарственное средство объёмом 0.05 мл, 0.1 мл или 0.2 мл в зависимости от рекомендаций врача. Побочное явление данного метода — резкое повышение внутриглазного давления пациента, что может привести к повреждению глазного нерва и внутренних оболочек глаза даже при непродолжительном воздействии. В данной работе рассмотрена модель корнеосклеральной оболочки глаза для исследования зависимости «объем–давление». Представляя, в первом приближении, глаз в виде сферической или эллипсоидальной оболочки, получено соотношение между изменением внутреннего радиуса, давлением, приложенным к внутренней поверхности оболочки, и изменением внутреннего объема. Проведён глобальный анализ чувствительности трёх моделей изменения внутриглазного давления при интравитреальных инъекциях. На основе предоставленных экспериментальных данных вычислены плотности вероятности случайных величин, характеризующих геометрию оболочки глаза и, затем, проведено сравнение плотности вероятности изменения давления, найденной на основе рассматриваемых моделей, с данными эксперимента. Используя введенный валидационный критерий, вычислены параметры распределения характеристик материала модели (модулей упругости).

Общественные науки в целом

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d363a5f1be77c40d589f7

UUID: c3a3585d-5b9c-4221-b891-7ef98d3ff27f

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 31

Журавлева Дарья Ильинична

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 1930134 bytes

Поделиться работой

Санкт-Петербургский государственный университет Кафедра математического моделирования энергетических систем Журавлева Дарья Ильинична Выпускная квалификационная работа бакалавра Анализ чувствительности математических моделей биомеханики глаза человека CB.5005.2012 Прикладная математика, фундаментальная информатика и программирование Научный руководитель, кандидат физ.-мат. наук, доцент Воронкова Е.Б. Санкт-Петербург 2016

Содержание Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Обзор литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Глава 1. Обзор методов глобального анализа чувствительности . . . 1.1 Метод элементарных эффектов (Elementary Effect Method — EEM) . . . . . . . . . . . . . . . . 1.2 Метод Соболя (Sobol Method) . . . . . . . . . . . . . . . . . . 1.3 Метод разложения в полиномиальный хаос (PCE Method) . . 1.4 Метод ПАВН (PAWN Method) . . . . . . . . . . . . . . . . . . Глава 2. Модели изменения внутриглазного давления при интравитреальных инъекциях . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Моделирование глаза как сферической оболочки . . . . . . . 2.2 Моделирование глаза как эллипсоидальной оболочки . . . . Результаты . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. О методе валидации . . . . . . . . . . . . . . . . . . . . . . . . 2. Обработка экспериментальных данных . . . . . . . . . . . . . 3. Анализ чувствительности . . . . . . . . . . . . . . . . . . . . . 4. Калибровка параметра E1 (E) . . . . . . . . . . . . . . . . . . 5. Калибровка параметра E3 . . . . . . . . . . . . . . . . . . . . 6. Нахождение валидационной ошибки модели (14) . . . . . . . . 7. Анализ чувствительности моделей с откалиброванными параметрами . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . Приложение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 3 5 6 8 8 9 12 15 16 17 20 22 22 23 24 26 27 29 29 32 33 37

Введение Биомеханика глаза человека на сегодняшний день является быстро развивающейся областью. За последние несколько десятилетий множество исследований показало, что методы биомеханики предоставляют клинически надёжные результаты, объясняющие механизмы многих офтальмологических заболеваний, например, невропатию головки зрительного нерва, развивающуюся с возрастом дальнозоркость, катаракту, различные патологии роговицы, отслоение сетчатки и др. [1]. Сотрудничество специалистов разных областей науки (биологии, математики, механики) принесло свои плоды в таких отраслях медицины как лазерная коррекция зрения, контроль развития глаукомы, катаракты и пресбиопии. Интравитриальные инъекции являются эффективным способом лечения (или частью лечения) различных заболеваний глаза [2–5], например, отёки и воспаления внутренних оболочек глаза, диабетический макулярный отёк, возрастная макулярная дегенерация и т. д., и количество применений данной метода лечения растёт с каждым годом [6]. При интравитреальных инъекциях непосредственно в стекловидное тело вводится лекарственное средство объёмом 0.05 мл, 0.1 мл или 0.2 мл в зависимости от рекомендаций врача. Побочное явление данного метода — резкое повышение внутриглазного давления пациента, что может привести к повреждению глазного нерва и внутренних оболочек глаза даже при непродолжительном воздействии. В случаях, когда внутриглазное давление пациента ещё до введения инъекции высокое, врачи рекомендуют уменьшить дозу инъекционного препарата. С помощью методов биомеханики, например, конечно-элементного моделирования, можно разработать модели изменения внутриглазного давления при интравитреальных инъекциях, но параметры таких моделей, в связи с природой своего происхождения, не всегда могут быть вычислены точно (как, например, модуль упругости склеры, [7]), поэтому возникает необходимость качественного изучения таких моделей, например, при помощи валидации модели и последующего анализа чувствительности. Анализ чувствительности — качественный анализ модели и её параметров, помогающий ответить на вопрос: какие параметры оказывают наи3

большее влияние на результат модели. Разделяют локальный и глобальный анализ чувствительности. Локальный исследует поведение модели только в одной точке области изменения параметров, в то время как глобальный охватывает всю данную область, и именно поэтому в данной работе рассматриваются методы глобального анализа чувствительности. Применяя методы анализа чувствительности, можно фиксировать невлиятельные параметры и проводить более точные численные эксперименты. Валидация — процесс, посредством которого результаты расчётов на основе модели сравниваются с экспериментальными данными для оценки ошибки моделирования. Он состоит из трёх этапов: 1) численные эксперименты на основе полученной модели 2) выбор метрики валидации 3) проведение оценки по выбранной метрике. Главная задача исследователя — верно выбрать метрику валидации. Комбинация методов валидации и анализа чувствительности позволяет получать больше информации о параметрах, которые сложно или невозможно вычислить точно. 4

Постановка задачи I Провести анализ чувствительности разработанных ранее [8, 9] математических моделей изменения внутриглазного давления. II На основе предоставленных экспериментальных данных вычислить плотности вероятности случайных величин, характеризующих геометрию оболочки глаза. III Используя введенный валидационный критерий, вычислить параметры распределения характеристик материала модели (модулей упругости). IV Провести сравнение плотностей вероятности изменения давления, найденных на основе рассматриваемых моделей, с данными эксперимента. 5

Обзор литературы При качественном анализе модели один из главных вопросов — анализ чувствительности, выявление влиятельных и невлиятельных параметров. Для моделей биомеханики этот вопрос ещё более актуален, т. к. параметры таких моделей имеют широкий диапазон вариации вследствие индивидуальных особенностей или патологий пациентов, например, в моделях биомеханики глаза параметры, описывающие геометрические размеры, можно определить вживую у пациента, однако, свойства оболочек неизмеримы таким образом, а данные литературы разнятся. Одними из первых анализ чувствительности в биомеханику глаза человека ввели Sigal I. A. и др. [10] в 2005 году. Предметом их изучения была биомеханика головки зрительного нерва. Они выделили 21 параметр, и, систематически варьируя каждый параметр, определяли степень его влияния на отклик модели. По результатам исследования, наибольшими показателями чувствительности обладали: модуль упругости склеры, радиус глаза, модуль упругости решётчатой пластинки, внутриглазное давление и толщина склеральной оболочки. В 2011 году такой же метод анализа чувствительности был применён Sigal. I. A. и др. [11] к модели деформации решётчатой пластинки и модели расширения склерального канала у обезьян, в результате повышения внутриглазного давления. Подводя итоги исследования, наибольшими показателями чувствительности обладали: жёсткость и расположение решётчатой пластинки, в первой модели, и жёсткость и толщина склеры, во второй модели. В том же 2011 году Nguyen T. D. и Boyce B. L. [12] исследовали модель деформации роговицы свиней, и её чувствительность к анизотропным характеристики роговицы. По результатам численных экспериментов, наибольшим влиянием на отклик модели обладали: матрица модулей сдвига, параметр упругости коллагеновых волокон и степень анизотропии лимба роговицы. Результаты были получены методом систематического варьирования параметров. Sigal I. A. и др. в 2014 году изучали характеристики тканей головки зрительного нерва [13]. Авторы предположили, что они могут быть числен6

но оценены через параметры, которые можно вживую померить у человека. В доказательство этой теории авторы вывели статистические модели, прогнозирующие характеристики тканей, и исследовали все модели, описывающие биомеханику головки зрительного нерва. Далее их интересовал вопрос чувствительности статистических моделей. В качестве метода анализа чувствительности авторы выбрали: матрицу корреляции параметров и систематическое варьирование. Наибольшими показателями чувствительности обладали: размер глаза и шлеммов канал, параметры, представляющие влияние повышения внутриглазного давления на склеру, а также динамический эффект внутриглазного давления. Saltelli A. и др. [14] показали, что такой принцип исследования чувствительности параметров, локальный, т.е. исследование влияние параметра, основываясь на его изменениях при фиксации всех остальных параметров, не является надёжным и может давать ошибочные результаты, поэтому необходимо применять методы глобального анализа чувствительности. Подробный обзор таких методов приведён в Главе 1. 7

Глава 1. Обзор методов глобального анализа чувствительности Задачи глобального анализа чувствительности можно сформулировать следующим образом: деление параметров на группы (factor screening), ранжирование параметров (factor ranking) и фиксирование параметров (factor fixing). 1.1 Метод элементарных эффектов (Elementary Effect Method — EEM) Метод элементарных эффектов был предложен Morris M. D. [15] в 1991 году. Данный метод предназначен для деления параметров на группы (для так называемой задачи factor screening). Его основным назначением служит получение информации о чувствительности модели с количеством параметров большим 100. Рассмотрим математическую модель с k независимыми параметрами, заданную в виде функции y = f (x1 , . . . , xk ), y ∈ R, xi ∈ [0, 1] = I, i = 1, k, x ∈ Ω = Ik . (1) Следуя [15], введём регулярную k-мерную p-уровневую сетку ω = {0, 1/(p − 1), 2/(p − 1), ..., 1} и множитель ∆, кратный 1/(p − 1). Элементарным эффектом i-ого параметра в точке X ∈ ω называется: EEi (X) = f (X + ei ∆) − f (X) , ei = {0, ...0, 1, 0, ..., 0}, i ∆ | {z } k причём X выбирается так, чтобы (X + ei ∆) ∈ ω. Таким образом, для каждого i-ого параметра можем посчитать pk−1 (p−∆(p−1)) элементарных эффектов. Положим, что для каждого параметра вычислено r элементарных эффектов. Для каждого такого набора элементарных эффектов можем посчитать выборочное математическое ожидание µi и выборочную дисперсию σi . Тогда деление параметров на группы происходит следующим образом: 8

1. параметры с незначительным эффектом на отклик модели; им соответствует малое значение µi и малое значение σi , 2. параметры с линейным эффектом; им соответствует большое значение µi и малое значение σi , 3. с нелинейным эффектом; им соответствует большое значение σi . Таким образом, для разделения k параметров на группы необходимо совершить N = 2rk вычислений. Уменьшить количество вычислений можно, если выбрать p чётным, а ∆ = p/(2(p − 1)). Кроме того можно воспользоваться предложенным в [15] траекторным алгоритмом генерации точек X (trajectory-based) (см. Приложение). Результаты исходного алгоритма ненадёжны в случаях немонотонных функций [16], поэтому в 2007 году Campolongo F. и др. предложили вместо µi использовать µ?i , абсолютное выборочное математическое ожидание, посчитанное для выборки |EEi | размером r. В этой же статье авторы предложили усовершенствование алгоритма генерации точек X, однако такое улучшение оказалось не оптимальным по времени. Оптимальный, равномерно-траекторный, алгоритм генерации точек X, предложили Khare Y. P. и др. в 2015 году [17] ( Sampling for Uniformity). (см. Приложение) 1.2 Метод Соболя (Sobol Method) На сегодняшний день является наиболее широко используемым методом. С помощью него могут быть решены все задачи анализа чувствительности: деление параметров на группы (factor screening), ранжирование параметров (factor ranking) и фиксирование параметров (factor fixing). Метод был представлен профессором, доктором физико-математических наук, Ильёй Мееровичем Соболем, российским и советским математиком в 2001 году. [18] Рассмотрим модель в виде (1), будем полагать, что f (x) интегрируемая c квадратом в Ω функция. ANOVA-разложением функции f (x) 9

называется следующее представление: f = f0 + k X X fi1 ...is (xi1 , . . . , xis ), если s=1 i1 <...<is Z 1 (2) fi1 ...is (xi1 , . . . , xis )dxi1 . . . dxis = 0. 0 ANOVA-разложение единственно [18]. Из соотношения (2) следует алгоритм вычисления компонентов разложения: Z 1 f (x)dx = f0 0 Z 1 f (x) 0 Z dxp = f0 + fi (xi ) p6=i 1 f (x) 0 Y Y dxp = f0 + fi (xi ) + fj (xj ) + fij (xi , xj ) . . . p6=i,j . . . fi1 ...ik = f (x) − f (0) − k X fi (xi ) − . . . . i=1 Пример. ANOVA-разложение f = xR + xy 1 f0 = 0 (x + xy)dxdy = 3/4 R1 R1 f1 = 0 (x + xy)dy − f0 = 23 x − 34 , f2 = 0 (x + xy)dx − f0 = 21 y − f12 = f (x) − f0 − f1 − f2 = xy − 21 y − 21 x + 41 Z 1 4 1 Следующие константы назовём дисперсиями: D = f 2 (x)dx − f0 , 0 Z 1 Di1 ...is = fi21 ...is dxi1 . . . dxis . Выделим набор параметров xZ , Z = {i1 , . . . , im }. 0 Тогда xT = x \ xZ . Главным (main) (SxZ ) и полным (total) (SxTZOT ) эффектами или 10

показателями Соболя набора параметров xZ называются величины [18] m P SxZ = P s=1 (i1 <...<is )∈Z SxZ , SxTZOT D ∈ [0, 1] Di1 ...is SxTZOT , DxTZOT D − DxT = = . D D Пример. Главные и полные показатели Соболя. x = (x1 , x2 , x3 ), x = (xZ , xT ) xZ = x1 , xT = (x2 , x3 ) Sx1 = Sx1 , SxT1OT = Sx1 + Sx1 ,x2 + Sx1 ,x3 + Sx1 ,x2 ,x3 = 1 − Sx2 ,x3 xZ = (x1 , x2 ), xT = x3 Sx1 ,x2 = Sx1 + Sx2 + Sx1 ,x2 SxT1OT ,x2 = Sx1 + Sx2 + Sx1 ,x2 + Sx1 ,x3 + Sx2 ,x3 + Sx1 ,x2 ,x3 = 1 − Sx3 Полный показатель чувствительности, в отличие от главного, учитывает взаимодействие параметров. Главные показатели используются для задачи ранжирования параметров (factor ranking). Полные показатели используются для задач деления на группы и фиксирования параметров (factor screening & factor fixing): если полный показатель чувствительности параметра близок к нулю (< 10−3 ), то эффект параметра на результат модели можно считать незначительным и принять значение параметра за константу, если же полный показатель параметра близок к единице, то можно сделать вывод, что результат модели наибольшим образом зависит от этого параметра. На практике, однако, не все модели представимы в виде аналитических функций, как это бывает, например, при конечно-элементном моделировании в задачах биомеханики, поэтому в 2010 году Saltelli A. и др. [19] предложили вычислительный эквивалент показателей Соболя. Сначала генерируются независимо друг от друга матрицы точек ∈ Ik : A = [N × k] и B = [N × k], и формируются матрицы Ci = [A1 ...Ai−1 B i Ai+1 ...Ak ], i = 1, k, где Ai — i-ый столбец матрицы . Главный и полный эффект i-ого параметра вычисляются: N N 2 1 X j 1 X j j j j T OT Si = y y − yA , S i =1− y − yC i , DN j=1 B Ci 2DN j=1 A 11

N P !2 N P 1 1 yAj − yAj , а yA = y(A), yB = y(B) и yCi = y(Ci ) N j=1 N j=1 векторы [N × 1]. Рекомендуется выбирать N в пределах 102 < N < 105 [19]. где D = 2 1.3 Метод разложения в полиномиальный хаос (PCE Method) Не для всех моделей можно посчитать аналитически показатели чувствительности по методу Соболя, поэтому Sudret B. [20] в 2008 году предложил использовать разложение модели в полиномиальный хаос в качестве приближения модели, и вычислять показатели Соболя аналитически, но для приближения. В модели (1) будем полагать: x ∈ Ω ∈ Rk . Разложением в полиномиальный хаос будем называть: y = f (x) = X aα Ψα (x) = aT Ψ(x), α∈N k ye ≈ X aα Ψα (x) α∈A Ψα (x) = k Y φαi i (xi ), α = [α1 , . . . , αk ] – мультииндекс. i=1 φαi i - полином степени αi , тип которого выбирается в соответствие с плотностью вероятности параметра xi по схеме Аски (см. таблицу 1) [21]. Так же требуется, чтобы параметры xi , i = 1, k были независимы и соответствовали одному и тому же закону распределения, если это не так, то ищется преобразование ξi = T (xi ), i = 1, k, в котором ξi уже соответствуют требованию. Здесь и далее будем полагать, что условие выполнено для xi , i = 1, k Полагается, что все параметры имеют одну и ту же плотность вероятности. Если же это не так, то найдётся преобразование [22], приводящее к новым параметрам xbi , i = 1, k, которые соответствуют данному требованию. Далее будем считать, что утверждение верно для xi , i = 1, k. Существует множество алгоритмов, определяющих количество и значения коэффициентов aα .Часто применяется алгоритм LARS, Adaptive 12

Таблица 1: Примеры соответствий между семейством полиномов и законом распределения случайной величины ξi . Сем-во пол-ов φ Полиномы Хана Полиномы Якоби Полиномы Мейкснера Полиномы Кравчука Полиномы Лагерра Полиномы Шарлье Полиномы Эрмита Полиномы Лежандра Распределение xi Гипергеометрическое Бета-распределение Отрицательное биномиальное Биномиальное Гамма-распределение Распределение Пуассона Стандартное нормальное Равномерное Пар-ры р-ия HG(D, N, n) Be(α, β) N B(r, p) B(n, p) Γ(k, θ) P (λ) N (0, 1) U (a, b) Sparse Least Angle Regression, предложенный в 2011 году Blatman G. и Sudret B. [22]. В каждом алгоритме, строящем разложение в полиномиальный хаос, есть три принципиальных вопроса: ? выбор мультииндексов α, ? вычисление коэффициентов разложения a, ? вычисление погрешности (Err) приближения ye к y. В алгоритме LARS предлагается решать эти вопросы следующим образом k A = Ak,p | ||α||q 6 p}, q = {α ∈ N !1/q k X ||α||q = αiq , 0 < q 6 1, p ∈ N. i=1 Введём следующее обозначение: f P C (x) = что имеем N результатов работы модели: X = {x(1) , . . . , x(N ) }, Y = {y (1) , . . . , y (N ) }. a= !−1 N X 1 Ψ(x(i) )ΨT (x(i) ) N i=1 13 P a∈A aα Ψα (x). Положим, N 1 X Ψ(x(i) )y (i) N i=1 ! (3)

(4) Err = E ∗ T (N, p), 2 (i) PN y − f P C (x(i) ) i=1 1 − hi , E= y (i) − µ by −1 T T h = diag Q Q Q Q , Qij = Ψj (x(i) ), j = 0, |A|, i = 1, N , N 1 X (i) µ by = y , N i=1 N N −p T (P, N ) = Cemp = −1 tr(Cemp ) 1+ N ! , 1 T A A, tr − след матрицы. N Выбор составляющих базиса (т.е. Ψα ) таких, которые имеют наибольшее влияние на результат модели, обеспечивает часть алгоритма LARS – функция LAR, least angle regression. Далее будем полагать: Ψ = Ψ?min – полученный из функции LAR min – соответствующий набор мультииндексов. Тогда, слебазис, A = Ak,p q дуя [20], можем аналитически вычислить показатели чувствительности. Введём следующее обозначение: Li1 ,...,is = {α ∈ A | αj = 0 ∀j 6= i1 , . . . , is }. Перепишем полиномиальное разложение в следующем виде: f PC (x) = a0 + k X X aα Ψα (xi )+ i=1 α∈Li + X aα Ψα (xi1 , xi2 ) + . . . 16i1 <i2 6k ... + X aα Ψα (x1 , . . . , xk ). α∈Li1 ,...,ik Пользуясь тем свойством, что ANOVA-разложение единственно [18], и учитывая свойство ортонормальности базиса [21], показатели чувствительности могут быть вычислены: 14

yα2 P D= X a2α , Si1 ,...,is = α∈Λi1 ,...,is D α∈A , OT SiT1 ,...,i =1− s Sk\{i1 ,...,is } . D 1.4 Метод ПАВН (PAWN Method) До сих пор полагалось, что дисперсия наилучшим образом отражает чувствительность модели. Однако, не для всех моделей это верно, например, в моделях с сильно асимметричной или многомодальной плотностью вероятности [23] дисперсия не может полностью характеризовать их чувствительность к параметрам. Возникает необходимость в методе, который бы не зависел от моментов и исследовал бы функцию распределения модели целиком. В 2015 году Pianosi F. и Wagener T. [24] и предложили метод ПАВН, в рамках которого в модели (1) будем полагать: x ∈ Ω ∈ Rk . PAWN 1: Определим Nu - количество точек, необходимое чтобы построить выборочную функцию распределения, Nc - количество точек, необходимое, чтобы построить условную выборочную плотность распределения, и n - количество условных значений. 2: Генерируются Nu точек, строится выборочная функция распределения Fy (y). 3: Для кажого параметра xi , i = 1, k выбираются n условных значений. Для каждого такого значения определяются Nc точек ∈ Rk−1 , строится выборочная условная функция распределения Fy|jx , j = 1, n, вычисляi ется статистика Колмогорова-Смирнова: KSij = max|Fy (y) − Fy|jx | y 4: i Для каждого параметра считается показатель чувствительности: Ti = stat KSij , {1,...,j,...,n} stat - либо выборочная медиана, либо выборочное математическое ожидание, либо максимум. В исходной статье [24] авторы использовали максимум. 15

Глава 2. Модели изменения внутриглазного давления при интравитреальных инъекциях Рис. 1: Строение глаза человека и роговицы. Видоизменено с https://resources.emaze.com/vbscenes/layoutimages/635127658825199996_eye. png, http://www.garyfostermd.com/wp-content/uploads/2013/11/17754256_ml.jpg. Исследование, проведённое в 2007 году Котляром К. Е. и др. [2], показало, что интравитреальные инъекции (а, следовательно, и изменение объёма стекловидного тела) оказывают значительный эффект на внутриглазное давление. В 2015 году Бауэр С. М. и Воронкова Е. Б. представили подробный обзор неклассических теорий оболочек в биомеханике глаза человека [8], позволяющих описывать различные виды деформации. Глаз человека – сложная биологическая структура, поэтому для построения моделей авторами [8] были сделаны некоторые предположения. Склера составляет бóльшую часть оболочки глаза, поэтому будем полагать, что оболочка глаза неделима и имеет толщину и характеристики склеры. Строма составляет более 90% всей толщины роговицы, поэтому положим, что оболочка — один слой. Будем считать, что на этот слой действует только внутриглазное давление, а снаружи нет никаких воздействий. Наконец, будем полагать, что внутри оболочки находится однородная несжимаемая жидкость [9]. Если пациент не страдает миопией или гиперопией, то его глаз по форме близок к шару. В случаях, когда пациент страдает такими патологиями, форма глаза близка к эллипсоиду. [8, 9]. Самая простая модель механических свойств слоя — изотропный слой, 16

но такое предположение плохо согласуется с клиническими результатами [9], поэтому будем полагать, что слой является анизотропным, а в первом приближении — трансверсально-изотропным (т.е. модуль упругости в направлении толщины на два порядка меньше модуля упругости в тангенциальном направлении). 2.1 Моделирование глаза как сферической оболочки В [8] глазное яблоко моделируется трансверсально-изотропной однослойной оболочкой, заполненной однородной несжимаемой жидкостью. Рис. 2: Модель трансверсально-изотропного сферического слоя Следуя [8], рассмотрим однослойную сферическую оболочку толщины h = R1 − R2 , где R1 — внутренний радиус оболочки, R2 - внешний; (ρ, ϑ, φ) — сферическая система координат (рис. 2). На внутреннюю и внешнюю поверхность слоя действуют давления p1 и p2 , соответственно. Введём следующие обозначения: x1 = ϑ, x2 = φ, x3 = ρ. Пусть справедливы закон Гука и условия симметрии σ11 ν12 ν13 − σ22 − σ33 , E1 E1 E1 ν21 σ22 ν23 − σ33 , e22 = − σ11 + E2 E2 E2 ν31 ν32 σ33 e33 = − σ11 − σ22 + , E3 E3 E3 νji νij = i 6= j, i, j = 1, 3, Ei Ej e11 = σ12 , G12 σ23 = , G23 σ13 = , G13 e12 = e23 e13 (5) где e – тензор деформаций, σ – тензор напряжений, Ei – модуль Юнга в i-ом направлении, Gij – модуль сдвига в (i − j) направлении, νij – коэффициент Пуассона при растяжении в i-ом направлении и сжатии в j-ом направлении. Полагая в каждой точке (x1 − x2 ) поверхностью изотропии, т.е. E1 = E2 , ν13 = ν23 , перепишем (5): 17

(6) σ11 = E11 e11 + E12 e22 + E13 e33 , σ22 = E12 e11 + E11 e22 + E13 e33 , σ33 = E13 e11 + E13 e22 + E33 e33 , σ12 = G12 e12 , G13 = G23 , E11 E13 σ23 = G23 e23 , σ13 = G13 e13 , E1 , 2(1 + ν12 ) 2 ) E1 (E1 ν12 + E3 ν13 = 2 ), (1 + ν12 )(E1 (1 − ν12 ) − 2E3 ν13 E1 E3 (1 − ν12 ) = 2 . E1 (1 − ν12 ) − 2E3 ν13 G12 = 3 ) E1 (E1 − E3 ν13 = 2 ) , E12 (1 + ν12 )(E1 (1 − ν12 ) − 2E3 ν13 E1 E3 ν13 = E33 2 , E1 (1 − ν12 ) − 2E3 ν13 Будем считать, что после инъекции глаз сохраняет форму шара. В силу этого предположения и симметрии, деформация (перемещение) происходит только вдоль x3 = ρ, и единственная из трёх отличная от нуля компонента вектора перемещений – радиальная, u3 (ρ). Тогда из трёх уравнений равновесия остаётся одно: dσ33 2 + (σ33 − σ22 ) = 0, dρ ρ σ33 (Ri ) = −pi , (7) i = 1, 2. Компоненты тензора деформации примут следующий вид: e11 = e22 = u3 (ρ) , ρ e33 = du3 (ρ) , dρ e12 = e13 = e23 = 0. (8) Из (6) можем получить компоненты тензора напряжений: σ11 = σ22 = (E11 + E12 )e22 + E13 e33 , σ12 = σ13 = σ23 = 0. 18 σ33 = 2E13 e22 + E33 e33 , (9)

Таким образом, подставляя (8) и (9) в (7), получаем: du3 (ρ) 2 du3 (ρ) u3 (ρ) + − M = 0, dρ2 ρ dρ ρ2 u3 (Ri ) du3 (ρ) σ33 (Ri ) = 2E13 + E33 = −pi , Ri dρ Ri E11 + E12 − E13 . M =2 E33 (10) Решая дифференциальное уравнение (10), получаем радиальную компоненту перемещения как функцию от ρ: A1 m A2 −(m+1) ρ + ρ , µ1 µ2 √ 1 − 1 + 4M m=− , 2 µ1 = mE33 + 2E13 , u3 (ρ) = где (11) µ2 = 2E13 − (m + 1)E33 , p1 R1m+2 − p2 R2m+2 A1 = − 2m+1 , R1 − R22m+1 m+2 A2 = (R1 R2 ) p1 R2m−1 − p2 R1m−1 . R22m+1 − R12m+1 Согласно сделанным ранее предположениям, p1 = −∆P, p2 = 0. Таким образом, перемещение вдоль радиальной компоненты может быть найдено по формуле (11): u(R1 ) = ∆P C1 R1m + −(m+1) C2 R1 = ∆P ∗ u e(R1 , R2 , E1 , E3 , ν12 , ν13 ), где −R1m+2 , C1 = R12m+1 − R22m+1 (mE33 + 2E13 ) R1m+2 R22m+1 C2 = (2E13 − (m + 1)E33 ) R12m+1 − R22m+1 В рассматриваемой модели после деформации слой сохраняет сферическую форму. Это позволяет вычислить изменение внутреннего радиуса оболочки, зная изменение объёма: ∆R = 3 4π 4π 3 R + ∆V 3 1 1/3 − R1 , ∆V – объём инъекции. 19 где

Тогда зависимость между изменением внутриглазного давления ∆P и изменением внутреннего объёма ∆V : ∆P = ∆R1 u e(R1 , R2 , E1 , E3 , ν12 , ν13 ) (12) В случае изотропного слоя выражения (12) упрощается. M = 2, m = 1, E1 = E3 = E, ν12 = ν13 = ν, E −2E µ1 = , µ2 = 1 − 2ν 1+ν 3 R13 R23 (1 + ν) −R1 (1 − 2ν) , C2 = , C1 = (R13 − R23 ) E −2E (R13 − R23 ) C2 u e(R1 , R2 , E, ν) = C1 R1 + 2 , R1 ∆P = ∆R1 u e(R1 , R2 , E, ν) (13) 2.2 Моделирование глаза как эллипсоидальной оболочки Для глаз с нормальным зрением характерна форма, близкая к сферической. Близорукие и дальнозоркие глаза, обычно, имеют форму вытянутого или сплюстнутого эллипсоида, соответственно [2]. Рассмотрим эллипсоид, образованный при вращении эллипса с полуосями Ra , Rb вокруг вертикальной полуоси Rb . Объем эллипсоида равен 4 V0 = π(Ra + ua )2 (Rb + ub ). Пусть в результате деформации длина верти3 кальной оси изменилась на ub = ∆P ũb , а горизонтальной на ua = ∆P ũa , и после деформации тело также является эллипсоидом. Тогда изменение объем тела в результате деформации, например, при изменении внутреннего давления (∆P ) будет равно 4 4 ∆V = V − V0 = π(Ra + ∆P ũa )2 (Rb + ∆P ũb ) − πRa2 Rb . 3 3 Решая обратную задачу, зная изменение внутреннего объема ∆V = V − V0 , 20

можно найти зависимость ∆P (∆V ) ∆p = Φ(∆V, V0 , E, ν), где E, ν — параметры, характеризующие свойства материала оболочки. Для определения перемещений ua , ub воспользуемся теорией анизотропных оболочек, разработанной Родионовой В.А., Титаевым Б.Ф., Черныхом К.Ф. [25]. Эта теория позволяет учесть при построении модели деформации анизотропных оболочек учитывать поперечные сдвиги, а также поперечные нормальные напряжения, повороты волокон и изменения их длины в процессе деформации. При построении теории полагается, что касательные напряжения распределены по толщине по закону полинома второй сте- пени, а нормальное напряжение — третьей степени. Тангенциальные и нормальная составляющие вектора перемешения распределены по толщине оболочки по кубической и квадратичной параболы соответсвенно [25]. В [8] проведено сравнение решений, полученны в рамках трехмерной теории упругости, с решениями на основе теорий анизотропных оболочек и сделан вывод о применимости уточненных теорий оболочек к задачам о дефорамции оболочек эллипсоидальной формы. Для нахождения определяюших соотношений для ua и ub можно воспользоваться, частично, результатами работы [26]. Таким образом, в предположении, что материал оболочки подчиняется закону (5), разрешающая система уравнений для нахождения ∆P (∆V ) имеет вид 4 4 ∆V = π(Ra + ua )2 (Rb + ub ) − πRa2 Rb , 3 3 2 h R0 (1 − ν12 ) R0 ν13 + + , (14) ua = ∆p 2E1 h 2E1 4E33 2 2 R0 ν13 (2 − 1/κ2 ) h 2 R0 (2 − ν12 − 1/κ ) ub = ∆p κ +κ + , 2E1 h 2E1 4E33 R0 = Ra2 /Rb , κ = Rb /Ra . Здесь Ra , Rb — радиусы полуосей срединной поверхности эллипсоидальной оболочки, моделирующей глазное яблоко. 21

Результаты 1. О методе валидации Валидация – это процесс, помогающий исследователю оценить, насколько хорошо представленная модель описывает исследуемое явление. Этот процесс состоит из трёх этапов [27]: 1. численные эксперименты на основе полученной модели, 2. выбор метрики валидации, 3. проведение оценки по выбранной метрике. Главная задача исследователя — верно выбрать метрику. Метрика валидации — это математическая мера соответствия между численными прогнозами и экспериментальными данными. Различают • детерминистические метрики: сравнение результатов происходит на основе их графического представления, • метрики экспериментальной чувствительности: оценивается точность измерения входных и выходных параметров. • метрики численной ошибки представляют собой оценку стабильности численных вычислений (например, устойчивость к выбору метода интегрирования), • недетерминистические: метрика выбирается на основе сравнения статистических показателей (например, выборочного математического ожидания или выборочной функции распределения) численных расчётов и экспериментальных данных. Один из недетерминистических методов, представленный в [28], используется в данной работе. Будем полагать, что имеется выборка реальных экспериментальных данных результатов процесса, который описывает исследуемая модель. Для этих данных можно вычислить: Eµ — выборочное математическое ожидание и EP10 , EP90 — два персентиля. Процесс валидации проходит следующим образом: 22

Валидация 1: Делается предположение о функциях распределения параметров модели и об области определения характеристик θ этих функций: θ ∈ Θ (например, функция распределения, соответствующая нормальному закону распределения, имеет две характеристики: математическое ожидание и среднеквадратичное отклонение). 2 2 2 1/2 2: θ = argmin (fµ (x|θ) − Eµ ) + (fP10 (x|θ) − EP10 ) + (fP90 (x|θ) − EP90 ) θ∈Θ где fµ (x|θ), fP10 (x|θ), fP90 (x|θ) – выборочное математическое ожидание и два персентиля выборки выхода модели с характеристиками функций распределения параметров равными θ. Следующее выражение будем называть валидационной ошибкой (Error): Error = (fµ (x|θ) − Eµ )2 + (fP10 (x|θ) − EP10 )2 + (fP90 (x|θ) − EP90 )2 1/2 (15) 2. Обработка экспериментальных данных Университет прикладных наук города Аахен, Германия (Aachen University of Applied Sciences) совместно с Университетской глазной клиной (University Eye Clinic RWTH Aachen University) проводили измерения ряда геометрических параметров глаза, а также значения внутриглазного давления до и после введения инъекции [29]. На основе полученных данных в таблице 9 (см. Приложение), были найдены (см. Приложение, код MATLAB) плотности вероятности параметров h (толщина склеры) и R1 (внутренний радиус оболочки). Было рассмотрено четыре закона распределения: бетараспределение (fB ), гамма-распределение (fG ), логнормальное распределение (fL N ) и нормальное распределение (fN ). Ниже приведены плотности вероятности этих распределений в общем виде. fB (x) = 1 xα−1 (1 − x)β−1 , B(α, β)  −x/θ  xk−1 e , x>0 θk Γ (k) fG (x) =  0, x<0 23 α > 0, β > 0, k > 0, θ > 0, x ∈ [0, 1] x ∈ [0, +∞)

1 √ e−(ln(x)−µ) 2 /2σ 2 , σ > 0, µ ∈ R, x ∈ (0, +∞) xσ 2π 1 2 2 fN (x) = √ e−(x−µ) /2σ , σ > 0, µ ∈ R, x ∈ R σ 2π fLN (x) = Параметры плотностей вероятности находились с помощью комбинации встроенных функций MATLAB fitdist + ksdensity + pdf (код приведён в Приложении). Затем вычислялась среднеквадратичная ошибка по следующей формуле qX Errdistr = (Distr(xi ) − f (xi ))2 , где Distr(xi ) — значение плотности вероятности предполагаемого закона распределения в xi , f (xi ) — значение плотности вероятности, построенной по экспериментальным данным, в xi . Значения ошибки для четырёх законов распределений представлены в табл. 2. Логнормальное распределение ближе других к плотностям вероятностей параметров: h соответствует функция распределения LN (−7.4910, 0.0545), а R1 — LN (−4.4714, 0.0513). Графики сравнения плотностей вероятности h и R1 и плотностей вероятности четырёх распределений изображены на рис. 3. Таблица 2: Среднеквадратичная ошибка аппроксимации плотностей вероятности параметров h и R1 . h R1 Beta Gamma Lognormal Normal 8899.9518 8899.9786 8452.8076 8932.3892 737.0450 736.0295 712.4760 837.8101 3. Анализ чувствительности Границы изменения и единицы измерения параметров моделей (13), (12) и (14) указаны в табл. 3. Положим, что параметры равномерно распределены в указанных интервалах. Построим плотности вероятности соответствующих моделей. Они не являются ни сильно асимметричными, ни многомодальными (рис. 4), следовательно, нет необходимости применять метод PAWN для анализа чувствительности. В этих моделях количество параметров не большое (меньше десяти), таким образом, можно предварительно не разделять их на группы влиятельных и нет. Метод разложения в 24

Рис. 3: Сравнение плотностей вероятности толщины h и внутреннего радиуса R1 и плотностей вероятности четырёх распределений. Синим цветом изображены плотности вероятности параметров, полученные по данным эксперимента, красным – плотности вероятности четырёх рассматриваемых распределений полиномиальный хаос является более расчётно-устойчивым по сравнению с методом Соболя, так как аналитически вычисляет показатели чувствительности, поэтому именно этот метод будем использовать в дальнейшем. Таблица 3: Границы изменения и единицы измерения параметров ν12 (ν) ν13 Ra (R1 ) Rb E1 (E) E3 h ∆V Н. граница 0.4 0.01 10 10 3 0.03 0.5 0.05 В. граница 0.49 0.1 14 14 10 0.1 1 0.2 Ед. изм. мм мм МПа МПа мм мл В предположении, что все параметры равномерно распределены в указанных в таблице 3 пределах, и проведя анализ чувствительности по методу PCE (рис. 5). Параметры ν12 (ν) и ν13 обладают очень малыми показателями чувствительности. Следовательно, в дальнейших расчётах они могут быть зафиксированы. Результаты представлены в таблице 4. 25

Рис. 4: Плотности вероятностей откликов моделей (13), (12) и (14) в предположении, что все параметры равномерно распределены в указанных в таблице 3 интервалах. Таблица 4: Результаты анализа чувствительности моделей (13), (12) и (14) в предположении, что все параметры равномерно распределены в указанных в таблице 3 интервалах. Модель (13) Модель (12) Модель (14) П-ры П-ры S ST OT S ST OT S ST OT ν 0.0026 0.0036 ν12 0.0030 0.0035 Ra 0.0631 0.0782 h 0.0677 0.0933 ν13 0.0001 0.0002 Rb 0.1431 0.1738 E 0.2151 0.2791 h 0.0222 0.0277 h 0.0295 0.0392 R1 0.3255 0.4093 E1 0.1607 0.2007 E1 0.1993 0.2416 ∆V 0.2663 0.3429 E3 0.0291 0.0441 E3 0.0357 0.0562 R1 0.3334 0.3943 ν12 0.0054 0.0089 ∆V 0.3589 0.4254 ν13 0.0001 0.0003 ∆V 0.4354 0.4944 П-ры 4. Калибровка параметра E1(E) Будем рассматривать модель (13), т. е. сферический изотропный случай. Наша задача минимизировать валидационную ошибку (Error), заданную формулой (15). Экспериментальные данные были получены именно для инъекциях объёма ∆V = 0.05 мл. Положим, что ν = 0.4, h распределён по закону LN (−7.4910, 0.0545) и R1 распределён по закону LN (−4.4714, 0.0513). Будем рассматривать ошибку (15) как функцию, зависящую от параметров функции распределения E. Считая, что E распределён последовательно равномерно, нормально и логнормально (см. Приложение, код MATLAB), было вычислено, что наименьшее значение ошибки соответствует логнормальному закону распределения параметра E с 26

Рис. 5: Результаты анализа чувствительности моделей (13), (12) и (14) в предположении, что все параметры равномерно распределены в указанных в таблице 3 интервалах. параметрами LN (15.8677, 0.4191). Результаты представлены в табл. 5. Зависимость ошибки от параметров нормального и логнормального распределения E можно посмотреть на рис. 6. Сравнение плотности вероятности отклика модели (13) изменения внутриглазного давления и плотности вероятности, полученной по данным эксперимента, при соответствующих законах распределения параметра E можно увидеть на рис. 7. Таблица 5: Валидационная ошибка при соответствующем законе распределения параметра E. Unif Normal Lognormal Error 851.7857 895.5186 177.7951 Рис. 6: Зависимость валидационной ошибки от параметров нормального и логнормального распределениий параметра E 5. Калибровка параметра E3 Действуя аналогично предыдущему пункту 4 и рассматривая модель (12), а так же полагая, что E1 распределён по закону распределения E, найденному в предыдущем пункте 4, будем теперь рассматривать валидационную ошибку как функцию, зависящую от параметров функции распределения E3 . Были рассмотрены последовательно равномерное, нормальное 27

Рис. 7: Сравнение плотности вероятности отклика модели (13) изменения внутриглазного давления и плотности вероятности, полученной по данным эксперимента, при соответствующих законах распределения параметра E и логнормальное распределение, и было получено, что наименьшее значение ошибки соответствует логнормальному закону распределения параметра E3 с параметрами LN (12.9613, 0.3156). Результаты представлены в табл. 6. Зависимость ошибки от параметров нормального и логнормального распределения E3 можно посмотреть на рис. 8. Сравнение плотности вероятности отклика модели (12) и плотности вероятности, построенной по данным эксперимента, при соответствующих законах распределения параметра E3 можно увидеть на рис. 9. Таблица 6: Валидационная ошибка при соответствующем законе распределения параметра E. Unif Normal Lognormal Error 843.6345 201.5694 139.3117 Рис. 8: Зависимость валидационной ошибки от параметров нормального и логнормального распределениий параметра E3 28

Рис. 9: Сравнение плотности вероятности отклика модели (12) и плотности вероятности, построенной по данным эксперимента, при соответствующих законах распределения параметра E3 6. Нахождение валидационной ошибки модели (14) На данном этапе нам известны функции распределения всех параметров модели (14), поэтому, полагая, что параметры соответствуют функциям распределения, указанным в таблице 7, вычислим валидационную ошибку. Сравнение плотностей вероятности модели (14) и эксперимента представлено на рис. 10. Валидационная ошибка равна = 128.9785. Таблица 7: Функции распределения параметров Пар-р ν12 (ν) ν13 h Ra (R1 ) Rb E1 (E) E3 ∆V Закон распр-ия Constant Constant Lognormal Lognormal Lognormal Lognormal Lognormal Constant Характеристики 0.4 0.01 -7.4910, 0.0545 -4.4714, 0.0513 -4.4714, 0.0513 15.8677, 0.4191 12.9613, 0.3156 0.05 мл 7. Анализ чувствительности моделей с откалиброванными параметрами Основываясь на ранее полученных результатах, будем полагать, что параметры соответствуют функциям распределения, указанным в таблице 29

Рис. 10: Сравнение плотностей вероятности модели (14) и эксперимента 7. Прежде всего исследуем плотности вероятности моделей (13), (12) и (14) (рис. 11). Рис. 11: Плотности вероятностей моделей (13), (12) и (14) в предположении, что все параметры распределены как указано в таблице 7 Плотности вероятности моделей 13 и 12 являются сильно асимметричными, поэтому для анализа чувствительности в этих случаях необходим метод PAWN; для анализа чувствительности модели (14) по-прежнему будет использоваться метод разложения в полиномиальный хаос (PCE). Результаты представлены на рис. 12 и в таблице 8. На рис. 13 можно посмотреть сравнение плотностей вероятности моделей (13), (12) и (14) с плотностью вероятности эксперимента. 30

Таблица 8: Результаты анализа чувствительности моделей (13), (12) и (14) PCE S Ra 0.0322 Rb 0.0567 h 0.0104 E1 0.8793 E3 0.0010 П-ры PAWN П-ры PAWN П-ры h 0.1600 E 0.9033 R1 0.3983 h E1 E3 R1 0.1617 0.8667 0.0883 0.4683 ST OT 0.0389 0.0677 0.0123 0.8994 0.0021 Рис. 12: Результаты анализа чувствительности моделей (13), (12) и (14) Рис. 13: Cравнение плотностей вероятности моделей (13), (12) и (14) с плотностью вероятности эксперимента 31

Заключение Результаты анализа чувствительности моделей (13), (12) и (14) представлены в табл. 8 и на рис. 12. На основе предоставленных экспериментальных данных (Приложение, табл. 9) вычислены плотности вероятности случайных величин, характеризующих геометрию оболочки глаза (рис. 3, табл. 2). С использованием введённого формулой (15) валидационного критерия вычислены параметры распределения характеристик материала моделей, модулей упругости (табл. 5, 6). Проведено сравнение плотностей вероятности изменения давления, найденных на основе рассматриваемых моделей, с данными эксперимента (Приложение, табл. 9), результаты представлены на рис. 13. В моделях (13), (12) и (14) глаз человека моделируется как шар или эллипсоид, заполненный несжимаемой жидкостью и покрытым оболочкой, состоящей из одного слоя. На основе результатов, представленных в таблицах 5 и 6, на рис. 7 и 9, можно заключить, что данный слой является трансверсально-изотропным. По результатам, представленным в таблицах 5, 6 и значении валидационной ошибки, полученной в пункте 6, можно сделать вывод, что изменение внутриглазного давления при интравитреальных инъекциях наилучшим образом описывается математической моделью в предположении, что форма глаза человека близка к эллипсоиду. По результатам анализа чувствительности, представленных в таблицах 4, 8 и на рис. 5, 12, можно заключить, что как на модели (13), (12), так и на модель (14) наибольшим влиянием обладают: объём инъекции ∆ V , внутренний радиус оболочки R1 и модуль упругости E(E1 ), а наименьшим — коэффициенты Пуассона ν12 (ν), ν13 . Следовательно, именно эти параметры должны быть определяющими в решении врачей о способе лечения и дозировке препарата. Законы распределения всех параметров моделей представлены в таблице 7. 32

Список литературы 1. Girard M. J. A., Dupps W. J., Baskaran M., Scarcelli G., Yun S. H., Quigley H. A., Sigal I. A., Strouthidis N. G. Translating ocular biomechanics into clinical practice: current state and future prospects // Current Eye Research, 2015. Vol. 40., No. 1, P. 1 – 18. 2. Kotliar K., Maier M., Bauer S., Feucht N., Lohmann C., Lanz I. Effect of intravitreal injections and volume changes on intraocular pressure: clinical results and biomechanical model // Acta Ophthalmologica Scandinavica, 2007. Vol. 85, P. 777 -– 781. 3. Kim J. E., Mantravadi A. V., Hur E.Y., Covert D. J. Short-term intraocular pressure changes immediately after intravitreal injections of anti–vascular endothelial growth factor agents // American Journal of Ophthalmology, 2008. Vol. 146, No. 6, P. 930 – 934. 4. Wu L., Evans T. Immediate changes in intraocular pressure after an intravitreal injection of 2.5 mg of bevacizumab // Archivos de la Sociedad Española de Oftalmologı́a (English Edition), 2010. Vol. 85, No. 11, P. 364 – 369. 5. Ozgur O. R., Ozkurt Y., Kulekci Z., Evciman T. The combination of phacoemulsification surgery and intravitreal triamcinolone injection in patients with cataract and diabetic macular edema // Saudi Journal of Ophthalmology, 2016. Vol. 30, P. 33 -– 38. 6. Avery R. L., Bakri S. J., Blumenkranz M. S., Brucker A. J., Cunningham E. T., D’Amico D. J., Dugel P. U., Flynn H. W., Freund K. B., Haller J. A., Jumper J. M., Liebmann J. M., McCannel C. A., Mieler W. F., Ta C. N., Williams G. A. Intravitreal injection technique and monitoring: Updated guidelines of an expert panel // Retina, The journal of retinal and vitreous diseases, 2014. Vol. 34, No. 12, P. S1 – S18. 7. Hamilton K., Pye D. C. Young’s modulus in normal corneas and the effect on applanation tonometry // Optometry and Vision Science, 2008. Vol. 85, No. 6, P. 445 – 450. 33

8. Bauer S. M., Voronkova E. B. Nonclassical shell theories in ocular biomechanics // Altenbach H., Mikhasev G.I., Shell and Membrane Theories in Mechanics and Biology. Cham: Springer International Publishing Switzerland, 2015. P. 81 – 97. 9. Иомдина Е. Н., Бауэр С. М., Котляр К. Е. Биомеханика глаза: теоретические аспекты и клинические приложения./ Под ред. В. В. Нероева. М.: Реал Тайм, 2015. 208 C. 10. Sigal I. A., Flanagan G. J., Ethier C. R Factors influencing optic nerve head biomechanics // Investigative Ophthalmology & Visual Science, 2005. Vol. 46, No. 11, P. 4189 – 4199. 11. Sigal I. A., Yang H., Roberts M. D., Burgoyne C. F., Downs J. C. IOPinduced lamina cribrosa displacement and scleral canal expansion: an analysis of factor interactions using parameterized eye-specific models // Investigative Ophthalmology & Visual Science, 2011. Vol. 52, No. 3, P. 1896 – 1907. 12. 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. Vol. 10, No. 3, P. 323 – 337. 13. Sigal I. A., Grimm J. L., Schuman J. S., Kagemann L., Ishikawa H., Wollstein G. A method to estimate biomechanics and mechanical properties of optic nerve head tissues from parameters measurable using optical coherence tomography // IEEE Transactions on medical imaging, 2014. Vol. 33, No. 6., P. 1381 – 1389. 14. Saltelli A., Annoni P. How to avoid a perfunctory sensitivity analysis // Environmental Modelling & Software, 2010. Vol. 25, No. 12, P. 1508 – 1517. 15. Morris M. D. Factorial sampling plans for preliminary computational experiments // Technometrics, 1991. Vol. 33, No 2., P. 161 – 174. 16. Campolongo F., Cariboni J., Saltelli A. An effective screening design for sensitivity analysis of large models // Environmental Modelling & Software, 2007. Vol. 22, P. 1509 – 1518 34

17. Khare Y. P, Munoz-Carpena R., Rooney R. W, Martinez C. J. A multicriteria trajectory-based parameter sampling strategy for the screening method of elementary effects // Environmental Modelling & Software, 2015. Vol. 64, P. 230 – 239. 18. Sobol’ I. M. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates // Mathematics and Computers in Simulation, 2001. Vol. 55, P. 271 – 280. 19. Saltelli A., Annoni P., Azzini I., Campolongo F., Ratto M., Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index // Computer Physics Communications, 2010. Vol. 181, P. 259 – 270. 20. Sudret B. Global sensitivity analysis using polynomial chaos expansions // Reliability Engineering & System Safety, 2008. Vol. 93, P. 964 -– 979. 21. Koekoek R., Lesky P. A., Swarttouw R. F. Hypergeometric orthogonal polynomials and their q-analogues. Berlin: Springer-Verlag Berlin Heidelberg, 2010. 578 с. 22. Blatman G., Sudret B. Adaptive sparse polynomial chaos expansion based on least angle regression // Journal of Computational Physics, 2011. Vol. 230, P. 2345 — 2367. 23. Borgonovo E., Castaings W., Tarantola S. Moment independent importance measures: new results and analytical test cases // Risk Analysis, 2011. Vol. 31, No. 3, P. 404 – 428 24. Pianosi F., Wagener T. A simple and efficient method for global sensitivity analysis based on cumulative distribution functions // Environmental Modelling & Software, 2015. Vol. 67, P. 1 – 11. 25. Родионова В. А., Титаев Б. Ф., Черных К. Ф. Прикладная теория анизотропных пластин и оболочек. СПб.: Изд-во С.-Петерб. ун-та, 1996. 280 C. 26. Ермаков А. М. Напряженно-деформированное состояние трансверсально-изотропных сопряженных эллиптических оболочек, находящих35

ся под действием внутреннего давления // Вестник С-Петерб. ун-та. Сер.1. Математика. Механика. Астрономия, 2009. Вып. 3. С. 110 – 118. 27. Anderson A. E., Ellis B. J., Weiss J. A. Verification, validation and sensitivity studies in computational biomechanics // Computer Methods in Biomechanics and Biomedical Engineering, 2007. Vol. 10, No. 3, P. 171 – 184. 28. Lee D., Ahn T. Statistical calibration of a finite element model for human middle ear // Journal of Mechanical Science and Technology, 2015. Vol. 29, No. 7, P. 2803 – 2815. 29. Fuest M., Kotliar K., Walter P., Plange N. Monitoring intraocular pressure changes after intravitreal Ranibizumab injection using rebound tonometry // Ophthalmic and Physiological Optics, 2014. Vol. 34, No. 4, P. 438 – 444. 30. Журавлева Д. И., Воронкова Е. Б. О методах анализа глобальной чувствительности математических моделей // Процессы управления и устойчивость: Труды 47-й международной научной конференции аспирантов и студентов. (В печати) 31. Efron B., Hastie T., Johnstone I., Tibshirani R. Least angle regression // The Annals of Statistics, 2004. Vol. 32, No. 2, P. 407 – 499. 36

Приложение Trajectory-based, алгоритм выбора точек X для метода элементарных эффектов 1: Создаётся матрица B = [m × k], состоящая из 0 и 1, такая, что для каждого столбца j = 1, k существуют две строки, отличающиеся друг от друга только в j-ом элементе, (траектория). ? 2: Генерируется диагональная матрица D = [k × k] такая, что каждый диагональный элемент равен +1 или -1 с равной вероятностью. 3: Задаётся единичная матрица Jmk = [m × k]. ? ? ? 4: Определяется X = (X1 , ..., Xk ) - базовое значение для X (начало траектории); каждый элемент принимает значение из {0, 1/(p − 1), 2/(p − 1), ..., 1 − ∆} с равной вероятностью. ? 5: Создаётся матрица P = [k × k] такая, что каждый столбец содержит один элемент равный 1, остальные - 0, причём одинаковых столбцов нет. ∆ ? 6: Вычисляется B = Jm1 X ? + [(2B − Jmk ) D? + Jmk ] P ? . Такая 2 матрица обладает теми же свойствами, что и матрица B, но её строки являются элементами ω. ? ? T 7: Повторяются этапы 1.–6. r раз, X = [B1 , ..., Br ] Матрица X задаёт r траекторий, в каждой из которых m точек, причём каждая последняя отличается от предыдущей лишь в одной компоненте. В итоге общее количество вычислений модели N = rm. В исходной статье [15] было предложено m = (k + 1). Пример расположения четырёх траекторий приведён на рис. 14 для функции Ishigami Homma (16). [30]. f (x) = sin(x1 ) + a sin2 (x2 ) + bx43 sin(x1 ), a = 7, b = 0.1, xi ∼ U nif orm[−π, π], i = 1, 3. 37 (16)

Sampling for Uniformity, алгоритм выбора точек X для метода элементарных эффектов 1: Зафиксируем p = 4 (количество уровней), ∆ = 2/3, r (количество траекторий) будем полагать чётным числом. Число точек в каждой траектории равно (k + 1). Также зададим Q - число повторений (в [17] предлагают 300 < Q < 500). 2: Устанавливаем начало и конец траекторий так, чтобы выполнялось: • на каждом уровне каждого параметра лежит начало и конец ровно r/4 или (r/2 − 1)/2 траекторий, • начало и конец траектории выбираются с учётом значения ∆. 3: 4: Произвольным образом, но с учётом значения ∆, генерируются промежуточные точки траекторий. За одно перемещение от точки к точке меняется только одна координата. Вычисляется Евклидово расстояние (ED) для набора траекторий. v u k k+1 k+1 X X uX 2 t Xia (z) − Xjb (z) . da,b = i=1 j=1 z=1 a, b - две траектории; Xia (z), Xjb (z) - z-ая координата точек Xia и Xjb , соответственно. v uX r u r X d2i,j . ED = t i=1 j=1 5: Q раз повторяются этапы 1.–3.. Оптимальным является набор траекторий с максимальным значением ED. Ψ? , A = LAR(p,q,X,Y,Ψ), выбор базиса в методе разложения в полиномиальный хаос k,p 1: aα = 0 ∀α ∈ A = Aq Ψ - базис-кандидат Ψ? = ∅ - оптимальный базис r0 = Y 2: Находим Ψαj , который наибольшим образом коррелирует с Y. ? 3: Добавляем Ψαj в Ψ , обновляем значения aα , α ∈ A: aα = aα + γw. Значения γ и w вычисляются, как показано в [31]. T ? 4: Обновляем r = a Ψ (X). Вычисляем Error. 5: Повторяем этапы 2.–3., пока размер оптимального базиса не будет равен min(p, N − 1), где N - размер выборки X. 38

LARS 1: Выбирается q, pmax . Генерируется X, вычисляется Y . ? 2: p = 0; Errormin = 1; Ψmin = ∅; pmin = 0 (опционально) Задаётся Errortgt . 3: while p 6= pmax & Errormin > Errortgt do 4: p=p+1 5: Применяется LAR(p,q,Y,Ψ = {Ψα |α ∈ Ak,p q }). ? 6: Результат предыдущего шага = Ψ . Используя данный результат, вычисляется Error по формуле (4). 7: if Error < Errormin then 8: Errormin = Error 9: Ψ?min = Ψ? 10: pmin = p 11: end if 12: if Error увеличивается уже два раза подряд (случай чрезмерной аппроксимации) then 13: p = pmin 14: Генерируется новый X, вычисляется Y . 15: end if 16: end while 17: По формуле (3) вычисляются коэффициенты a полиномиального разложения, соответствующие базизу Ψ?min . 39

Рис. 14: Пример расположения четырёх траекторий для функции Ishigami Homma. 1 %% F i n d i n g h p r o b a b i l i t y d i s t r i b u t i o n f u n c t i o n 2 % h - g i v e n data 3 b e t a = f i t d i s t ( h , 'Beta' ) ; gamma = f i t d i s t ( h , 'Gamma' ) ; 4 l o g n o r m a l = f i t d i s t ( h , 'Lognormal' ) ; normal = f i t d i s t ( h , 'Normal' ) ; 5 [ f , x i ]= k s d e n s i t y ( h ) ; 6 beta_f = pdf ( 'Beta' , xi , b e t a . a , b e t a . b ) ; 7 gamma_f = pdf ( 'Gamma' , xi , gamma.a , gamma.b ) ; 8 lognormal_f = pdf ( 'Lognormal' , xi , lognormal.mu , l o g n o r m a l . s i g m a ) ; 9 normal_f = pdf ( 'Normal' , xi , normal.mu , n o r m a l . s i g m a ) ; 10 beta_m = s q r t ( sum ( ( beta_f - f ) . ^2) ) ; gamma_m = s q r t ( sum ( ( gamma_f - f ) . ^2) ) ; 11 lognormal_m = s q r t ( sum ( ( lognormal_f - f ) . ^2) ) ; 12 normal_m = s q r t ( sum ( ( normal_f - f ) . ^2) ) ; 13 %% F i n d i n g R_1 p r o b a b i l i t y d i s t r i b u t i o n f u n c t i o n 14 % R_1 - g i v e n data 15 b e t a = f i t d i s t (R1 , 'Beta' ) ; gamma = f i t d i s t (R1 , 'Gamma' ) ; 16 17 l o g n o r m a l = f i t d i s t (R1 , 'Lognormal' ) ; normal = f i t d i s t (R1 , 'Normal' ) ; [ f , x i ]= k s d e n s i t y (R1) ; 18 beta_f = pdf ( 'Beta' , xi , b e t a . a , b e t a . b ) ; 19 gamma_f = pdf ( 'Gamma' , xi , gamma.a , gamma.b ) ; 20 lognormal_f = pdf ( 'Lognormal' , xi , lognormal.mu , l o g n o r m a l . s i g m a ) ; 21 normal_f = pdf ( 'Normal' , xi , normal.mu , n o r m a l . s i g m a ) ; 22 beta_m = s q r t ( sum ( ( beta_f - f ) . ^2) ) ; gamma_m = s q r t ( sum ( ( gamma_f - f ) . ^2) ) ; 23 lognormal_m = s q r t ( sum ( ( lognormal_f - f ) . ^2) ) ; 24 normal_m = s q r t ( sum ( ( normal_f - f ) . ^2) ) ; 40

1 %% SEARCH optim E_1 (E) v a l u e 2 3 x0 = [ ' f i g u r e g u e s s p o i n t ' ] ; options = optimset ('Display' ,' i t e r ') ; 4 [ x , f v a l , e x i t f l a g , output ] = f m i n s e a r c h (@SPHEREISO_005 , x0 , o p t i o n s ) ; 1 %% SPHEREISO_005 2 %% You w i l l need SAFE Toolbox t o run t h i s code 3 f u n c t i o n [ Err ] = SHEREISO_005( x ) 4 mu=x ( 1 ) ; sigma=x ( 2 ) ; 5 f u n _ t e s t = 'SPHEREISO' ; 6 GoldMU = 3543 . 9 0 9 2 6 8 4 2 1 1 ; Gold10 = 1546 . 5 3 5 2 ; Gold90 = 5414 . 2 0 6 4 2 ; 7 M=5; N=500; dV = 0 . 0 5 * ( 1 0 ^ ( - 6 ) ) ; nu = 0 . 4 ; 8 d i s t r _ f u n {1} = ' l o g n ' ; 9 d i s t r _ p a r {1} = [ -4 . 4 7 1 4 5 7 3 4 2 2 6 2 8 1 10 d i s t r _ f u n {2} = ' l o g n ' ; 11 d i s t r _ p a r {2} = [ -7 . 4 9 0 9 9 6 9 2 1 4 6 3 5 2 12 d i s t r _ f u n {3} = ' u n i f /norm/ l o g n ' ; 13 d i s t r _ p a r {3} = [mu sigma ] ; 0 .0513971925020194 ] ; 0 .0545691728478291 ] ; 14 X = AAT_sampling ( ' l h s ' ,M- 2 , d i s t r _ f u n , d i s t r _ p a r ,N) ; X = [X ... 15 Y = model_evaluation ( fun_test ,X) ; 16 Err = s q r t ( ( mean (Y) -GoldMU) ^2 + ( p r c t i l e (Y, 1 0 ) - Gold10 ) ^2 + ... repmat ( nu , N, 1 ) repmat (dV , N, 1 ) ] ; ( p r c t i l e (Y, 9 0 ) - Gold90 ) ) ; 1 %% SPHEREISO 2 f u n c t i o n P = SPHEREISO( x ) 3 R1 = x ( 1 ) ; h = x ( 2 ) ; E = x ( 3 ) ; nu = x ( 4 ) ; dV = x ( 5 ) ; 4 R2 = R1+h ; 5 dR = ( ( 3 / ( 4 * p i ) ) * ((4 * p i / 3 ) * (R1^3)+dV) ) ^ ( 1 / 3 ) - R1 ; 6 C1 = ( - ( R1^3) * (1 -2 * nu ) ) / (E * (R1^3 - R2^3) ) ; 7 C2 = ( ( R1^3) * (R2^3) *(1+nu ) ) / ( ( - 2 * E) * (R1^3 -R2^3) ) ; 8 u = C1*R1 + C2/ (R1^2) ; P = dR/u ; 9 41

Таблица 9: Экспериментальные данные 39.9966 3653.0228 4359.6294 6386.1238 3986.3278 4826.2564 3946.3312 2213.1452 11.238 10.786 11.86 11.2615 11.191 11.644 10.586 11.701 11.331 11.518 11.5465 11.4085 11.89 0.60415 0.60815 0.54415 0.59715 0.63015 0.52915 0.54115 0.61115 0.56715 0.57115 0.57515 0.58915 0.52015 5506.1986 3226.3924 3386.3788 1039.9116 2599.779 3972.9956 3666.355 3453.0398 ∆P, Па 6106.1476 2973.0806 3679.6872 4639.6056 3853.0058 3039.7416 5199.558 1306.5556 5906.1646 4172.9786 3319.7178 3119.7348 5106.2326 1893.1724 1733.186 2399.796 1466.542 4012.9922 2639.7756 3533.033 5186.2258 3119.7348 11.3685 11.659 11.78 11.4865 11.203 11.2655 10.678 10.953 11.1235 11.1535 11.8575 13.7295 11.9485 R1 , мм 10.9235 11.6345 11.928 11.5215 17.17 11.2005 11.8045 10.9095 11.0565 11.0845 11.2425 13.7305 10.35 10.9495 12.7355 11.849 11.3025 11.3175 11.5555 11.883 11.223 10.8955 11.2495 11.2685 16.1795 10.3705 10.754 12.4425 11.318 10.9975 11.764 11.5105 11.6915 11.25 11.4955 11.235 11.586 15.2465 h, мм 0.57315 0.56115 0.54415 0.57515 0.58215 0.52915 0.55715 0.58015 0.53915 0.53015 0.54815 0.52715 0.57115 0.56915 0.55715 0.53815 0.58215 0.53815 0.52815 0.49815 0.54415 0.50715 0.57315 0.52015 0.59215 0.58515 0.59715 0.54415 0.53915 0.54415 0.56415 0.51315 0.53315 0.50515 0.59115 0.52315 0.59315 0.55215 0.53215 0.57515 0.64515 0.54915 0.54415 0.57415 0.55915 0.54115 0.58315 0.57115 0.53915 42

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв