CC BY
31
Наука к Образование
МГТУ им. Н.Э. Баумана
Сетевое научное издание
ISSN 1994-0408
Ссылка на статью:
// Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. №11. С. 406-425.
БОТ: 10.7463/1115.0826730
Представлена в редакцию: 6.12.2015 © МГТУ им. Н.Э. Баумана
УДК 51.76: 517.9 : 57.085.23
Получение точечных оценок и законов распределения вероятностей параметров математической модели развития клеточной популяционной системы с учетом контактного торможения
Виноградова М. С.*' *mmvinogradova@rambler.ru
Для математической модели, описывающей динамику развития клеточной популяционной системы, состоящей из нормальных и аномальных стволовых клеток человека, предложена методика построения оценок и законов распределения параметров модели на основе ограниченных выборок экспериментальных данных. Модель представляет собой нелинейную систему обыкновенных дифференциальных уравнений первого порядка, определяющих изменения численности популяций нормальных и аномальных клеток. Переменными модели служат численности популяций. Модель учитывает ограниченность ресурсов, является нелинейной относительно переменных модели, но линейной относительно неизвестных параметров и рассматривается в области, где развиваются обе популяции. Для решения задачи параметрической идентификации математической модели был получен ее дискретный аналог. Задача идентификации векторов параметров модели была разбита на две независимые подзадачи. На первом этапе найдены оценки параметров, определяющих динамику развития популяции нормальных клеток, на втором ? оценки параметров, определяющих динамику развития популяции аномальных клеток. Особенности предложенной методики заключаются в выборе семейства непересекающихся шаблонов, которые позволяют получить как точечные оценки, так и апостериорные законы распределения параметров указанной математической модели при малых выборках. Получены соотношения, позволяющие вычислить точечные оценки и построить функции плотности распределения вероятностей векторов параметров модели на ограниченных выборках экспериментальных данных. При построении законов распределения параметров модели использовался байесовский подход и теория инвариантности Джеффриса. Полученные законы распределения параметров являются обобщенными распределениями Стьюдента.
Ключевые слова: математическая модель; стволовая клетка; параметрическая идентификация; байесовский подход; точечная оценка; функция плотности распределения
Для исследования динамики развития клеточных популяций в лабораторных условиях (in vitro) наряду с экспериментальными исследованиями активно используются методы математического моделирования [1, 2, 3].
1МГТУ им. Н.Э. Баумана, Москва, Россия
Введение
При экспериментальных исследованиях процесса культивирования обычно ограничиваются несколькими точками, характеризующими популяционную систему в начале культивирования, в середине и при завершении процесса [4, 5, 6], поскольку более детальные исследования обычно требует существенных затрат. Прогнозировать развитие клеточной популяции по таким данным достаточно затруднительно и в связи с этим математическое моделирование может играть значительную роль в анализе процессов, происходящих в клеточных популяциях in vitro.
Интерес к культивированию in vitro клеточных популяций и особенно к культивированию популяций стволовых клеток в последние десятилетия обусловлен активным развитием клеточной терапии [7, 8, 9]. Клеточный материал в необходимых количествах получают путем культивирования клеток, взятых у пациента или у донора. В силу естественной изменчивости клеток в процессе культивирования в клеточной популяции могут появиться клетки с геномными мутациями, в частности анеуплоидные клетки. Интерес к динамике развития популяционной системы включающей в себя анеуплоидные клетки вызван теорией [10, 11, 12, 13], связывающей развитие рака с наличием в организме значительного количества анеуплоидных клеток.
В работах [14, 15, 16, 17] были предложены и исследованы линейная и нелинейная модели развития клеточных популяций, состоящих нормальных (здоровых) и аномальных (анеуплоидных) клеток. Переменными модели являлись приведенные численности нормальных и аномальных клеток, в параметры модели входят естественные характеристики развития клеточной популяции, такие как средняя продолжительность клеточного цикла в популяции, доли клеток популяции, которые в течении этого цикла разделились в результате митоза, погибли в результате апоптоза или перешли из популяции здоровых клеток в популяцию аномальных клеток. Также в математических моделях учтено, что средние продолжительности клеточного цикла в популяции нормальных и аномальных клеток разные, то есть процессы развития популяций нормальных и аномальных клеток являются разнотемпо-выми. Линейная модель развития популяционной системы позволяет описывать поведение системы на ранних этапах культивирования, в условия неограниченности ресурсов. Методика идентификации параметров линейной модели описана в работе [18], там же был предложен подход к определению вероятностей реализации возможных сценариев развития клеточной популяционной системы. Базовая нелинейная модель учитывает ограниченность ресурсов.
Целью исследования является построение методики идентификации параметров нелинейной модели, которая позволяет не только получить точечные оценки параметров модели, но и законы их распределения.
Материал статьи организован следующим образом. Работа состоит из пяти разделов и заключения. В первом разделе приведена модель динамики развития клеточных популяций, учитывающая ограниченность ресурсов. Во втором обсуждается задача параметрическая идентификация математической модели, в третьем разделе описан процесс нахождения то-
чечных оценок параметров модели, в четвертом разделе приведена процедура проверки статистической гипотезы о виде закона распределения случайный возмущений модели, в пятом описан процесс построения функций плотности распределения вероятностей параметров математической модели.
1. Базовая нелинейная модель развития популяционной системы
Культивирования на поздних этапах происходит в условиях нехватки места на подложке. Известно, что при увеличении количества клеток на единице площади подложки скорость роста культуры клеток уменьшается, в случае возникновения контакта между клетками рост популяции может прекратиться [19]. Данный эффект получил название "контактного торможения". Применительно к исследуемой популяционной системе этот факт находит свое отражение в первую очередь в наличии функциональной зависимости доли клеток, разделившихся за среднее время клеточного цикла от численности популяций нормальных и аномальных клеток. Предложенная в [20] модель имеет вид
(1)
^ = ((! - Ао)(1 + Мо(х,у) - 2Мо(х,у)7о) - 1) х,
= Ц (-А + М 1(х,у) - АгИ 1(х,у)) у + 2Мо(х,у)(1 - Ао)7ох, Здесь:
М-7 (х, у), з = 1, 2, — функции, вычисляемые по формулам
0, ¡7(х,у) < 0; М7(х,у) = { 1, ¡7(х,у) > 1; з €{0, 1}; (2)
¡7(х, у), 0 <¡7(х, у) < 1,
где
¡7(х у) = а7 - в70х - 7у; (3)
т 0
ц = —т — параметр, введен для учета разнотемповости процессов деления в популяциях;
т1
т7 — средняя продолжительность клеточного цикла в з -й популяции клеток; А7 — доля клеток з-й популяции клеток, погибающих в результате апоптоза на временном интервале длительности т7, 0 < А7 < 1;
7о — доля здоровык клеток популяции нормальный клеток в рассматриваемой популяционной системе, переходящих на временном интервале длительности тх в процессе деления в популяцию аномальный клеток У, 0 < 7о < 1;
Рг7 — коэффициенты, характеризующие уменьшение доли делящихся клеток на временном интервале длительности т7 как за счет внутрипопуляционной (з = г), так и за счет межпопуляционной (з = г) конкуренции, > 0, г, з = 0, 1;
а7 — доля клеток з -й популяции, разделившихся в результате митоза на временном интервале длительности т7 в условиях отсутствия конкуренции, 0 < а7 < 1.
Модель представляет собой нелинейную систему обыкновенных дифференциальных уравнений первого порядка, определяющих изменения суммарной численности популяций нормальных и аномальных клеток. Несмотря на то, что формально полученная система определена на всей плоскости хОу, переменные х и у не могут принимать отрицательные значения, согласно биологическому смыслу моделируемых процессов.
В работе [17] было показано, что множество С = {(х, у) : х > 0, у > 0} является инвариантным компактом для динамической системы (1). Из-за кусочно-линейного характера соотношения (2) множество С необходимо разбить на четыре подмножества, на каждом из которых система (1) имеет разный вид.
Из четырех подмножеств множества С наибольший интерес представляет множество С1 = {(х, у): х > 0; у > 0; /о(х, у) > 0; Д(х, у) > 0}, так как в этом множестве размножаются оба типа клеток популяционной системы. С учетом (2) и (3) система (1) на множестве С1 имеет вид
^ = Л-01Ж + ^02х2 + Л-озху,
= кпх + Л,12у + Л-13 х2 + Л-14 ху + к^у2
(4)
где
йо1 = (1 - Ао)(1 - 27о)«о - Ао,
^02 = -воо(1 - Ао)(1 - 27о),
(5)
коз = -во1 (1 - Ао)(1 - 27о), к- = 2ао(1 - ЛЬ,
^12 = М(1 - А1)«1 - А1), к1з = -2воо(1 - Ао)7о,
к14 = -(^Ао(1 - А1) + 2во1 (1 - АоЫ, к15 = -^ви(1 - А1). Система (4) записана в модельном времени. Переход к модельному времени задавался
, Т
соотношением г = -о, где то — средняя продолжительность клеточного цикла в популяции
То
нормальных клеток, а т — реальное время. Переменными модели являются нормированные величины х = Х/Хо, у = У/Хо, где Хо численности нормальных клеток в начальный момент времени.
2. Постановка задачи параметрической идентификации математической модели
Перейдем в системе (4) к реальному времени т и ненормированным численностям, т.е. выполним обратные замены. После соответствующих преобразований система (4) примет вид
Г ЛХ -г
¿о! X + ко2Х + козх
— = к о1 X + ко2Х2 + козХУ, ат
ОТ
Лт
= ких + к 12У + к 1з х2 + кмху + к^г2,
где
к
г _ -о1 г
кот = —к|
то
ко1
то-, з =1;
з = 2 3
тоХо , 3 = 2, 3;
к
^ ■ 3 = 1.2;
У
У
т оХо
3 = 3, 4, 5;
параметры ку определены в (5).
(6)
Для решения задачи параметрической идентификации математической модели (6) необходимо получить ее дискретный аналог. Обозначим через —к, Ук численности нормальных и аномальных клеток в момент времени ¿к, а через —к+ь У+ — численность этих клеток в момент времени ¿к+ь Будем предполагать, что измерения численности нормальных и аномальных клеток проводятся через равные промежутки времени Д*, поскольку в результате эксперимента можно оценить количество клеток только в фиксированнные моменты времени [7, 8, 5]. Запишем разностные уравнения, соответствующие уравнениям (6):
-+ДД- = Л01—к + Л02—к + Л03—к Ук,
(7)
— у. ~ ~ ~ „ „
-Д- = Л11 — к + Л12Ук + Л13—к + Л14— кУк + Л15Ук .
В момент времени ¿к+1 численности нормальных и аномальных клеток задаются следующими соотношениями:
—к+1 = (ЛиД* + 1)—к + '02 Д*—2 + '03 Д*—к Ук, Ук+1 = ЛИД*—к + (Л^Д* + 1)Ук + ^1зД*—2 + ЛмД*—к Ук + ^15 Д*Ук2
Запишем математическую модель (8) в виде:
—к+1 = Ло1—к + 2 + Л°3—кУЬ
Ук+1 = Л11 —к + Л12Ук + Л13— + Л14—кУк + ^15УА2,
(8)
(9)
где
Л01 = '01Д* + 1, '12 = '12 Д* + 1,
ло Л02 = Л^Д*, Л13 = Л13 Д*,
ло Л03 = '03 Д*, Л14 = Л14 Д*,
Л11 = ЛПД*, Ло Л15 = Л15 д*.
(10)
Параметрам Л- можно дать следующую содержательную интерпретацию. Параметры определяют динамику развития популяции нормальных клеток, а параметры — динамику развития популяции аномальных клеток. Рост численности популяций нормальных и аномальных клеток характеризуют параметры Л°1, ЛЦ и Л°2, параметры Л°1 и Л12 определяют рост численности популяций за счет процессов митоза, параметр ЛЦ характеризует рост числа аномальных клеток за счет перерождения нормальных клеток в аномальные. Остальные параметры характеризуют уменьшение численности популяций за счет внутрипопуляцион-ной и межпопуляционной конкуренции.
Обозначим через — и У наблюдаемые значения компонент вектора состояния системы, а через —к, Ук — сами эти компоненты.
В итоге получим соотношения
"-01 "11
"12
"02 "13
"03 "14
"15
П \ /
Относительно оцениваемых параметров "1,. .., "(5 полученная модель является линейной. Для решения задачи идентификации и одновременного получения законов распределения параметров математической модели воспользуемся методом, основанным на использовании байесовского подхода [21] и теоремы об эквивалентности [22], предложенным в [18], [23]. Решение данной задачи удобно разделить на три этапа.
На первом этапе найдем точечные оценки "(, параметров модели (11).
На втором этапе проверим статистическую гипотезу о виде закона распределения случайных возмущений исходной модели состояния.
В случае, если проверка статистической гипотезы о виде закона распределения случайных возмущений показала, что гипотеза о нормальном распределении случайных возмущений должна быть отклонена, решить задачу параметрической идентификации математической модели (11) предложенными методами не удается. В противном случае, если гипотеза о нормальном распределении случайных возмущений принята, можно приступать к третьему этапу.
На третьем этапе решения задачи идентификации параметров находим законы распределения параметров нелинейной модели.
3. Нахождение точечных оценок параметров нелинейной модели
Воспользуемся методом, предложенным в [22], для нахождения точечных оценок элементов матрицы параметров модели (11). Математическую модель представим в виде
01
"02 "
03
1,1 1,(1 1,(1 1,(1 "11 "12 "13 "14 "15
/
Х22-1
\
К
22-1 Х 2
V 2 \ Г22-1
+
(12)
где матрицы-строки
£ = (£1,
£м) е М^м(Я), п = (П1, ..., Пм) е М^м(Я)
(13)
являются реализациями случайных возмущений (случайных ошибок наблюдения), обусловленных случайными возмущениями самой математической модели (11) и погрешностями в экспериментальных данных.
0
0
£
п
Задачу идентификации векторов параметров модели (12) можно разбить на две независимые подзадачи, как это было сделано в работе [23]. На первом этапе найдем оценки параметров , определяющих динамику развития популяции нормальных клеток, на втором — оценки параметров , определяющих динамику развития популяции аномальных клеток. Это возможно, поскольку обратная связь между реализацией случайного возмущения е в первом из уравнений состояния в математической модели (12) и реализацией случайного возмущения п во втором из уравнений состояния отсутствует и, следовательно, их ковариационная матрица имеет диагональный вид.
Построим п-й двухточечный шаблон первого уравнения состояния в (12) {¿2п-1; *2п}. Пусть —2п-1, —2п — экспериментальные значения численностей популяций нормальных клеток в моменты фиксации ¿2п-1 и *2п соответственно, где ¿2п-1 = (2п — 1)Д*, *2п = 2пД*.
Сформируем массивы экспериментальных данных для нахождение оценки параметров, определяющих динамику развития популяции нормальных клеток, используя указанный шаблон:
—ч — (—2, —4, . . . , —2^) е М^(Я), —нч = (—1, —3, . . . , —2^-1) е М1ХМ(Я). (14) Согласно (12), (14)
—ч = Л01—нч + Л02—нЧ + Л03—нчУнч + е = + е (15)
где
/
Дг —
—1 —3
—12 —1У1
—32
—3 У3
Ло , Ло , Ло Л01, Л02, Л03
—2М -1 —2
— 2М-1 —2М-1У2М— 1
—нч —н2ч
—У
у —нчУ нч у
(16)
(17)
матрица-строка е определена в (13).
Сформируем массивы экспериментальных данных для нахождение оценки параметров, определяющих динамику развития популяции аномальных клеток. Для этого построим п-й двухточечный шаблон второго уравнения состояния в (12) ({¿2п-1; ¿2п}), в узлах которого известны экспериментальные значения У2п-1, У2п в моменты времени ¿2п-1 и *2п соответственно, где ¿2п-1 = (2п — 1)Д*, *2п = 2пД*. Тогда
Уч — (У2, У4, ..., У^) е М^*(Я), Унч — (У1, У3, ..., У2^-1) е М^*(Я). (18) Из (12), (14) и (18) получим
У — Л11—нч + Л12Унч + Л13—нч + Л14—нчУнч + Л15Унч + П = Н1 + п
(19)
Здесь
Н? = (Л?!, Л?2, Ь?з, Л?4, Л
—
V
Х1 Хз . . Х2М -1 \ Х Хнч
К1 Кз . . К2М -1 К нч
X2 Хз2 . . . Х2 = Хн2ч
Х1К1 Х3У3 . . Х2М-1К2М- 1 ХнчКНч
К2 К2 . Кз . . . К2 / К2
(20)
(21)
матрица-строка п определена в (13).
Решая матричные уравнения (15) и (19), получаем оценки вектора параметров первого Н? (16) и второго Н? (20) уравнений модели:
Н?
Н?
(22) (23)
где — матрица, псевдообратная по отношению к матрице определенной соотношением (17); — матрица, псевдообратная по отношению к матрице , определенной соотношением (21).
4. Проверка статистической гипотезы о виде закона распределения случайных возмущений математической модели
На втором этапе решения задачи идентификации параметров модели (11) необходимо проверить статистическую гипотезу нормальном распределении случайный возмущений модели (невязок). Запишем невязки в виде
£ = £2,..., ] — Хч - Н?Дв, п = [п1, П2,..., Пм] — К, -
(24)
где матрицы экспериментальны« даннык Хч, Кч и определены равенствами (14), (17), (18) и (21) соответственно, Н?, Н? — оценки вектора параметров первого и второго уравнений модели. В качестве критерия согласия будем использовать критерий Колмогорова адаптированный Тюриным [24, 25] к выборкам малого обыема. Алгоритм проверки статистической гипотезы о нормальном распределении невязок хорошо изложен в [18, 23]. Была выбрана гипотезы: основная — невязки являются реализациями нормальных случайных величин, и стандартная альтернативная — невязки не представляют собой реализации нор-мальнык случайный величин. Уровень значимости а был выбран равным 0, 05, пороговое значение квантиля — Аа = 0,895.
Если нет оснований для отклонения основной гипотезы, то, как было показано в работах [21, 26], байесовские апостериорные законы распределения параметров Н? и Н?, построенные с учетом теории инвариантности Джеффриса, являются обобщенными распределениями Стьюдента.
5. Построение функций плотности распределения вероятностей
Третий этап решения задачи параметрической идентификации математической модели (11) состоит в построении функций плотности распределения вероятностей векторов параметров и указанной математической модели. Параметры соответствующих законов распределения необходимо знать для построения интервальные оценки для элементов векторов параметров математической модели. Воспользуемся результатами полученными в [23] для построения функций плотности распределения вероятностей параметров рассматриваемой модели.
Построение функции плотности распределения вероятностей вектора параметров модели Н^. Характер динамики численности популяции нормальных клеток в модели (11) определяется выражением (15) и имеет вид
Хч = + е
где матрица-строка е определена в (13), матрица определена в (17), матрица-строка Но = (^01, Л-°2, определяет характер динамики суммарной численности популяции нормальных клеток.
Считаем, что все возмущения ек, к = 1, ..., N, есть нормально распределенные случайные величины. Для величин ек, к = 1, ..., N, математические ожидания М(ек) равны нулю, дисперсии Д(ек) равны а2, случайные величины е^ и е^ не коррелированы. Заметим, что зависимые параметры жЧ, учитывая допущения о нестохастической природе жНЧ, также являются нормально распределенными случайными величинами.
В первом уравнении модели (12) воздействие неучтенных случайных факторов и ошибок определяется с помощью дисперсии случайных возмущений а2. Несмещенной оценкой дисперсии случайных возмущений а2 является выборочная остаточная дисперсия Я2
2 _ (Хч - Дс)(Хч - Д^хГ _ еет _ Е^=1 е!
Я2 =
N - 1 N - 1 N - 1
где е = Хч — — выбочная оценка реализации случайного возмущения ек.
0
0
Для оценки вектора параметров модели Д применим метод наименьших квадратов. Обозначим величину, равную сумме квадратов отклонений расчетных значений состояния системы от их экспериментальных значений Хч, через Я2. Тогда
'0
N
Я2 = = £ ек. (25)
к=1
Величина Яд пропорциональна выборочной дисперсии Я2.
Для нахождения оценки вектора параметров методом наименьших квадратов минимизируем величину
Я2 = (Хч — )(Хч — . (26)
Заметим, что
\ч(Н0Дв) = хтдт (Но ) = н0дТхт
Хч(Н?Дж) = ХчДв (Н?) = Н?ДжХч, (27)
так как произведения матриц Н?ДТХТ и ХчДХ(Н?)т есть скалярные величины, поскольку Н? е М1 хз, Дв е Мзх^, Хчт е х 1 и, следовательно, , ХчДТ(Н?)т е М1Х1.
После раскрытия скобок в (26) с учетом (27) получим
£2 = ХчХч — 2Н?ДжХч + Н?(Н?) . (28)
Продифференцируем выражение (28) по Н? и запишем необходимое условие экстремума
— 2ДЖХЧТ + 2ДЛТ (Я?)т = 0.
С учетом равенства (ДТДХ)т = ДТДХ из последнего соотношения находим оценку вектора Н ? в явном виде:
Н? = ХчДТ (ДтДТ )-1 = ХчД+. (29)
Функция £2, определенная равенством (28), является сильно выпуклой, поскольку матрица является положительно определеннной. Полученное необходимое условие минимума является и достаточным [29]. Найденное значение Н? соответствует минимуму функции £02(Н?). Отметим, что полученное представление (29) совпадает с (22) [27, 28].
Согласно формуле Байеса [21], совместная апостериорная функция плотности распределения вероятностей вектора параметров Н? и дисперсии а2 может быть представлена в виде
/(Н?, а2|Хч,Дт) ~ /(Н?,а2)/(Хт|Дх, Н?, а2), (30)
где / (ХТ | ДТ, Н?, а2) — условная плотность распределения вероятностей новык наблюдений ХТ при определенный значениях вектора параметров Н?, а2 и ДТ, /(Н?, а2) — совместная априорная функция плотности распределения вероятностей для вектора параметров Н? и а2. Функцию /(ХТ|ДТ, Н?, а2) будем рассматривать как функцию правдоподобия для параметров Н? и а2, обозначим ее через /(Н?, а2|ХТ, ДТ).
г? „ ^2
При нормальном распределении ошибок е и заданный значениях ДТ, Н? и а2 апостери-ную функцию плотности распределения вероятностей /(ХТ|ДТ, Н?, а2) можно записать в виде
, „ „ Г 1 , , ,л
(31)
/(Хч|Дт,Н?,а2) - (а2)-М/°ехр
— (Хч — ДтНО)(Хт — ДтН?)
2а2
Выражение (ХТ — ДТН?)(ХТ — ДТН?)т с учетом (25) и (29) преобразуем к виду
(ХТ — ДТН0)(ХТ — ДТН0) = £ 2 + (Н0 — Н0)ДТДТ (Н0 — Н0 ) . (32)
Функция правдоподобия для параметров Н? и а2 задается как плотность распределения
ГТ ¿т
вероятностей совместного появления результатов выборки в виде хТ, ¿Х, к = 1, N,
N
/(Н?,а2|Хч,Дт) = П /(хТ|Н?,а2).
к=1
Получим функцию правдоподобия для параметров и а2, рассматривая правую часть соотношения (31) как функцию от и а2 при фиксированных Хч и Дт и учитывая выражение (32):
/(Д0°,а2|Хч,^) ~ (а2)-^ ехр
1 с<2 1 ( ттО 7 ут / ттО т
— 2а2 Я0 — 202 (Н0 — Н0 )ДтАт (Но — Но)
(33)
Совместную априорную функцию плотности распределения вероятностей f (Н°, а2) будем считать расплывчатой [26, 21]. Примем допущение, что элементы вектора и дисперсия
а2 независимо распределены, т.е.
f «а2) = f (НО) f (а2).
Положим константой априорную функцию плотности распределения вероятностей вектора [26], поскольку о значении элементов вектора параметров ничего сказать нельзя, и формально элементы вектора могут принимать любые значения от —то до то. При этом
f (Д0) = 0СП81.
(34)
f (а2) —.
(35)
Поскольку а2 есть положительная величина, априорная функция плотности распределения вероятностей а2 будет иметь вид [26]
1 а
С учетом (33), (34) и (35) соотношение (30), определяющее совместную апостериорною функцию плотности распределения вероятностей неизвестного вектора параметров и выборочной дисперсии а2, примет вид
N+ 1
f (Д0>2|Хч,^) - (а2)-^ ехр
— 2а2 (Я2 + (Н° — Но)^ ДТ(Но — #0)т)
(36)
Для упрощения записи положим х = а 2 и перепишем выражение (36) следующим образом:
N+1
f (Я0°,х|Хч,Дж) - ехр
—X (Я2 + (До — яО)дХ (Н — НО)т)х
(37)
Проинтегрируем (37) по х для получения апостериорной функции плотности распределения вероятностей вектора параметров [21]. В результате получим
f (До°|Хч, Дх) - (Яо2 + (Но0 — Н0)Дх(Но0 — НО)т)-т X
ехр
х
—X (Я02 + (НО — Я0)ДДТ (НО — Н0)т)
N+1
X 2
Я2 + (НО — яО)ДтДт(яО — НО)т
т " 2
-¿X. (38)
Можно показать, что интеграл в выражении (38) сходится и его значение равно константе, поэтому плотность совместного распределения вероятностей вектора параметров модели Н может быть представлена так:
f (НО|Хч, Дт) - я2 + (НО — яО)(ДтДтт)(НО — я0О)
О
т
т " 2
(39)
Здесь
£2 — (ХТ — Н0ДТ)(ХТ — Н0ДТ) = ее ,
матрицы экспериментальны« данных ХТ и определены равенствами (14) и (17), оценка Н? вектора параметров Н? — равенством (22).
Можно виде ь, ч о полученный апос ериорный закон распределения век ора параметров Н?,построенный на основании байесовского подхода с учетом теории инвариантности Джеффриса [21, 26], является обобщенным распределением Стьюдента.
Построение функции плотности распределения вероятностей вектора параметров модели НВ математической модели (11) характер динамики суммарной численности популяции аномальных клеток полностью определяется выражением (15):
Уч = Н? Ду + п,
где матрица-строка п определена в (13); матрица определена в (21); матрица-строка Н? = (к^, Л,?2, Л?, , определяет характер динамики суммарной численности популяции аномальных кле ок.
Как и при построении функции плотности распределения вероятностей параметров Н?, считаем все ошибки наблюдения пк, к = 1, ..., N, нормально распределенными случайными величинами с математическим ожиданием, равным нулю, и дисперсией, равной а2. Случайные величины п» и п не коррелированы.
По аналогии с методикой построения функции плотности распределения вероятностей параметров Н?, несмещенной оценкой дисперсии случайных возмущений а2 является оста-очная выборочная дисперсия
= (Уч — Н?Ду)(Уч — )т = = " ,2
N — 1 N — 1 N — 1 ¿1Пк,
где п = УТ — Н?Ду — выбочная оценка реализации случайного возмущения пк.
Обозначим через величину, пропорциональную выборочной дисперсии возмуще-
2
ний а2:
N
£2 = Е п2. (40)
к=1
Применяя метод наименьших квадратов для оценки вектора Н?, получим
Я? = УтДут (Ду )-1 = УчД+. (41)
Представим совместную апостериорную функцию плотности распределения вероятностей вектора параметров Н? и дисперсии а2 в виде
/(Н?,а2|Ут,Ду) - /(Н?,а2) /(У,, Н?, X2), (42)
где f , Н, а2) —условная плотность распределения вероятностей новых наблюдений
Уч при определенный значениях вектора параметров Н, а2 и ; f (Н, а2) — совместная априорная функция плотности распределения вероятностей для вектора параметров Н и а2.
Запишем апостериорную функцию плотности распределения вероятностей при нормальном распределении ошибок п и заданный значениях , Н и а2
f (Y4|Zy,Hf,a2) - (а2)-ж exp
1
■d\т
- (Уч - ZyH?)(Y4 - ZyHf)
(43)
Далее, по аналогии с методикой построения функции плотности распределения вероятностей параметров Hf, рассмотрим f(Y4|Zy, Hf,a2) как функцию правдоподобия и положим, что априорные функции плотности распределения вероятностей вектора Hf постоянны (f (Hf) = const), а дисперсии а2 имеют вид f (а2) — ^). Тогда плотносты совместного распределения вероятностей вектора параметров модели Hf можно представиты следующим образом:
тИ\т'
f (Hf|Y4, Zy) - S2 + (Hf - Hf)(ZyZ/)(Hf - Hf)
ж " 2
(44)
Здесь = (Уч — )(Уч — )т = ппт, матрицы экспериментальных даннык Хч, Zz определены равенствами (18) и (21), а оценка Н вектора параметров Н — равенством (23).
Апостериорный закон распределения вектора параметров Н, также как и (39), является обобщенным распределением Стьюдента.
Заключение
Разработанная методика получения точечных оценок параметров математической модели, описывающей динамику селективного размножения клонообразующей популяции аномальный: клеток в культуре стволовык клеток человека при стандартных лабораторных условиях культивирования с учетом контактного торможения. На основе байесовского подхода получены функции плотности распределения вероятностей параметров математической модели. Знание функций плотности распределения вероятностей параметров математической модели необходимо для нахождения их маргинальны« функций плотности распределения и построения на их основе доверительных интервалов для параметров модели. Доверительные интервалы могут служить основой для оценки вероятности реализации различных сценариев развития клеточнык популяций.
Оценка вероятности реализации различнык сценариев развития клеточнык популяций может быть проведена, например, методами математического моделирования, основываясь на идеях предложенных в работе [18].
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 13-07-00720.
Список литературы
1. Kresnowati M.T., Forde G.M., Chen X.D. Model-based analysis and optimization of bioreactor for hematopoietic stem cell cultivation // Bioprocess and Biosystems Engineering. 2011. Vol. 34, no. 1. P. 81-93. DOI: 10.1007/s00449-010-0449-z
2. Winkler D.A., Burden F.R. Robust, quantitative tools for modelling ex-vivo expansion of haematopoietic stem cells and progenitors//Molecular BioSystems. 2012. Vol. 8, no. 3. P. 913920. DOI: 10.1039/c2mb05439f
3. Ducrot A., Le foll F., Magal P., Murakawa H., Pasquier J., Webb G.F. An in vitro cell population dynamics model incorporating cell size, quiescence, and contact inhibition // Mathematical Models and Methods in Applied Sciences. 2011. Vol.21, iss. supp01. P. 871-892. DOI: 10.1142/S0218202511005404
4. БочковН.П., НикитинаВ.А., Буяновская O.A., ВоронинаЕ.С., ГольдштейнД.В., Кулешов Н.П., Ржанинова А.А., Чаушев И.Н. Анеуплоидия в стволовых клетках, выделенных из жировых тканей человека // Бюллетень экспериментальной биологии и медицины. 2008. Т. 146, №9. C. 320-323.
5. Бочков Н.П., Никитина В.А., Воронина Е.С., Кулешов Н.П. Методическое пособие по тестированию клеточных трансплантатов на генетическую безопасность // Клеточные технологии в биологии и медицине. 2009. №4. C. 183-189.
6. Бочков Н.П., Виноградова М.С., Волков И.К., Воронина Е.С., Кулешов Н.П. Статистический анализ клонообразования в культурах стволовых клеток человека // Клеточные технологии в биологии и медицине. 2011. №2. C. 63-66.
7. Бочков Н.П., Никитина В.А. Цитогенетика стволовых клеток человека // Молекулярная медицина. 2008. №3. С. 40-47.
8. Бочков Н.П., Никитина В.А., Рослова Т.А., Чаушев И.Н., Якушина И.И. Клеточная терапия наследственных болезней // Вестник РАМН. 2008. № 10. С. 20-28.
9. Осипова Е.Ю., Шаманская Т.В., Пурбуева Б.Б., Устюгов А.Ю., Астрелина Т.А., Яковлева М.В., Румянцев С.А. Культивирование мезенхимальных стволовых клеток ex vivo в различных питательных средах (обзор литературы и собственный опыт) // Онкогематология. 2010. №3. C. 65-71.
10. Duesberg P., Mandrioli D., McCormack A., Nicholson J.M. Is carcinogenesis a form of speciation? // Cell Cycle. 2011. No. 10, P. 2100-2114.
11. Duesberg P., Li R., Fabarius A., Hehlmann R. Aneuploidy and Cancer: From Correlation to Causation // Infection and Inflammation: Impacts on Oncogenesis / Dittmar T., Zaenker K.S., Schmidt A., eds. Basel: Karger. 2006. P. 16-44 (ser. Contributions to Microbiology, vol. 13). DOI: 10.1159/000092963
12. Duesberg P., Fabarius A., Hehlmann R. Aneuploidy, the Primary Cause of the Multilateral Genomic Instability of Neoplastic and Preneoplastic Cells // IUBMB Life. 2004. Vol. 56, no. 2. P. 65-81. DOI: 10.1080/15216540410001667902
13. Тимошевский B.A., Назаренко C.A. Биологическая индикация мутагенных воздействий и генетической нестабилыности у человека путем учета числовых хромосомных нарушений//Информационный вестник ВОГиС. 2006. Т. 10, № 3. C. 530-539.
14. Бочков Н.П., Виноградова М.С., Волков И.К., Кулешов Н.П. Математическая моделы суммарных численостей взаимодействующих клеточных популяций // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2011, № 1. С. 18-24.
15. Виноградова М.С. Качественный анализ модели функционирования взаимодействующих клеточных популяций // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2011. № 11. С. 1-20. Режим доступа: http://technomag.edu.ru/doc/251409.html (дата обращения: 20.08.2015).
16. Виноградова М.С. Динамическая моделы клеточной популяционной системы // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2013. №12. С. 175-192. DOI: 10.7463/1213.0646463
17. Виноградова М.С. Анализ сценариев развития клеточной популяционной системы // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2014. № 11. С. 607-622. DOI: 10.7463/1114.0735732
18. Бочков Н.П., Виноградова М.С., Волков И.К. Оценка вероятности реализации вариантов развития взаимодействующих клеточных популяций // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2011, № 3. С. 31-43.
19. Abercrombie M. Contact inhibition in tissue culture // In Vitro. 1970. Vol. 6, no. 2. P. 128-142. DOI: 10.1007/BF02616114
20. Виноградова М.С. Исследование нелинейной модели развития клеточной популяцион-ной системы // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2014. № 8. С. 123-138. DOI: 10.7463/0814.0720269
21. Зелынер А. Байесовские методы в эконометрии. М.: Статистика, 1980. 440 с.
22. Волков И.К. Условия идентифицируемости математических моделей эволюционных процессов по дискретным косвенным измерениям вектора состояния // Известия РАН. Техническая кибернетика. 1994. №6. C. 55-72.
23. Виноградова М.С. Параметрическая идентификация модели взаимодействующих клеточных популяций на основе байесовского подхода // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журнал. 2012. № 11. C. 155-182. DOI: 10.7463/1112.0490900
24. Тюрин Ю.Н. Непаметрические методы статистики. М.: Знание, 1978. 64 с.
25. Тюрин Ю.Н. О пределыном распределении Колмогорова — Смирнова для сложной гипотезы//Известия АН СССР. Математика. 1984. Т. 48, №6. C. 1314-1343.
26. Jeffreys H. Theory of probability. Oxford: Clarendon Press, 1983. 459 с.
27. Алберт А. Регрессия, псевдоинверсия и рекурентное оценивание. М.: Наука, 1977. 224 с.
28. Гантмахер Ф.Р. Теория матриц. М.: Наука, 1988. 552 с.
29. Беллман Р. Введение в теорию матриц. М.: Наука, 1969. 368 c.
Science and Education of the Bauman MSTU,
Science S Education ™ ¡^P™
of the Bauman MSTU Received: 6.12.2015
ElectroTLIC joumdl ® Bauman Moscow State Technical University
Point Estimates and Probability Distribution of Mathematical Model Parameters of Evolving Cell Population System Taking into Account Contact Inhibition
Vinogradova M. S.*'* *mrnvinogradova@rambler.ru
1 Bauman Moscow State Technical University, Russia
Keywords: mathematical model, stem sell, parameter identification, Bayesian inference, point estimation, probability density function
The paper considers a mathematical model describing dynamics of the cell population system consisting of normal and abnormal stem cells. It offers a technique to provide the estimates and a distribution of the model parameters on the basis of limited sampling of experimental data. The model is a non-linear system of ordinary differential equations of the first order to determine the changing populations of normal and abnormal cells. The number of populations serves as the model variables. The model takes into account the limited resources and is nonlinear regarding the model variables, but it is linear with respect to unknown parameters and is considered in the domain where both populations evolve. To solve the problem of parametric identification of mathematical model was obtained its discrete analogue. The task of identifying vectors of the model parameters was divided into two independent sub-tasks. The first stage estimated parameters, which determined dynamics of the normal cell population evolution. The second one estimated parameters defining dynamics of the abnormal cell population evolution. The features of the proposed technique are to select a family of disjoint templates that allows us to have both the point estimates and the posterior distributions of the mathematical model parameters for small samples. Relationships are obtained to allow us to calculate point estimates and construct a probability density function of the vectors of the model parameters using the limited samples of experimental data. When building a distribution of the model parameters a Bayesian approach and a Jeffrey's theory of invariance were used. The obtained distributions of parameters are the generalized t-distribution.
References
1. Kresnowati M.T., Forde G.M., Chen X.D. Model-based analysis and optimization of bioreactor for hematopoietic stem cell cultivation. Bioprocess andBiosystems Engineering, 2011, vol. 34, no. 1, pp. 81-93. DOI: 10.1007/s00449-010-0449-z
2. Winkler D.A., Burden FR. Robust, quantitative tools for modelling ex-vivo expansion of haematopoietic stem cells and progenitors. Molecular BioSystems, 2012, vol. 8, no. 3, pp. 913— 920. DOI: 10.1039/c2mb05439f
3. Ducrot A., Le foll F., Magal P., Murakawa H., Pasquier J., Webb G.F. An in vitro cell population dynamics model incorporating cell size, quiescence, and contact inhibition. Mathematical Models and Methods in Applied Sciences, 2011, vol.21, iss. supp01, pp. 871—892. DOI: 10.1142/S0218202511005404
4. BochkovN.P., Nikitina V.A., Buyanovskaya O.A., VoroninaE.S., Gol'dshteynD.V., Kuleshov N.P., Rzhaninova A.A., Chaushev I.N. Aneuploidy of stem cells isolated from human adipose tissue. Byulleten' eksperimental'noy biologii i meditsiny, 2008, vol. 146, no. 3, pp. 320—323. (English translation: Bulletin of Experimental Biology and Medicine, 2008, vol. 146, iss. 3, pp. 344—347. DOI: 10.1007/s10517-008-0293-1).
5. Bochkov N.P., Nikitina V.A., Voronina E.S., Kuleshov N.P. Methodological Guidelines for Genetic Safety Testing of Cell Transplants. Kletochnye tekhnologii v biologii i meditsine, 2009, no. 4, pp. 183—189. (English translation: Bulletin of Experimental Biology and Medicine, 2009, vol. 148, iss. 4, pp. 677—683. DOI: 10.1007/s10517-010-0793-7).
6. Bochkov N.P., Vinogradova M.S., Volkov I.K., Voronina E.S., Kuleshov N.P. Statistical Analysis of Clone Formation in Cultures of Human Stem Cells. Kletochnye tekhnologii v biologii I meditsine, 2011, no. 2, pp. 63—66. (English translation: Bulletin of Experimental Biology and Medicine, 2011, vol. 151, iss. 4, pp. 498—501. DOI: 10.1007/s10517-011-1366-0).
7. Bochkov N.P., Nikitina V.A. Cytogenetics of human stem cells. Molekulyarnaya meditsina = Molecular medicine, 2008, no. 3, pp. 40—47. (in Russian)
8. Bochkov N.P., Nikitina V. A., Roslova T.A., Chaushev I.N., Yakushina I.I. Cellular therapy of hereditary diseases. Vestnik RAMN = Annals of the Russian Academy of Medical Sciences, 2008, no. 10, pp. 20—28. (in Russian)
9. Osipova E.Yu., Shamanskaya T.V., Purbueva B.B., Ustyugov A.Yu., Astrelina T.A., Yakovleva M.V., Rumyantsev S.A. Ex vivo expansion of mesenchymal stem cells in different culture conditions (the literature review and own experience). Onkogematologiya = Oncohematology, 2010, no. 3, pp. 65—71. (in Russian)
10. Duesberg P., Mandrioli D., McCormack A., Nicholson J.M. Is carcinogenesis a form of speciation? Cell Cycle, 2011, no. 10, pp. 2100—2114.
11. Duesberg P., Li R., Fabarius A., Hehlmann R. Aneuploidy and Cancer: From Correlation to Causation. In: Dittmar T., Zaenker K.S., Schmidt A., eds. Infection and Inflammation: Impacts on Oncogenesis. Basel, Karger, 2006, pp. 16—44. (Ser. Contributions to Microbiology; vol. 13.). DOI: 10.1159/000092963
12. Duesberg P., Fabarius A., Hehlmann R. Aneuploidy, the Primary Cause of the Multilateral Genomic Instability of Neoplastic and Preneoplastic Cells. IUBMB Life, 2004, vol. 56, no. 2, pp. 65—81. DOI: 10.1080/15216540410001667902
13. Timoshevsky V.A., Nazarenko S.A. Biological Indication of the Mutagenic Influences and Genetic Instability in Human Using Evaluation of Numerical Chromosome Aberrations. In-formatsionny Vestnik VOGiS = VOGiSHerald, 2006, vol. 10, no. 3, pp. 530—539 (in Russian).
14. Bochkov N.P., Vinogradova M.S., Volkov I.K., Kuleshov N.P. Mathematical Model of Dynamics of Total Quantities of Interacting Cell's Populations. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki = Herald of the Bauman MSTU. Ser. Natural science, 2011, no. 1, pp. 18—24. (in Russian)
15. Vinogradova M.S. The qualitative analysis of model of functioning cooperating cellular populations. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2011, no. 11, pp. 1—20. Available at: http://technomag.edu.ru/doc/251409.html, accessed 07.07.2014. (in Russian)
16. Vinogradova M.S. A dynamic model of the cellular population system. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2013, no. 12, pp. 175—192. DOI: 10.7463/1213.0646463 (in Russian)
17. Vinogradova M.S. Analysing Scenarios of Cell Population System Development. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2014, no. 11. pp. 607—622. DOI: 10.7463/1114.0735732 (in Russian)
18. Bochkov N.P., Vinogradova M.S., Volkov I.K. Estimating the Probability of Implementation of Variants of Development of Interacting Cell's Populations]. Vestnik MGTU im.N.E. Baumana. Ser. Estestvennye nauki = Herald of the Bauman MSTU. Ser. Natural science, 2011, no. 3, pp. 31—43. (in Russian)
19. Abercrombie M. Contact inhibition in tissue culture. In Vitro, 1970, vol. 6, no. 2, pp. 128—142. DOI: 10.1007/BF02616114
20. Vinogradova M.S. Investigation of the Nonlinear Model of the Cellular Population System Development. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2014, no. 8, pp. 123—138. DOI: 10.7463/0814.0720269 (in Russian)
21. Zel'ner A. Bajesovskie metody v ekonometrii [Bayesian methods in econometrics]. Moscow, Statistika publ., 1980. 440 p. (in Russian)
22. Volkov I.K. Identifiability conditions of mathematical models of evolutionary processes over discrete indirect measurements of the state vector. Izvestija Rossiiskoi akademii nauk. Tehnich-eskajakibernetika = Bulletin of the Russian Academy of Sciences. Technical Cybernetics, 1994, no. 6. pp. 55—72. (in Russian)
23. Vinogradova M.S. Parametrical identification of a model of cooperating cellular populations on the basis of Bayesian approach. Nauka i obrazovanie MGTU im. N.E. Baumana = Science
and Education of the Bauman MSTU, 2012, no. 11, pp. 155-182. DOI: 10.7463/1112.0490900 (in Russian)
24. Tyurin Yu.N. Nepametricheskie metody statistiki [Nonparametric methods of statistics]. Moscow, Znanie publ., 1978. 64 p.
25. Tyurin Yu.N. ON the limit distribution of Kolmogorov-Smirnov statistics for a composite hypothesis Izvestija AN SSSR. Matematika, 1984, vol.48, no. 6. pp. 1314-1343. (English Translation: Mathematics of the USSR — Izvestiya, 1985, vol. 25, no. 3, pp. 619-646 DOI: 10.1070/IM1985v025n03ABEH001311).
26. Jeffreys H. Theory of probability. Oxford, Clarendon Press, 1983. 459 p.
27. Albert A. Regression, And The Moor-Penrose Pseudoinverse. New York and London, Academic Press, 1972. 180 p. (Mathematics in Science and Engineering, vol.94). (Russ. ed.: —authAlbert A. Regressiia, psevdoinversiia i rekurentnoe otsenivanie. Moscow, Nauka publ., 1977. 224 p.).
28. Gantmacher F.R. Teorija matrits [The theory of matrices]. Moscow, Nauka publ., 1988. 552 p. (in Russian)
29. Bellman R. Introduction to Matrix Analysis. 2nd ed. Philadelphia, SIAM, 1997. 428 p. (Ser. Classics in applied mathematics; vol. 19). (Russ. ed.: Vvedenie v teoriju matrits. Moscow, Nauka publ., 1969. 368 p.)
