1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 60, № 4, 2020

Журнал вычислительной математики и математической физики, 2020, T. 60, № 4, стр. 590-600

Оптимизационные методы решения обратных задач иммунологии и эпидемиологии

С. И. Кабанихин 1, 2, О. И. Криворотько 1, 2, *

1 Институт вычислительной математики и математической геофизики СО РАН
630090 Новосибирск, пр-т Акад. Лаврентьева, 6, Россия

2 Новосибирский государственный университет
630090 Новосибирск, ул. Пирогова, 2, Россия

* E-mail: krivorotko.olya@mail.ru

Поступила в редакцию 10.10.2019
После доработки 09.12.2019
Принята к публикации 16.12.2019

Полный текст (PDF)

Аннотация

В работе исследованы обратные задачи для систем нелинейных обыкновенных дифференциальных уравнений, которые заключаются в определении неизвестных коэффициентов и начальных данных по дополнительной информации о решении соответствующих прямых задач, измеренной в заданные моменты времени. Приведены примеры обратных задач иммунологии и эпидемиологии, возникающих при исследовании развития инфекционных заболеваний, динамики ВИЧ и распространения туберкулеза в особо эндемичных районах с учетом лечения. В случае, когда решение обратной задачи неединственно, рассматриваются три подхода к исследованию идентифицируемости математических моделей. Предложен алгоритм численного решения, основанный на минимизации квадратичного целевого функционала. На первом этапе находятся окрестности точек глобального минимума, на втором применяются градиентные методы. Градиент целевого функционала вычисляется через решение соответствующей сопряженной задачи. Приведены результаты численных расчетов. Библ. 31. Фиг. 1.

Ключевые слова: обратные задачи, ОДУ, идентификация параметров, градиентный метод, градиент функционала, иммунология, эпидемиология.

1. ВВЕДЕНИЕ

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений (ОДУ):

(1)
$\frac{{dy(t)}}{{dt}} = \Phi (y(t),t,\varphi ),\quad t \in ({{t}_{0}},T),$
(2)
$y({{t}_{0}}) = \psi .$

Предположим, что некоторые компоненты векторов $\varphi \in {{\mathbb{R}}^{m}}$ и $\psi \in {{\mathbb{R}}^{n}}$ неизвестны. Эти неизвестные будем обозначать через $q$. В обратной задаче требуется определить $q$ в случае, когда решение прямой задачи (1), (2) задано в некоторых точках интервала $({{t}_{0}},T)$:

(3)
${{y}_{i}}({{t}_{k}}) = {{f}_{{ik}}},\quad i \in I \subseteq {{I}_{n}},\quad k \in {{I}_{K}}.$
Здесь и далее ${{I}_{n}} = {\text{\{ }}1,2,\; \ldots ,\;n{\text{\} }}$.

Обратные задачи (1)–(3) возникают в биологии и медицине (включая иммунологию [1]–[8], эпидемиологию [9], [10], фармакокинетику [11]). Вектор $y(t;q) = y(t) = ({{y}_{1}}(t),\; \ldots ,\;{{y}_{N}}(t))$ характеризует концентрацию антигенов и иммунных составляющих (иммунология), количество заболевших в процессе распространения эпидемии (эпидемиология). Коэффициенты модели $\varphi $ описывают скорость воспроизводства антител, антигенов, регенерации клеток (иммунология), развития заболевания в популяции (эпидемиология). Решение обратной задачи позволяет получить информацию о характере заболевания, об иммунном ответе, о скорости развития эпидемии.

В работе [12] для построения решения обратных задач для динамических систем используются идеи метода позиционного управления из теории дифференциальных игр [13] и показывается, что этот подход применим для устойчивого восстановления $q$, реализующих невозмущенное измерение движения. В работе [14] использованы градиентные методы и получены критерии разрешимости задачи о гарантированном получении целевого множества решений обратной задачи. В данной работе предлагаются алгоритмы решения коэффициентных обратных задач для систем ОДУ, описывающих процессы в биологии.

В разд. 2 приведены примеры математических моделей инфекционного заболевания [1], динамики ВИЧ [5] и эпидемии туберкулеза [9]. В разд. 3 приведены теоремы о корректности прямой задачи. В разд. 4 обсуждены вопросы идентифицируемости исследуемых обратных задач. В разд. 5 обратная задача (1)–(3) сведена к вариационной постановке, получено явное выражение градиента целевого функционала, связанное с решением соответствующей сопряженной задачи. Результаты численных экспериментов обсуждаются в разд. 6.

2. ПРИМЕРЫ ЗАДАЧ

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

2.1. Математическая модель инфекционного заболевания Г.И. Марчука

Задача Коши для математической модели Г.И. Марчука [1] простейшего инфекционного заболевания (течение острой пневмонии), характеризующая изменение концентраций патогенных размножающихся антигенов $V(t)$, антител $F(t)$, плазматических клеток $C(t)$ и относительной характеристики пораженного органа $m(t)$, имеет вид

(4)
$\begin{gathered} \dot {V} = (\beta - \gamma F(t))V(t), \\ \dot {F} = \rho C(t) - ({{\mu }_{f}} + \eta \gamma V(t))F(t), \\ \dot {C} = \xi (m)\alpha \left( {F(t)V(t) - \tau F(t)\dot {V}(t) - \tau V(t)\dot {F}(t)} \right) - {{\mu }_{c}}(C(t) - C{\kern 1pt} {\text{*}}), \\ \dot {m} = \sigma V(t) - {{\mu }_{m}}m(t), \\ V(0) = {{V}^{0}},\quad C(0) = C{\kern 1pt} *,\quad F(0) = \rho C{\kern 1pt} {\text{*/}}{{\mu }_{f}},\quad m(0) = 0,\quad t \in (0,T). \\ \end{gathered} $
Здесь $\tau = 1$ сутки – время, в течение которого осуществляется формирование каскада плазматических клеток (параметр запаздывания), $\alpha = 5 \times {{10}^{{ - 11}}}$ [мл/(моль · сут)] – коэффициент, учитывающий вероятность встречи антиген-антитело, $a = 0.35$ [1/сут] – коэффициент размножения антигенов, $\gamma = 8.5 \times {{10}^{{ - 14}}}$ [мл/(моль · сут)] – коэффициент, связанный с вероятностью нейтрализации антигена антителами при встрече с ним, ${{\mu }_{c}} = 0.5$ [1/сут] и ${{\mu }_{f}} = 0.05$ [1/сут] – коэффициенты, обратно пропорциональные времени жизни плазматических клеток и антител, соответственно $\rho = 7000$ [1/сут] – скорость производства антител одной плазматической клеткой, $\eta = 20$ – количество антител, необходимых для нейтрализации одного антигена, $\sigma = 9 \times {{10}^{{ - 9}}}$ [мл/(моль · сут)] – скорость поражения органов-мишеней антигеном, ${{\mu }_{m}} = 0.4$ [1/сут] – скорость регенерации органа-мишени и $C{\kern 1pt} * = 2850$ [кл/мл] – концентрация плазмоклеток в здоровом организме 18-летнего инфицируемого. Задана функция, отвечающая за производительность выработки антител $\xi (m) = 1$ при $m \in [0,0.3]$ и $\xi (m) = 1 - (m - 0.3){\text{/}}0.7$ при $m \in (0.3,1]$. Начальная доза поражения ${{V}^{0}} = 1000$ [моль/мл].

В случае $y = (V,F,C,m)$ ($n = 4$), $\varphi = \left( {\tau ,\alpha ,\beta ,\gamma ,\rho ,\eta ,{{\mu }_{c}},{{\mu }_{f}}} \right)$ ($m = 8$) математическая модель (4) сводится к виду (1). Остальные параметры модели $C{\kern 1pt} *$, $\sigma $, ${{\mu }_{m}}$ в дальнейшем считаются известными.

2.2. Модель динамики ВИЧ

Математическая модель динамики ВИЧ-инфекции [5] описывается системой нелинейных дифференциальных уравнений:

$\begin{gathered} {{{\dot {T}}}_{1}} = {{\lambda }_{1}} - {{d}_{1}}{{T}_{1}} - {{k}_{1}}V{{T}_{1}}, \\ {{{\dot {T}}}_{2}} = {{\lambda }_{2}} - {{d}_{2}}{{T}_{2}} - {{k}_{2}}V{{T}_{2}}, \\ \end{gathered} $
(5)
$\begin{gathered} \dot {T}_{1}^{ * } = {{k}_{1}}V{{T}_{1}} - \delta T_{1}^{ * } - {{m}_{1}}ET_{1}^{ * }, \\ \dot {T}_{2}^{ * } = {{k}_{2}}V{{T}_{2}} - \delta T_{2}^{ * } - {{m}_{2}}ET_{2}^{ * }, \\ \end{gathered} $
$\begin{gathered} \dot {V} = {{N}_{T}}\delta (T_{1}^{ * } + T_{2}^{ * }) - cV - [{{\rho }_{1}}{{k}_{1}}{{T}_{1}} + {{\rho }_{2}}{{k}_{2}}{{T}_{2}}]V, \\ \dot {E} = {{\lambda }_{E}} + {{b}_{E}}{{\alpha }_{b}}E + {{d}_{E}}{{\alpha }_{d}}E - {{\delta }_{E}}E,\quad t \in (0,T). \\ \end{gathered} $
Здесь ${{T}_{1}}$, $T_{1}^{ * }$ – Т-лимфоциты не инфицированные и инфицированные, соответственно, ${{T}_{2}}$, $T_{2}^{ * }$ – не инфицированные и инфицированные макрофаги, $V$ – вирус, не проникший в клетку, $E$ – иммунные эффекторы, ${{\alpha }_{b}} = (T_{1}^{ * } + T_{2}^{ * }){{(T_{1}^{ * } + T_{2}^{ * } + {{K}_{b}})}^{{ - 1}}}$, ${{\alpha }_{d}} = (T_{1}^{ * } + T_{2}^{ * }){{(T_{1}^{ * } + T_{2}^{ * } + {{K}_{d}})}^{{ - 1}}}$. В качестве начальных данных для математической модели (5) был выбран вектор
(6)
${{T}_{1}}(0) = 500\,000,\quad {{T}_{2}}(0) = 4800,\quad T_{1}^{ * }(0) = 5000,\quad T_{2}^{ * }(0) = 10,\quad V(0) = 10\,000,\quad E(0) = 15,$
характеризующий состояние слабой инфекции пациента.

Математическая модель (5) содержит $19$ параметров, которые характеризуют особенности иммунитета и заболевания пациента [3], [5]. В статье [3] показано, что большинство из этих параметров определены достаточно точно по статистической информации, и только четыре из них $\varphi = ({{k}_{1}},{{k}_{2}},{{\lambda }_{1}},{{\lambda }_{2}})$ ($m = 4$) нуждаются в уточнении. Здесь ${{k}_{1}} = 8 \times {{10}^{{ - 7}}}$, ${{k}_{2}} = {{10}^{{ - 4}}}$ [мл/вирионы $ \cdot $ день] – скорости инфицирования Т-лимфоцитов и макрофагов, соответственно, ${{\lambda }_{1}} = {{10}^{4}}$, ${{\lambda }_{2}} = 31.98$ = 31.98 [кл/мл · день] – скорости производства Т-лимфоцитов и макрофагов (указаны усредненные значения параметров обратной задачи).

2.3. Модель распространения туберкулеза в особо эндемичных районах Российской Федерации с учетом лечения

Математическая модель распространения туберкулеза, основанная на законе баланса масс, учитывает динамику девяти групп в популяции [10]

$\begin{gathered} \dot {S} = \Pi + \varphi \mathcal{T} + {{\varphi }_{m}}{{\mathcal{T}}_{m}} - ({{\lambda }_{d}} + {{\lambda }_{{dm}}} + \mu )S, \\ {{{\dot {L}}}_{A}} = {{\lambda }_{d}}(S + {{L}_{B}} + {{L}_{{Bm}}}) - (\varepsilon + k + \mu ){{L}_{A}}, \\ {{{\dot {L}}}_{{Am}}} = {{\lambda }_{{dm}}}(S + {{L}_{B}} + {{L}_{{Bm}}}) - (\varepsilon + k + \mu ){{L}_{{Am}}}, \\ \end{gathered} $
(7)
$\begin{gathered} {{{\dot {L}}}_{B}} = k{{L}_{A}} + \gamma I - ({{\lambda }_{d}} + {{\lambda }_{{dm}}} + \nu + \mu ){{L}_{B}}, \\ {{{\dot {L}}}_{{Bm}}} = k{{L}_{{Am}}} + \gamma {{I}_{m}} - ({{\lambda }_{d}} + {{\lambda }_{{dm}}} + \nu + \mu ){{L}_{{Bm}}}, \\ \dot {I} = \varepsilon {{L}_{A}} + \nu {{L}_{B}} + (1 - \eta )\omega \mathcal{T} - (\gamma + \delta + {{\mu }_{i}})I, \\ \end{gathered} $
$\begin{gathered} {{{\dot {I}}}_{m}} = \varepsilon {{L}_{{Am}}} + \nu {{L}_{{Bm}}} + \eta \omega \mathcal{T} + \omega {{\mathcal{T}}_{m}} - (\gamma + {{\delta }_{m}} + {{\mu }_{i}}){{I}_{m}}, \\ \dot {\mathcal{T}} = \delta I - (\varphi + \omega + {{\mu }_{t}})\mathcal{T}, \\ {{{\dot {\mathcal{T}}}}_{m}} = {{\delta }_{m}}{{I}_{m}} - ({{\varphi }_{m}} + \omega + {{\mu }_{t}}){{\mathcal{T}}_{m}},\quad t \in ({{t}_{0}},T). \\ \end{gathered} $
Здесь $S$ – чувствительные, вакцинированные от туберкулеза, $L$ – носители латентной инфекции (около 90% всего населения),  больные с открытой  формой болезни: $I$ – без лечения, $\mathcal{T}$ – находящиеся на лечении. Индекс $A$ характеризует быстрое развитие активной формы болезни, а индекс $B$ – медленное, индекс $m$ обозначает наличие штамма множественной лекарственной устойчивости к микобактериям туберкулеза (МЛУ-ТБ), ${{\lambda }_{d}} = \chi \beta (I(t) + o\mathcal{T}(t)){\text{/}}N$, ${{\lambda }_{{dm}}} = \chi {{\beta }_{m}}({{I}_{m}}(t) + o{{\mathcal{T}}_{m}}(t)){\text{/}}N$, $N$ – вся популяция. В условиях Сибирского федерального округа начальные данные рассматривались за ${{t}_{0}} = 2009$ (в тыс. человек):
$\begin{gathered} S = 249,\quad {{L}_{A}} = 8500,\quad {{L}_{{{{A}_{m}}}}} = 1150,\quad {{L}_{B}} = 8656,\quad {{L}_{{{{B}_{m}}}}} = 1189, \\ I = 34,\quad {{I}_{m}} = 4,\quad \mathcal{T} = 51,\quad {{\mathcal{T}}_{m}} = 6, \\ \end{gathered} $
естественный прирост $\Pi = 266$ тыс. чел/год, средняя продолжительность жизни $1{\text{/}}\mu = 62.5$ лет, $\beta = 0.025$.

Здесь в рамках модели (1) $y = (S,{{L}_{A}},{{L}_{{{{A}_{m}}}}},{{L}_{B}},{{L}_{{{{B}_{m}}}}},I,{{I}_{m}},\mathcal{T},{{\mathcal{T}}_{m}})$ ($n = 9$), $\varphi = \left( {\varepsilon ,k,\nu ,\delta ,{{\delta }_{m}},\chi ,o} \right)$ ($m = 7$), усредненные значения которых следующие: $\varepsilon = 0,129$ и $k = 0.821$ [год] – скорости раннего и перехода к позднему прогрессированию болезни, соответственно, $\nu = 0.075$ [год] – скорость развития активной формы болезни при эндогенной активации, $\delta = 0.72$ и ${{\delta }_{m}} = 0$ [чел/год] – скорости обнаружения индивидов с активной формой туберкулеза без и со штаммом МЛУ-ТБ, соответственно, $\chi = 0.343$ – параметр частичного иммунитета и $o = 0.6$ – параметр заразности индивидов во время лечения. Значения остальных параметров считаются известными [10].

3. О КОРРЕКТНОСТИ ПРЯМОЙ ЗАДАЧИ

Прямая задача (1), (2) сводится к системе нелинейных интегральных уравнений

(8)
$y(t) = \psi + \int\limits_{{{t}_{0}}}^t {({{\Phi }_{\tau }}y)} d\tau .$
Здесь ${{\Phi }_{\tau }}y = \Phi (y(\tau ),\tau ,q)$. Будем говорить, что $y \in {{L}_{2}}(S;X)$, $S = [{{t}_{0}},T]$, если для всех $j \in {{I}_{n}}$ функция ${{y}_{j}} \in {{L}_{2}}(S)$ со значениями в некотором вещественном банаховом пространстве $X$.

Определение 1. Оператор ${{\Phi }_{t}} \in \left( {{{L}_{2}}(S;X) \to {{L}_{2}}(S;X)} \right)$ называется оператором Вольтерра, если для любой пары $g,r \in {{L}_{2}}(S;X)$ из соотношения $g(\lambda ) = r(\lambda )$, выполняющегося для всех $\lambda \in {{S}_{\tau }}: = [{{t}_{0}},\tau ]$, $\tau \in S$, следует, что для всех $\lambda \in {{S}_{\tau }}$ имеет место $({{\Phi }_{t}}g)(\lambda ) = ({{\Phi }_{t}}r)(\lambda )$.

Определение 2. Семейство операторов ${{\Phi }_{t}} \in \left( {{{L}_{2}}(S;X) \to {{L}_{2}}(S;X)} \right)$, $t \in S$, называется ограниченно липшиц-непрерывным, если существует вещественная функция $\mu ({{y}_{1}},{{y}_{2}})$, возрастающая по ${{y}_{1}}$ и ${{y}_{2}}$ и такая, что для всех $t \in S$ и ${{y}_{1}},{{y}_{2}} \in {{L}_{2}}(S;X)$ справедлива оценка

${{\left\| {{{\Phi }_{t}}{{y}_{1}} - {{\Phi }_{t}}{{y}_{2}}} \right\|}_{{{{L}_{2}}(S;X)}}} \leqslant \mu (y_{1}^{ * },y_{2}^{ * }){{\left\| {{{y}_{1}} - {{y}_{2}}} \right\|}_{{{{L}_{2}}(S;X)}}}.$
Здесь $y_{j}^{ * } = {{\left\| {{{y}_{j}}} \right\|}_{{{{L}_{2}}(S;X)}}}$, $j = 1,2$.

Определение 3. Семейство операторов ${{{\text{\{ }}{{\Phi }_{t}}{\text{\} }}}_{{t \in S}}}$ принадлежит классу $\mathcal{K}(T,\mu )$, если для всех $t \in S$ оператор ${{\Phi }_{t}}$ есть оператор Вольтерра, ${{\Phi }_{t}}0 = 0$ и семейство ${{{\text{\{ }}{{\Phi }_{t}}{\text{\} }}}_{{t \in S}}}$ ограниченно липшиц-непрерывно.

Справедливы следующие утверждения относительно корректности задачи (8) [15], [16]:

Теорема 1 (о корректности в малом). Предположим, что семейство операторов ${{{\text{\{ }}{{\Phi }_{t}}{\text{\} }}}_{{t \in S}}}$ принадлежит классу $\mathcal{K}(T,\mu )$. Тогда существует ${{T}_{ \star }} \in ({{t}_{0}},T)$ такое, что для любого $t \in ({{t}_{0}},{{T}_{ \star }})$ решение операторного уравнения (8) существует, единственно в классе ${{L}_{2}}({{S}_{t}};X)$, а также непрерывно зависит от данных.

Теорема 2 (о корректности в окрестности точного решения). Предположим, что для некоторого $\psi \in {{\mathbb{R}}^{n}}$ существует $y \in {{L}_{2}}(S;X)$. Предположим, что семейство операторов ${{{\text{\{ }}{{\Phi }_{t}}{\text{\} }}}_{{t \in S}}}$ принадлежит классу $\mathcal{K}(T,\mu )$. Тогда найдется $\delta \in {{\mathbb{R}}_{ + }}$ такое, что для любого ${{\psi }^{\delta }} \in O(\psi ,\delta ): = {\text{\{ }}\phi :{{\left\| {\phi - \psi } \right\|}_{{{{\mathbb{R}}^{n}}}}} < \delta {\text{\} }}$ решение уравнения (8) в классе ${{L}_{2}}({{S}_{t}};X)$ существует, единственно и непрерывно зависит от данных.

4. ИДЕНТИФИЦИРУЕМОСТЬ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ И r-РЕШЕНИЕ

Отметим, что в некоторых случаях решение обратной задачи (1)–(3) неединственно, однако можно указать комбинации неизвестных параметров, которые определяются однозначно. С этой целью исследуется идентифицируемость математических моделей [11], [17].

Определение 4. Модель (1), (2) называется априорно идентифицируемой, если ее параметры $q \in {{\mathbb{R}}^{m}}$ можно однозначно определить по данным измерений (3).

Условно разделяют два основных вида идентифицируемости.

1. Априорная (структурная) идентифицируемость используется при изучении структуры модели без зашумленных данных. К методам структурной идентифицируемости относятся методы придаточной функции, дифференциальной алгебры, разложения в ряд и другие. Более детальное описание методов можно найти в работах [18], [19].

2. Практическая (апостериорная) идентифицируемость позволяет выявить идентифицируемые параметры относительно зашумленных данных методами Монте Карло, корреляционной матрицы, DAISY и других [20].

Для оценки чувствительности параметров модели относительно погрешностей в измерениях (данных модели) используется анализ чувствительности, основанный на изучении свойств матрицы чувствительности системы. Методы анализа чувствительности требуют знания точных параметров модели, расположение и количество измерений. К методам чувствительности относятся ортогональный метод, метод собственных значений, корреляционный метод и другие [21].

Понятие идентифицируемости тесно связано с существованием $r$-решения компактного оператора или системы линейных алгебраических уравнений [22], [23]. Действительно, применив к исследуемым обратным задачам (1)–(3) процедуру линеаризации относительно параметров q (для модели (5) она подробно описана в работе [24]), получим систему линейных алгебраических уравнений вида

$Aq = f.$
Здесь $f = {{{\text{\{ }}{{f}_{{ik}}}{\text{\} }}}_{\substack{ i \in I \\ k \in {{I}_{K}} } }}$ – вектор измерений, $A$ – матрица линеаризованной обратной задачи. Применяя метод сингулярного разложения к произвольным прямоугольным матрицам $A = U\Sigma V{\text{*}}$, где $U$ и $V$ – унитарные матрицы, получаем
$U\Sigma V{\text{*}}q = f \to \Sigma V{\text{*}}q = U{\text{*}}f,$
где $\Sigma = \operatorname{diag} \{ {{\sigma }_{1}},\; \ldots ,\;{{\sigma }_{p}}\} $, ${{\sigma }_{j}}$ – сингулярные числа матрицы $A$. Обозначим $g = U{\kern 1pt} {\text{*}}f$, $z = V{\text{*}}q$, тогда решение системы алгебраических уравнений с точностью до унитарного преобразования можно выписать в следующем виде:
(9)
$\begin{array}{*{20}{c}} {{{z}_{1}} = \frac{{{{g}_{1}}}}{{{{\sigma }_{1}}}},\quad {{z}_{2}} = \frac{{{{g}_{2}}}}{{{{\sigma }_{2}}}},\quad \ldots ,\quad {{z}_{r}} = \frac{{{{g}_{r}}}}{{{{\sigma }_{r}}}},} \\ {{{z}_{i}} = 0,\quad {\text{если}}\quad {{\sigma }_{i}} < \delta ,\quad i = r + 1,\; \ldots ,\;p.} \end{array}$
Здесь $r$ – натуральное число обрыва последовательности малых сингулярных чисел, т.е. тех ${{\sigma }_{i}}$, которые меньше некоторого заданного $\delta $ (в [22], [23] выбор $\delta $ связан с уровнем ошибки в измерениях $f$).

Полученная последовательность $q = Vz$ решения задачи $Aq = f$, в которой $r$ ненулевых компонент, аналогична последовательности идентифицируемых параметров для линеаризованных моделей.

5. ВАРИАЦИОННАЯ ПОСТАНОВКА ОБРАТНОЙ ЗАДАЧИ

Решение обратной задачи (1)–(3) будем искать, минимизируя целевой функционал:

(10)
$J(q) = \gamma \sum\limits_{k = 1}^K {\sum\limits_{i \in I} {{{{\left( {{{y}_{i}}({{t}_{k}};q) - {{f}_{{ik}}}} \right)}}^{2}}} } .$
Здесь $\gamma = (T - {{t}_{0}}){\text{/}}K$.

5.1. Локализация точек глобального минимума

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

В качестве альтернативного метода упомянутым экристическим методам глобальной многопараметрической оптимизации мы использовали работу [8], в которой предложен полуаналитический метод тензорного разложения, состоящий в представлении многопараметрической функции $J(q)$ в виде тензора с последующим разложением в тензорный поезд. Важное преимущество метода тензорного разложения заключается в возможности многоуровнего распараллеливания алгоритма.

5.2. Градиентный метод

Для минимизации функционала (10) используем метод минимальных ошибок, который заключается в определении приближенного решения по формуле [26]

(11)
${{q}_{{n + 1}}} = {{q}_{n}} - {{\alpha }_{n}}J{\kern 1pt} {\text{'}}({{q}_{n}}),\quad {{\alpha }_{n}} = 2J({{q}_{n}}){\text{/}}{{\left\| {J{\kern 1pt} {\text{'}}({{q}_{n}})} \right\|}^{2}},\quad {{q}_{0}} \in {{\mathbb{R}}^{{M + N}}}.$
Здесь ${{\alpha }_{n}}$ – параметр спуска, $J{\kern 1pt} {\text{'}} \in {{\mathbb{R}}^{{{{M}_{1}}}}}$ – градиент целевого функционала (10) [27].

Лемма. Градиент $J{\kern 1pt} {\text{'}}(q)$ целевого функционала (10) удовлетворяет соотношению

(12)
$J{\kern 1pt} {\text{'}}(q) = \mathop {\left( { - \int\limits_{{{t}_{0}}}^T {\Phi _{q}^{{\text{т}}}} (y(t),t,q)\Psi (t)dt, - \Psi ({{t}_{0}})} \right)}\nolimits^{\text{т}} .$
Здесь $\Psi (t) = {{({{\Psi }_{1}}(t), \ldots ,\;{{\Psi }_{N}}(t))}^{{\text{т}}}}$, ${{\Psi }_{j}}$ – кусочно-гладкие функции на интервалах $({{t}_{k}},{{t}_{{k + 1}}})$, $k = 0,\; \ldots ,\;K$, $j \in \mathcal{J}$, которые являются решением сопряженной задачи
$\dot {\Psi } = - \Phi _{y}^{{\text{т}}}(y(t),t,q)\Psi (t),\quad t \in \bigcup\limits_{k = 0}^K {({{t}_{k}},{{t}_{{k + 1}}})} ,{{t}_{{K + 1}}} = T;$
(13)
$\Psi (T) = 0,\quad {{[{{\Psi }_{i}}]}_{{t = {{t}_{k}}}}} = 2\gamma ({{y}_{i}}({{t}_{k}};q) - {{f}_{{ik}}}),\quad k \in {{I}_{K}},\quad i \in I,$
${{\Phi }_{y}}(y(t),t,q) \in {{\mathbb{R}}^{N}} \times {{\mathbb{R}}^{N}}\quad {\text{и}}\quad {{\Phi }_{q}}(y(t),t,q) \in {{\mathbb{R}}^{N}} \times {{\mathbb{R}}^{M}}$
суть матрицы Якоби, ${{[{{\Psi }_{i}}]}_{{t = {{t}_{k}}}}} = {{\Psi }_{i}}({{t}_{k}} + 0) - {{\Psi }_{i}}({{t}_{k}} - 0)$ – скачок функции ${{\Psi }_{i}}$ в точке ${{t}_{k}}$.

Доказательство. Согласно определению производной Фреше функции $J(q)$

$\delta J = \left\langle {J{\kern 1pt} {\text{'}}(q),\delta q} \right\rangle + o(\left| {\delta q} \right|).$
Здесь $\delta J = J(q + \delta q) - J(q)$, $\left\langle { \cdot , \cdot } \right\rangle $ – скалярное произведение в векторном пространстве.

Используя (10), запишем выражение

$\begin{gathered} \delta J = 2\gamma \sum\limits_{k = 1}^K {\sum\limits_{i \in I} {({{y}_{i}}(t;q) - {{f}_{{ik}}})\delta {{y}_{i}}({{t}_{k}};\delta q)} } + \gamma \sum\limits_{k = 1}^K {\sum\limits_{i \in I} {\delta y_{i}^{2}} } ({{t}_{k}};\delta q) = \\ = \;2\gamma \sum\limits_{k = 1}^K {\sum\limits_{i \in I} {\int\limits_{{{t}_{0}}}^T {({{y}_{i}}(t;q) - {{f}_{{ik}}})} } } {\kern 1pt} \delta {{y}_{i}}(t;\delta q)\delta (t - {{t}_{k}})dt + \gamma \sum\limits_{k = 1}^K {\sum\limits_{i \in I} {\delta y_{i}^{2}} } ({{t}_{k}};\delta q) = {{A}_{1}} + {{A}_{2}}. \\ \end{gathered} $

Введем задачу Коши для функции $\delta y(t;\delta q) = y(t;q + \delta q) - y(t;q)$

(14)
$\mathcal{L}\delta y = 0,\quad \delta y({{t}_{0}}) = \delta \psi ,\quad t \in ({{t}_{0}},T).$
Здесь

$\mathcal{L}\delta y = \tfrac{{d\delta y}}{{dt}} - {{\Phi }_{y}}(y(t),t,\phi )\delta y - {{\Phi }_{\phi }}(y(t),t,\phi )\delta \phi .$

Распишем скалярное произведение $\left\langle {\mathcal{L}\delta y,\tilde {\Psi }} \right\rangle = 0$ в векторном пространстве ${{L}_{2}}$, где $\tilde {\Psi }$ – пробная вектор-функция:

(15)
$0 = \left\langle {\mathcal{L}\delta y,\tilde {\Psi }} \right\rangle = \left\langle {\frac{{d\delta y}}{{dt}},\tilde {\Psi }} \right\rangle - \left\langle {{{\Phi }_{y}}(y(t),t,\phi )\delta y,\tilde {\Psi }} \right\rangle - \left\langle {{{\Phi }_{\phi }}(y(t),t,\phi )\delta \phi ,\tilde {\Psi }} \right\rangle .$

Введем обобщенную задачу Коши для сопряженной функции $\tilde {\Psi }(t)$:

$\begin{gathered} \frac{{d\tilde {\Psi }}}{{dt}} = - \Phi _{y}^{{\text{т}}}(y(t),t,\phi )\tilde {\Psi }(t) + 2\gamma \sum\limits_{k = 1}^K {\int\limits_{{{t}_{0}}}^T {(y(t;q) - {{f}_{k}})} } \delta (t - {{t}_{k}})dt,\quad t \in \mathbb{R}; \\ \tilde {\Psi }(T) = 0. \\ \end{gathered} $
Здесь ${{f}_{k}} \in {{\mathbb{R}}^{N}}$ состоит из 0 и компонент ${{f}_{{ik}}}$, $i \in I$. Отметим, что обобщенная задача Коши может быть сведена к задаче Коши (13) [27].

Тогда применяя к выражению (15) метод интегрирования по частям и учитывая, что

${{A}_{1}} = \left\langle {2\gamma \sum\limits_{k = 1}^K {\int\limits_{{{t}_{0}}}^T {(y(t;q) - {{f}_{k}})} \delta } (t - {{t}_{k}})dt,\delta y(t;\delta q)} \right\rangle ,$
получим

${{A}_{1}} = - \left\langle {\delta \phi ,\int\limits_{{{t}_{0}}}^T {\Phi _{\phi }^{{\text{т}}}} (y(t),t,\phi )\tilde {\Psi }dt} \right\rangle - \left\langle {\delta \psi ,\tilde {\Psi }({{t}_{0}})} \right\rangle .$

Покажем, что ${{A}_{2}} = o(\left| {\delta q} \right|)$. Введем обозначения:

$B(t) = exp\left\{ {\int\limits_{{{t}_{0}}}^t {{{\Psi }_{y}}} (y(\tau ),\tau ,\phi )d\tau } \right\},\quad D(t) = exp\left\{ { - \int\limits_{{{t}_{0}}}^t {{{\Psi }_{y}}} (y(\tau ),\tau ,\phi )d\tau } \right\}.$
Тогда решение задачи Коши (14) может быть представлено в виде

$\delta y = B(t)\left( {\delta \psi - \int\limits_{{{t}_{0}}}^t D (\tau ){{\Phi }_{\phi }}(y(\tau ),\tau ,\phi )\delta \phi d\tau } \right).$

Применяя неравенства треугольника и Гёльдера, получаем

${{\left| {\delta y} \right|}^{2}} \leqslant {{\left\| {B(t)} \right\|}^{2}}{{\left| {\delta \psi } \right|}^{2}} + {{\left\| {B(t)} \right\|}^{2}}\mathop {\left| {\int\limits_{{{t}_{0}}}^t D (\tau ){{\Phi }_{\phi }}(y(\tau ),\tau ,\phi )\delta \phi d\tau } \right|}\nolimits^2 \leqslant {{\left\| {B(t)} \right\|}^{2}}{{\left| {\delta \psi } \right|}^{2}} + C(t){{\left| {\delta \phi } \right|}^{2}}.$
Здесь
$C(t) = \sqrt {t - {{t}_{0}}} {{\left\| {B(t)} \right\|}^{2}}\mathop {\left( {\int\limits_{{{t}_{0}}}^t {{{{\left\| {D(\tau )} \right\|}}^{2}}} d\tau } \right)}\nolimits^{1/2} \mathop {\left( {\int\limits_{{{t}_{0}}}^t {{{{\left\| {{{\Phi }_{\phi }}(y(\tau ),\tau ,\phi )} \right\|}}^{2}}d\tau } } \right)}\nolimits^{1/2} .$
Тогда верны оценки для ${{A}_{2}}$:

${{A}_{2}} = \gamma \sum\limits_{k = 1}^K {\sum\limits_{i \in I} {\delta y_{i}^{2}} } ({{t}_{k}};\delta q) \leqslant \gamma \int\limits_{k = 1}^K {{{{\left| {\delta y({{t}_{k}};\delta q)} \right|}}^{2}}} \leqslant \gamma \sum\limits_{k = 1}^K {\left( {{{{\left\| {B({{t}_{k}})} \right\|}}^{2}}{{{\left| {\delta \psi } \right|}}^{2}} + C({{t}_{k}}){{{\left| {\delta \phi } \right|}}^{2}}} \right)} .$

Введем постоянные

$\alpha = \gamma \sum\limits_{k = 1}^K {{{{\left\| {B({{t}_{k}})} \right\|}}^{2}}} ,\quad \beta = \gamma \sum\limits_{k = 1}^K C ({{t}_{k}}).$
Тогда ${{A}_{2}} = o(\left| {\delta q} \right|)$, так как по определению малых величин справедливо

$\mathop {lim}\limits_{|\delta \phi |,|\delta \psi | \to 0} \frac{{\alpha {{{\left| {\delta \psi } \right|}}^{2}} + \beta {{{\left| {\delta \phi } \right|}}^{2}}}}{{\sqrt {{{{\left| {\delta \phi } \right|}}^{2}} + {{{\left| {\delta \psi } \right|}}^{2}}} }} = 0.$

Вспоминая определение производной Фреше и выражения для $\delta J$ и ${{A}_{1}}$, получим градиент $J{\kern 1pt} {\text{'}}(q)$ целевого функционала (10), удовлетворяющий выражению (12).

Отметим, что сходимость градиентного метода сильно зависит от выбора начального приближения ${{q}_{0}}$. Поэтому на первом этапе используем статистические методы (имитации отжига, генетического алгоритма и т.п.) для определения области глобального минимума [24], [25], [28], а затем в данной области уточняем решение обратной задачи (1)–(3). В приведенных численных расчетах начальное приближение ${{q}_{0}}$ находится в области решения обратной задачи.

6. ЧИСЛЕННЫЕ РАСЧЕТЫ

В следующих разделах приведены результаты численного решения обратной задачи (1)–(3) для моделей из разд. 2 в случае не зашумленных синтетических данных (3), которые получаются при решении прямой задачи (1), (2) при некоторых усредненных значениях параметров $q$, указанных в разд. 2 (считаем, что начальные данные $\psi $ известны). Задаем начальное приближение ${{q}_{0}}$ градиентного метода, находящееся в области глобального минимума задачи минимизации $J(q)$. Предполагая, что ${{q}_{n}}$ уже известно, решаем прямую задачу (1), (2) и определяем $y({{t}_{k}};{{q}_{n}})$. Затем численно решаем сопряженную задачу (13), вычисляем градиент (12) целевого функционала и определяем следующее приближение ${{q}_{{n + 1}}}$ согласно (11). Алгоритм повторяем до тех пор, пока не достигнуто условие остановки $J({{q}_{n}}) < {{10}^{{ - 3}}}$. Отметим, что все задачи определения вектора параметров по дополнительным измерениям некоторых функций в фиксированные моменты времени, изученные в данном разделе, являются априорно идентифицируемыми [19].

6.1. Математическая модель инфекционного заболевания

В численных расчетах использовались синтетические данные (3) при ${{I}_{4}} = {\text{\{ }}1,\;2,\;3,\;4{\text{\} }}$, $K = 5$, равномерно распределенные на отрезке $[0,T] = [0,40]$ дней, то есть ${{t}_{k}} = - 3 + hk$, $k = 1,\; \ldots ,\;K$, $h = 8$ сут. Шаг расчетной сетки при решении прямой и сопряженной задач ${{h}_{t}} = 0.01$. Для определения вектора восьми параметров $q = \left( {\tau ,\alpha ,\beta ,\gamma ,\rho ,\eta ,{{\mu }_{c}},{{\mu }_{f}}} \right)$ используется метод минимальных ошибок с градиентом целевого функционала $J{\kern 1pt} {\text{'}} \in {{\mathbb{R}}^{8}}$, удовлетворяющий (12), (13).

График зависимости относительной ошибки $e({{q}_{n}}) = \sqrt {{{{\left| {{{q}_{n}} - q} \right|}}^{2}}{\text{/}}{{{\left| q \right|}}^{2}}} $ от общего числа итераций $n = 5.6 \times {{10}^{6}}$ приведен на фиг. 1a. Погрешность восстановления восьми параметров составила 5.5%, т.е. $e({{q}_{n}}) < 0.055$. Для получения большей точности необходимо провести больше итерация градиентного метода.

Фиг. 1.

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

Отметим, что модель (4) получила развитие как для вирусных заболеваний (пневмония [4], гепатит В [6]), так и для бактериальных (туберкулез [9]) и вирусно-бактериальных [2].

6.2. Модель динамики ВИЧ

Для численного решения обратной задачи восстановления вектора $q = ({{k}_{1}},{{k}_{2}},{{\lambda }_{1}},{{\lambda }_{2}})$ для математической модели (5) использовались синтетические данные (3) о трех концентрациях:

(16)
${{T}_{1}}({{t}_{k}}) + T_{1}^{ * }({{t}_{k}}) = {{\Phi }_{1}}({{t}_{k}}),\quad V({{t}_{k}}) = {{\Phi }_{2}}({{t}_{k}}),\quad E({{t}_{k}}) = {{\Phi }_{3}}({{t}_{k}}),\quad k \in {{I}_{K}},$
взятые $K = 14$ раз равномерно в течение $T = 100$ дней (один раз в неделю). Расчетная область $[0,T]$ покрывалась равномерной сеткой с шагом 0.01.

График зависимости относительной ошибки $e({{q}_{n}})$ от общего числа итераций $n = 1400$ приведен на фиг. 1б. Погрешность восстановления четырех параметров составила 0.2%, т.е. $e({{q}_{n}}) < 0.002$.

6.3. Модель распространения туберкулеза в Сибирском федеральном округе с учетом лечения

Для численного решения обратной задачи для математической модели (7) использовались синтетические данные (3) о трех группах в популяции (${{I}_{3}} = {\text{\{ }}1,\;8,\;9{\text{\} }}$)

(17)
$S({{t}_{k}}) = {{S}_{k}},\quad \mathcal{T}({{t}_{k}}) = {{\mathcal{T}}_{k}},\quad {{\mathcal{T}}_{m}}({{t}_{k}}) = {{\mathcal{T}}_{{{{m}_{k}}}}},\quad k \in {{I}_{K}},$
за $K = 4$ года для определения семи параметров модели $q = \mathop {\left( {\varepsilon ,k,\nu ,\delta ,{{\delta }_{m}},\chi ,o} \right)}\nolimits^{\text{т}} $. Шаг расчетной сетки при решении прямой и сопряженной задач ${{h}_{t}} = 0.001$ в области $[{{t}_{0}},T] = [2009,2013]$.

В результате применения метода минимальных ошибок вектор параметров восстановлен с погрешностью 30%, т.е. $e({{q}_{n}}) < 0.3$, за $n = 1737440$ итераций метода минимальных ошибок (см. фиг. 1в). Такая погрешность объясняется некорректностью задачи определения семи параметров $q = \mathop {\left( {\varepsilon ,k,\nu ,\delta ,{{\delta }_{m}},\chi ,o} \right)}\nolimits^{\text{т}} $ по 12 точкам измерения некоторых функций (17) нелинейной модели. Для достижений лучших результатов необходимо сузить область поиска вектора параметров модели, а именно выбрать более точное начальное приближение ${{q}_{0}}$ градиентного метода с помощью глобальных методов оптимизации.

7. ЗАКЛЮЧЕНИЕ

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

Схожие структуры моделей и постановки обратных задач возникают в фармакокинетике при анализе распределения лекарственных препаратов в организме [29], экономике при описании поведения опционов и ценных бумаг [30], а также при анализе социальных процессов, описывающих распространение информации в онлайн социальных сетях [31].

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

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

  2. Бочаров Г.А., Марчук Г.И. Прикладные проблемы математического моделирования в иммунологии // Ж. вычисл. матем. и матем. физ. 2000. Т. 40. № 12. С. 1905–1920.

  3. Callaway D.S., Perelson A.S. HIV-I infection and low steady state viral loads // Bull. Math. Biol. 2001. V. 64. P. 29–64.

  4. Романюха А.А., Руднев С.Г. Вариационный принцип в исследовании противоинфeкционного иммунитета на примере пневмонии // Матем. моделирование. 2001. Т. 13. № 8. С. 65–84.

  5. Adams B.M., Banks H.T., Davidian M., Kwon H.-D., Tran H.T., Wynne S.N., Rosenberg E.S. HIV dynamics: Modeling, data analysis, and optimal treatment protocols // J. Comput. Appl. Math. 2005. V. 184. P. 10–49.

  6. Руднев С.Г., Романюха А.А., Яшин А.И. Моделирование развития Т-системы иммунитета и оценка эффективности распределения ресурсов // Матем. моделирование. 2007. Т. 19. № 1. С. 25–42.

  7. Bocharov G., Kim A., Krasovskii A., Chreshnev V., Glushenkova V., Ivanov A. An extremal shift method for control of HIV infection dynamics // Russ. J. Numer. Anal. Math. Modelling. 2015. V. 30. № 1. P. 11–25.

  8. Zheltkova V.V., Zheltkov D.A., Grossman Z., Bocharov G.A., Tyrtyshnikov E.E. Tensor based approach to the numerical treatment of the parameter estimation problems in mathematical immunology // J. Inverse Ill-Posed Probl. 2018. V. 26. № 1. P. 51–66.

  9. Авилов К.К., Романюха А.А. Математические модели распространения и контроля туберкулеза // Математическая биология и биоинформатика. 2007. Т. 2. № 2. С. 188–318.

  10. Trauer J.M., Denholm J.T., McBryde E.S. Construction of a mathematical model for tuberculosis transmission in highly endemic regions of the Asia-pacific // J. Theor. Biol. 2014. V. 358. P. 74–84.

  11. Cobelli C., Carson E. Introduction to Modeling in Physiology and Medicine. The Netherlands: Academic Press, 2008.

  12. Кряжимский А.В., Осипов Ю.С. Устойчивые решения обратных задач динамики управляемых систем // Тр. МИАН СССР. 1988. Т. 185. С. 126–146.

  13. Красовский Н.Н., Субботин А.И. Позиционные дифференциальные игры. М.: Наука, 1974.

  14. Кряжимский А.В., Осипов Ю.С. О разрешимости задач гарантирующего управления для частично наблюдаемых линейных динамических систем // Тр. МИАН СССР. 2012. Т. 277. С. 152–167.

  15. Кабанихин С.И. Регуляризация операторного уравнения Вольтерра первого рода с ограниченно липшиц-непрерывным ядром // Докл. АН СССР. 1989. Т. 39. № 3. С. 549–552.

  16. Kabanikhin S.I. Definitions and examples of inverse and ill-posed problems // J. Inverse Ill-Posed Probl. 2008. V. 16. P. 317–357.

  17. Bellman R., Astrom K. On structural identifiability // Math. Biosciences. 1970. V. 30. Issue 4. P. 65–74.

  18. Bellu G., Saccomani M., Audoly S., D’Angiò L. DAISY: a new software tool to test global identifiability of biological and physiological systems // Comput. Methods Programs Biomed. 2007. V. 88. Issue 1. P. 52–61.

  19. Кабанихин С.И., Воронов Д.А., Гродзь А.А., Криворотько О.И. Идентифицируемость математических моделей медицинской биологии // Вавиловский ж. генетики и селекции. 2015. Т. 19. № 6. С. 738–744.

  20. Miao H., Xia X., Perelson A.S., Wu H. On Identifiability of nonlinear ODE models and applications in viral dynamics // SIAM Rev. Soc. Ind. Appl. Math. 2011. V. 53. Issue 1. P. 3–39.

  21. Quaiser T., Mönnigmann M. Systematic identifiability testing for unambiguous mechanistic modeling – application to JAK-STAT, MAP kinase, and NF-B signaling pathway models // BMC Systems Biology. 2009. P. 3–50.

  22. Годунов C.К., Антонов А.Г., Кирилюк О.П., Костин В.И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. Новосибирск: Наука, 1992.

  23. Чеверда В.А., Костин В.И. $R$-псевдообратный для компактного оператора // Сиб. электрон. матем. изв. 2010. Т. 7. С. 258–282.

  24. Banks H.T., Kabanikhin S.I., Krivorotko O.I., Yermolenko D.V. A numerical algorithm for constructing an individual mathematical model of HIV dynamics at cellular level // J. Inverse Ill-Posed Probl. 2018. V. 26. № 6. P. 859–873.

  25. Kabanikhin S., Krivorotko O., Kashtanova V. A combined numerical algorithm for reconstructing the mathematical model for tuberculosis transmission with control programs // J. Inverse Ill-Posed Probl. 2018. V. 26. Issue 1. P. 121–131.

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

  27. Kabanikhin S.I., Krivorotko O.I. Identification of biological models described by systems of nonlinear differential equations // J. Inverse Ill-Posed Probl. 2015. V. 23. Issue 5. P. 519–527.

  28. Kabanikhin S.I., Krivorotko O.I., Ermolenko D.V., Kashtanova V.N., Latyshenko V.A. Inverse problems of immunology and epidemiology // Eurasian Journal of Math. and Computer Applications. 2017. V. 5. Issue 2. P. 14–35.

  29. Ilyin A., Kabanikhin S., Nurseitov D., Nurseitova A., Asmanova N., Voronov D., Bakytov D. Analysis of ill-posedness and numerical methods of solving a nonlinear inverse problem in pharmacokinetics for the two-compartmental model with extravascular drug administration // J. Inverse Ill-Posed Probl. 2012. V. 20. Issue 1. P. 39–64.

  30. Solow R.M. A contribution to the theory of economic growth // The Quarterly Journal of Economics. 1956. V. 70. Issue 1. P. 65–94.

  31. Kabanikhin S., Krivorotko O., Zhang S., Kashtanova V., Wang Yu. Tensor train optimization for mathematical model of social networks // arXiv:1906.05246 [math.OC]. 2019.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал