Сравнение градиентного и симплекс методов численного решения обратной задачи для простейшей модели инфекционного заболевания Текст научной статьи по специальности «Математика»

Математические заметки СВФУ Апрель—июнь, 2015. Том 22, № 2

УДК 519.622

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

C. И. Кабанихин, О. И. Криворотько, Д. В. Ермоленко, Д. А. Воронов

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

В работе численно исследована задача определения параметров инфекционного заболевания для простейшей математической модели «антиген-антитело» по измерениям концентраций антигенов и антител в фиксированные моменты времени. Исследован целевой функционал, описывающий отклонение экспериментальных данных от модельных. Получено явное выражение градиента целевого функционала, в котором используется решение соответствующей сопряженной задачи. Проведен сравнительный анализ численного решения обратной задачи, полученного градиентным методом (итерации Ландвебера) и симплекс-методом (методом Нелдера — Мида). Показано, что метод Нелдера — Мида в исследуемой математической модели определяет множество локальных приближенных значений скоростей распространения антигена, иммунного ответа и выработки специфичных антител с заданной точностью. Метод итерации Ландвебера находит ближайший к начальному приближению аргумент минимума целевого функционала за достаточное большое число итераций.

Ключевые слова: обратная задача, оптимизационный подход, метод итерации Ландвебера, метод Нелдера — Мида, моделирование в иммунологии.

S. I. Kabanikhin, O. I. Krivorotko, D. V. Yermolenko, D. A. Voronov. Comparison of gradient and simplex methods of the numerical solution of the inverse problem for the simplest model of infectious disease. Abstract: The organism releases antibodies after antigens are recognized. This antibodies help to cope with disease. Individual characteristics of immunity and disease are different. Therefore the reaction of different people to the disease is not similar. In spite of this doctors often make a standard treatment plan for patients. Therefore it is important to be able to identify individual parameters of immunity and disease for each patient by blood, urine tests and so on. In this paper the problem for determining

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 16—31—00382).

© 2015 Кабанихин C. И., Криворотько О. И., Ермоленко Д. В., Воронов Д. А.

of the immune and disease parameters by measurements of the antigen and antibody concentrations in fixed times is numerically investigated. We use the misfit functional that describes the deviation between experimental and model data. The explicit expression of the gradient of the misfit functional based on the adjoint problem solution is obtained. The results of numerical calculations of inverse problem is discussed. The results obtained using the methods of Landveber's iterations and Nelder-Mead. Comparative analysis of these methods is discussed. It was demonstrated that Nelder-Mead method is more sensitive to the choice of the initial guess than the method of Landweber iteration. If the initial guess has been chosen "good enough" than it needs less steps than Landweber iteration method. But there are some set of possible initial guesses when Landweber iteration method converges but the Nelder-Mead method gives no results.

Keywords: inverse problem, optimization approach, Landveber iterations, Nelder-Mead method, modeling in immunology.

Введение

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

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

Математические модели иммунологии, включая численное решение прямых и обратных задач, исследовали Г. И. Марчук [1], А. А. Романюха [2], S. M. Andrew [3], H. W. Engl [4], C. Molina-Paris, G. Lythe [5] и др. H. T. Banks, S. Hu [6] использовали прямые методы численного решения задачи наименьших квадратов со случайным распределением данных. Г. П. Кузнецова [7] в своей работе использовала метод численного интегрирования обратной задачи для простейшей модели инфекционного заболевания Г. И. Марчука. В [8] численно исследована обратная задача для простейшей математической модели взаимодействия антител и антигенов градиентным методом. Обоснованы оценки сходимости алгоритма, доказана теорема единственности и локальной устойчивости. Основная цель работы — анализ двух алгоритмов идентификации параметров простейшей математической модели для получения информации о характере заболевания и об иммунном ответе по данным анализа крови. Зада-

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

В работе исследовано численное решение обратной задачи для простейшей модели инфекционного заболевания (так называемой модели «антиген-антитело»), состоящей из двух нелинейных дифференциальных уравнений. Данная модель позволяет качественно описать взаимодействие антигенов и антител в организме. В работе численное решение обратной задачи получено с помощью методов итерации Ландвебера и Нелдера — Мида. Статья организована следующим образом. В разд. 1 представлена постановка обратной задачи для простейшей модели инфекционного заболевания. В разд. 2, 3 исследованы два метода решения обратной задачи. В частности, в разд. 2 описан градиентный метод (метод итерации Ландвебера) и представлены численные результаты, полученные с помощью данного метода, а в разд. 3 изучен метод Нелдера — Мида и приведено численное решение обратной задачи. В разд. 4 проведен сравнительный анализ методов итерации Ландвебера и Нелдера — Мида.

1. Постановка обратной задачи

Исследована задача Коши для простейшей модели инфекционного заболевания «антиген-антитело» [8,9]:

= КШРи - Р¿€(0,Г), ^ =/321^1 т2Ц), *€(0,Г), (1)

N1(0) = N10, N2(0) = N20, которая может быть записана в векторном виде:

^ =Р(Ж(4),/3), ¿€(0 ,Г), N (0) = N0.

Здесь N(4) = (N1(4), ^(¿))т — переменные системы (концентрации антигенов и антител в организме), в = (вш в12, в21)Т — вектор параметров, характеризующий особенности иммунитета: ви отражает рост числа антигенов в организме, в12 — скорость иммунного ответа, в21 — скорость выработки специфичных антител, Р — заданная вектор-функция. Задача моделирования (2) при заданных в и N° называется прямой задачей.

Пусть в фиксированные моменты времени ^, к = 1,..., К, измерены концентрации антигенов N1(4) и антител N2(4) (обозначим N¿(4) = N¿(4; в), г = 1, 2):

; в)= Фг(^), г = 1, 2; к =1,...,К. (3)

Обратная задача (2), (3) заключается в определении вектора параметров в по заданной функции Р, начальным данным N° и дополнительной информации вида (3). Определим оператор обратной задачи (2), (3):

А : & ^ Кк,

где & := {в € К3 : в^ > 0, г,^ = 1,2} — пространство рассматриваемых параметров.

(2)

Запишем обратную задачу (2), (3) в операторном виде:

А(в) = Ф, Ф =(Ф1(*1),...,Ф1(*К ),Ф2(*1),..., ))T. (4)

Вектор Ф определяется, например, по данным анализа крови и мочи в моменты времени tk, k = 1,...,K. Решение обратной задачи (4) будем искать, минимизируя целевой функционал J(в) = ||А(в) — Ф||2, определенный следующим образом:

к

J (в) = Е |N (tk; в) — ФЫ|2. (5)

k=0

Это означает, что решение прямой задачи (2) при оптимальном в в моменты времени tk, k = 1,..., K, будет максимально близко к измерениям состояний системы (концентраций антигенов N1(t) и антител N2(t)) в те же моменты времени tk.

2. Численное решение обратной задачи методом итерации Ландвебера

Используем метод итерации Ландвебера для решения задачи min J(в), ко-

ße-91

торый заключается в определении приближенного решения [10,11]:

вп+1 = вп — aJ'(вп), а > 0, во е ^, (6)

где а — параметр спуска, J'(в) е R3 — градиент целевого функционала (5), явное выражение которого определяется следующим образом [12]:

т

J'^) = —J Ф(t)TPß(N(t), в) dt. (7)

0

Здесь Ф (t) — решение сопряженной задачи

^ = -P£(JV(t), /J)*(t), t G U tfc+i), io = 0, ttf+1 = T,

k=0 Ф (T) = 0,

[Ф]t=tfc = 2(N(tk; в) — ФЫ), k = 1,..., K,

(8)

где PN(N(t), в) е R2 x R2 и Pe(N(t), в) е R2 x R3 — соответствующие матрицы Якоби:

в11 — вl2N2(t) —вl2Nl(t) А = ( N1(t) —N1(t)N2(t) 0 в21^С0 в21^(0 J ' в I 0 0 N1(t)N2(t)

[Ф]4=4к := Ф+ е) — Ф— е) — разрыв функции Ф в точке где 7 > 0 — сколь угодно малое число.

Для численного решения прямой (2) и сопряженной (8) задач используем метод Рунге — Кутты четвертого порядка аппроксимации. Построим равномерную сетку ш := ^^ = ¡ц = Т/А^, ] €= 0, А^}. Пусть время моделирования Т составляет 4 недели, N = 100 — количество точек сетки ш, а = 0.001, е3 = 10-6 — параметр остановки итерационного процесса, N0 = (1.8,1.8) —

начальные данные. Выберем вектор параметров в = (0.5, 0.5, 0.6)Т, описывающий иммунитет среднестатистического человека, который будем в дальнейшем называть точным решением обратной задачи (2), (3). Используем равномерно распределенные по сетке ш синтетические данные N , в) и N2 , в) в моменты времени к = 1,..., К, в качестве вектора данных Ф.

Алгоритм численного решения обратной задачи (2), (3) методом итерации Ландвебера состоит в следующем.

1. Задаем начальное приближение во = (0.1, 0.1, 0.2)Т, описывающее легкую форму инфекции. Решаем прямую задачу (2) при заданном во. Строим вектор N во), к = 1,...,К.

2. С помощью принципа математической индукции показываем, как вычислить приближение в«+1, зная в«.

3. Решаем прямую задачу (2) для набора параметров в«, т. е. находим N в), к = 1,...,К.

4. Если J(в«) < £в, то в« можно считать приближенным решением.

5. Если J(в«) > £в, то решаем сопряженную задачу (8), полагая в = в«.

6. Определяем J'(в«) по формуле (7).

7. Определяем в«+1 согласно соотношению (6).

Исследуем отклонение N¿(4; в«) — N¿(4; в)| расчетных кривых N¿(4; в«) от «экспериментальных» N¿(4; в) в зависимости от числа измерений К, г = 1, 2. Рис. 1 показывает, что при увеличении количества измерений отклонение расчетных кривых N¿(4; в«) от «экспериментальных» N¿(4; в) уменьшается. Однако сделать 35 измерений за 4 недели (т. е. 35 раз взять анализы) проблематично. Согласно [2] выберем максимальную меру такого отклонения, равную 7 • 10-4, т. е. в«) — N¿(4;в)| < 7 • 10-4. Таким образом, в дальнейшем для всех

расчетов будем использовать 20 измерений за 4 недели.

t t

Рис. 1. Графики |N¿(4; — N¿(4; в)| в зависимости от числа измерений К: слева г =1, справа г = 2.

Исследуем теперь относительную погрешность вычисления решения обратной задачи |в« — в|/|в| , которая является безразмерной величиной и равна

отношению абсолютной погрешности решения обратной задачи к точному решению обратной задачи. Рис. 2 показывает, что при увеличении числа измерений значение относительной погрешности уменьшается, т. е. приближенное решение обратной задачи становится ближе к точному. Отметим связь между выбранными мерами качества решения обратной задачи, а именно между отклонением расчетных кривых от «экспериментальных» |N¿(4; вп) — N¿(4; в)| (см. рис. 1), и относительной погрешностью |вп — в|/|в| (см. рис. 2). В дальнейшем в качестве критерия выбора оптимального количества измерений для численного решения обратной задачи (2), (3) будем использовать зависимость относительной погрешности от количества измерений К.

0.007 г-I-I-I-I-I-

0.0065 -\ -

0.006 - -

0.0055 - \ -

0.005 - ё\ -

0.0045 - \ч -

0.004 - -

0.0035 - -

0.003 - _____ -

0.0025 -1-1-1-1-1-^

5 10 15 20 25 30 35

К

Рис. 2. График зависимости относительной погрешности |вп — в|/|в| от числа измерений К.

На рис. 3 изображены графики зависимости целевого функционала J(вп) от числа итераций п и от абсолютной погрешности |вп — в|. На рис. 3 слева видно, что на первых двух итерациях функционал резко возрастает (ввиду слабой устойчивости обратной задачи (2), (3)), а начиная с третьей итерации монотонно убывает со скоростью 1/п, что свидетельствует о сходимости метода. Отметим, что при этом абсолютные погрешности |вп — в| монотонно убывают. На рис. 3 справа показано, что чем больше абсолютная погрешность |вп — в| (на первых итерациях), тем больше отклонение модельных данных от «экспериментальных» J (в„).

На рис. 4 приведены результаты численного решения обратной задачи (2), (3) при К = 20. Можно заметить, что в этом случае получаем приближенное решение: вп = 0.49675, вГ2 = 0.49903 и в^ = 0.60017 за 26836 итераций.

Отметим, что для любого начального приближения метод итерации Ланд-вебера сходится к нормальному решению обратной задачи (2), (3) [13]. В случаях а = 10-4 и а = 10-5 численное решение обратной задачи (2), (3) незначительно отличается от приведенных результатов, а время выполнения алгоритма на ЭВМ значительно больше, чем для случая а = 10-3.

3. Численное решение обратной задачи методом Нелдера — Мида

Метод Нелдера — Мида (симплекс-метод) [14] — метод безусловной оптимизации функционала, не использующий градиента. Ввиду этого метод Нелдера —

Рис. 3. График 3(вп) при общем числе итераций п = 26836, К = 20 (слева). График зависимости 3(вп) от абсолютной погрешности |вп — в| при общем числе итераций п = 26836, К = 20 (справа).

0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1

10000 20000 п

Рис. 4. Графики изменения параметров вп (слева), вп2 (в центре), вп1 (справа) в зависимости от числа итераций, п = 26836.

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

Решение обратной задачи, как и в случае градиентного метода, будем искать, минимизируя целевой функционал (5). Тем самым требуется найти безусловный минимум функции трех переменных J(вп, в12, в21).

В данном разделе представлены результаты численных расчетов, полученные с помощью метода Нелдера — Мида при таких же параметрах моделирования, как и в случае метода итерации Ландвебера (выбор сетки ш, начальных условий прямой задачи (2), точного вектора параметров в и параметра остановки ев).

В случае выбора начального симплекса в1 = (0.1, 0.1,0.2)т, в2 = (0.4,0.7,0.8)т, вз = (0.2,0.3,0.4)т, в4 = (0.9,0.7, 0.1)т

0

10000

20000

0

0

10000

20000

п

п

численное решение обратной задачи (2), (3) методом Нелдера — Мида при 10 измерениях (случай наименьшей относительной погрешности) следующее: вп сходится к 0.17, вп к 0.4, вП к 0.62. Можно заметить, что разность полученного решения обратной задачи и точного достаточна велика. Следовательно, метод Нелдера — Мида для данного начального симплекса определил локальный минимум.

Теперь проведем подобные вычисления для начального симплекса: в 1 = (0.15,0.2, 0.35)т, в2 = (0.05, 0.1, 0.3)т, вз = (0.05,0.1, 0.3)т, в4 = (0.3,0.3, 0.2)т.

На рис. 5 слева изображен график зависимости относительной погрешности |вп — в|/|в| от числа измерений К. По данному графику можно проследить, что при 25 измерениях относительная ошибка наименьшая, но сделать 25 измерений за 4 недели проблематично. Поэтому возьмем К = 20, как в разд. 2, и для дальнейших вычислений будем использовать это число измерений. Метод Нелдера — Мида сходится для данного набора параметров, т. е. функционал ■ (вп) стремится к нулю (см. рис. 5 справа).

20 К

25 30 35

0 10 20 30 40 50 60 70

5 10 15

п

Рис. 5. График зависимости относительной погрешности |вп — в|/|в| от числа измерений К (слева).

График J(вп) при общем числе итераций п = 72 (справа).

На рис. 6 представлено решение обратной задачи в зависимости от числа итераций п. Отметим, что полученные результаты близки к точному решению в = (0.5, 0.5, 0.6) . Следовательно, используя начальный симплекс в1 = (0.15, 0.2, 0.35)т, в2 = (0.05, 0.1, 0.3)т, вз = (0.05, 0.1, 0.3)т, в4 = (0.3,0.3,0.2)т, мы нашли минимум функционала ■ .

Представленные примеры показывают, что результаты, полученные с помощью метода Нелдера — Мида, зависят от выбора начального приближения. В зависимости от начального симплекса можно попасть как в локальный минимум, так и в глобальный. Это основная проблема данного метода. При работе с реальными данными метод Нелдера — Мида не гарантирует нахождения глобального минимума целевого функционала. Данную проблему помогают решить стохастические методы, такие как метод Монте-Карло [15], генетический алгоритм [16] и др.

4. Сравнительный анализ метода Нелдера — Мида с методом итерации Ландвебера

В табл. 1 представлен анализ численного решения обратной задачи (2),

0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2

10 20 30 40 50 60 70 п

10 20 30 40 50 60 70 п

10 20 30 40 50 60 70 п

Рис. 6. Графики изменения параметров в™1 (слева), в™2 (в центре), в21 (справа) в зависимости от числа итераций п = 72.

0

0

0

(3), полученный методами Нелдера — Мида и итерации Ландвебера для 20 измерений. Можно заметить, что метод Нелдера — Мида получает ответ гораздо быстрее, чем градиентный метод. Относительная ошибка для метода Нелдера — Мида на порядок ниже, чем для метода итерации Ландвебера. Метод итерации Ландвебера сходится к нормальному относительно начального приближения решению обратной задачи (2), (3) для любого начального приближения [13]. Однако вычислительное время метода итерации Ландвебера в несколько раз превышает время, требуемое для метода Нелдера — Мида ввиду сложности алгоритма (вычисление градиента целевого функционала).

Таблица 1. Анализ методов Нелдера — Мида и итерации Ландвебера решения обратной задачи (2), (3)

Метод Нелдера — Мида Метод итерации Ландвебера

К, число измерений 20 20

Начальное приближение /3 С15 = (0.15,0.2,0.35)т, /3(2) = (0.05, 0.1, 0.3)т, /3(3) = (0.2, 0.07,0.25)т, /3(4' = (0.3, 0.3,0.2)т /Зо = (0.1, 0.1, 0.2)т

е3, параметр остановки 10~6 10~6

/ЗГх 0.50059584 0.49675121

¿»5*2 0.50013173 0.49902571

/3?! 0.59992131 0.60017153

\Рп ~ Р\/\Р\, относительная погрешность 0.00066348 0.00366208

п, общее число итераций 72 26836

Время выполнения алгоритма (сек) 0.013 4.882

5. Заключение

В работе исследовано численное решение обратной задачи для простейшей математической модели «антиген-антитело», полученное с помощью двух итерационных методов: итерации Ландвебера и Нелдера — Мида. В рассмотренных численных экспериментах показано, что метод Нелдера — Мида определяет множество локальных приближенных значений скоростей распространения антигена, иммунного ответа и выработки специфичных антител с заданной точностью. Метод итерации Ландвебера находит ближайший к начальному приближению аргумент минимума целевого функционала. Таким образом, построен численный алгоритм, позволяющий уточнять параметры простейшей математической модели (скорости распространения антигенов, иммунного ответа и выработки специфичных антител) по 20 измерениям концентраций антигенов и антител в течение 4 недель за 5 с на персональном компьютере с процессором Intel (R) Core (TM) i3 2.13GHz и оперативной памятью 4 Гб.

ЛИТЕРАТУРА

1. Марчук Г. И. Математические модели в иммунологии. М.: Наука, 1980.

2. Романюха А. А., Руднев С. Г., Зуев С. М. Анализ данных и моделирование инфекционных заболеваний // Современные проблемы вычислительной математики и математического моделирования. Т. 2. Математическое моделирование (под ред. В. П. Дымникова). М.: Наука, 2005. С. 352-404.

3. Andrew S. M., Baker C. T. H., Bocharov G. A. Rival approaches to mathematical modeling in immunology // J. Comput. Appl. Math. 2007. V. 205. P. 669-686.

4. Engl H. W., Flamm C., Kugler P., Lu J., Muller S., Schuster P. Inverse problems in systems biology // Inverse Probl. 2009. V. 25. P. 123014.

5. Molina-Paris C., Lythe G. Mathematical models and immune cell biology. New York; Dordrecht; Heidelberg; London: Springer-Verl., 2011.

6. Banks H. T., Hu Sh., Clayton Thompson W. Modeling and inverse problems in the presence of uncertainty. Boca Raton; London; New York; Washington, DC: CRC Press, 2014. (Monogr. Res. Notes Math.).

7. Kuznetsova G. P. The inverse problem for the Marchuk immunologic "simplest model" // Дальневост. мат. журн. 2003. Т. 4, № 1. С. 134-140.

8. Afraites L., Atlas A. Parameters identification in the mathematical model of immune competition cells // J. Inverse Ill-posed Probl. 2015. V. 23, N 4. P. 323-337.

9. Bellouquid A., Delitala M. Modeling complex biological systems: A kinetic theory approach. Boston; Basel; Berlin: Birkhauser, 2006.

10. Алифанов О. М., Артюхин Е. А., Румянцев С. В. Экстремальные методы решения некорректных задач. М.: Наука, 1988.

11. Кабанихин С. И. Обратные и некорректные задачи. Новосибирск: Сиб. науч. изд-во, 2009.

12. Ильин А. И., Кабанихин С. И., Криворотько О. И. Об определении параметров моделей, описываемых системами нелинейных дифференциальных уравнений // Сиб. электрон. мат. изв. 2014. Т. 11. С. 1-14.

13. Васин В. В. О сходимости методов градиентного типа для нелинейных уравнений // Докл. РАН. 1998. Т. 359, № 1. С. 7-9.

14. Nelder J., Mead R. A simplex method for function minimization // Comput. J. 1965. V. 7. P. 308-313.

15. Михайлов Г. А., Войтишек А. В. Численное статистическое моделирование. Методы Монте-Карло. М.: Академия, 2006.

16. Гладков Л. А., Курейчик В. В., Курейчик В. М. Генетические алгоритмы. М.: Физматлит,

2006.

Статья поступила 28 августа 2015 г.

Кабанихин Cергей Игоревич, Криворотько Ольга Игоревна, Ермоленко Дарья Владимировна, Воронов Дмитрий Андреевич Институт вычислительной математики и математической геофизики СО РАН,

проспект Академика Лаврентьева, 6, Новосибирск 630090;

Новосибирский гос. университет,

ул. Пирогова, 2, Новосибирск 630090

ksi52@mail.ru, krivorotko.olya@mail.ru,

ermolenko.dasha@mail.ru, dmitriy.voronov@gmail.com