Оценка вероятности реализации вариантов развития взаимодействующих клеточных популяций Текст научной статьи по специальности «Математика»

УДК 62.505

Н. П. Бочков, М. С. Виноградова, И. К. Волков

ОЦЕНКА ВЕРОЯТНОСТИ РЕАЛИЗАЦИИ ВАРИАНТОВ РАЗВИТИЯ ВЗАИМОДЕЙСТВУЮЩИХ КЛЕТОЧНЫХ ПОПУЛЯЦИЙ

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

E-mail: mathmod@bmstu.ru

Ключевые слова: математическая модель, вероятность реализации, динамика численности, стволовые клетки, селективное размножение.

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

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

В [2] также предложен возможный путь практического использования полученной модели, проверка реализуемости которого методами имитационного моделирования [3] и является основной целью настоящих исследований.

Возможные варианты динамик численностей популяций нормальных и аномальных клеток. Согласно [2] при стандартных условиях культивирования в лабораторных условиях суммарные численности популяций X (t) нормальных и Y (t) аномальных клеток в момент

фиксации £ > 1 определяются как решения следующей задачи Коши:

X((п0 + 1)т0) = о°Х(п0т0), п0 е {1, 2,...}; У ((п1 + 1)т^ = а1^ (п1т1) + (п1тп1 е {1, 2,...}; (1)

X(0) = Х0; У(0) = 0.

Комплексные параметры

а0 = (1 - Мг0)(1 - АЕ0) + 2Мг0(1 - АЕ0)(1 - 70),

а1 = (1 - Мг1)(1 - АЕ1) + 2Мг1(1 - АЕ1), (2)

Ь1 = 2Мг0(1 - АЕ0)70

полностью определяют характер динамики сумммарной численности как популяции нормальных, так и популяции аномальных клеток, образующих изучаемую популяционную систему. В (1) и (2) использованы следующие обозначения: верхний индекс ] е {0,1} относится к характеристикам ]-й клеточной популяции, где значение ] =0 соответствует популяции нормальных, а ] = 1 — аномальных клеток; т^ — средняя продолжительность клеточного цикла в ]-й популяции, а п — номер среднего клеточного цикла в ней; X(п0т0) — численность популяции нормальных клеток в момент фиксации £ = п0т0; У (п1т1) — численность популяции аномальных клеток в момент фиксации £ = п1т1; Мг-7 — доля клеток в ] -й популяции, разделившихся на временном интервале длительности т^; АЕ-7 — доля клеток в ^'-й популяции, погибающих на временном интервале продолжительностью т^; 70 — доля нормальных клеток в изучаемой популяционной системе, переходящих на временном интервале длительностью т0 в процессе деления в популяцию аномальных клеток.

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

dX(t) _ 1 {

dt = Т0 "

^ = 4(а0 - 1)X(t), t> 0;

^ = ^ ((а1 - 1)У(£) + Ь1Х(£)), £> 0; (3)

X(0)= Х0; У(0) = 0,

где £ — текущее время, единица которого соответствует средней продолжительности клеточного цикла популяции нормальных клеток.

Для удобства анализа выполним в (3) замену времени £ = т0£, а также введем новые переменные х = X/X0, у = Y/X0, где X0 — число

посаженных клеток. Тогда система (3) примет вид

^ = (а0 - 1)х, Ь > 0;

| = М(а* - 1)у + Ь1 х),Ь> 0; (4)

х(0) = 1; у(0) = 0,

где коэффициент д = т0/т1 характеризует отношение скоростей размножения аномальных и нормальных клеток.

Решение задачи Коши (4) находим стандартными методами [4]:

х(Ь) = е(а°-1)*, Ь > 0;

дЬ1

(e(a°-1)t - -1)t), ^ = 0, t > 0;

У (t) = { У

дЬ1 te(a°-1)t, p = 0, t > 0,

(5)

где ^ = а0 — 1 — ^(а1 — 1).

Проведем параметрический анализ полученных решений.

На рис. 1 приведены графики возможного изменения суммарной численности популяции нормальных клеток (переменной х) от времени Ь. Здесь кривая 1 соответствует случаю, когда а0 > 1, кривая 2 — случаю, когда а0 = 1, а кривая 3 — случаю а0 < 1.

На рис.2 приведены графики возможного изменения суммарной численности популяции аномальных клеток (переменной у) от времени Ь при ^ = 0. Здесь кривая 1 соответствует случаю, когда а0 > 1, кривая 2 случаю, когда а0 = 1, а кривая 3 — случаю а0 < 1. На кривой 3 имеем точку Ь* = 1/(1 — а0), у * = Ь1 д/(1 — а0 )е.

-

2

2.5

Рис. 1. Изменение суммарной численности популяции нормальных клеток

1.5

У

Рис. 2. Изменение суммарной численности популяции аномальных клеток при ф = 0

На рис. 3 приведены графики возможного изменения суммарной численности популяции аномальных клеток (переменной у) от времени £ при ^ = 0. Здесь кривая 1 соответствует случаю, когда значение хотя бы одного из параметров, а0 или а1, больше 1, т.е. а0 > а1 > 1 или а0 > 1 > а1, или а1 > а0 > 1, или а1 > 1 > а0.

Кривая 2 соответствует случаю, когда а0 = 1 > а1 или а1 = 1 > а0.

Ь1 0 Ь1

Предельное значение у * = --— при а0 = 1 и у * = --— при

(1 - а1 )е (1 - а0)е

а1 = 1.

На рис. 4 приведен график возможного изменения суммарной численности популяции аномальных клеток (переменной у) от времени

Рис. 3. Изменение суммарной численности популяции аномальных клеток при У = 0_

34 ISSN 1812-3368. Вестник МГТУ им. Н.Э. Баумана. Сер. "Естественные науки". 2011. № 3

Рис. 4. Изменение суммарной численности популяции аномальных клеток при а0 < 1 и а1 < 1, а0 = а1

t при а0 < 1 и а1 < 1, а0 = а1. При этом

„ L (а1 - 1)д t* = -lnv 7 •

У * =

^ а0 - 1 ' Ь1 д ((а1 - 1)д (а°-1)/^ (а1 - 1)д

^ \ а0 - 1 а0 - 1

Таким образом, характер динамики суммарных численностей популяций нормальных и аномальных клеток, образующих изучаемую популяционную систему, полностью определен значениями параметров а0 и а1, входящих в математическую модель (1) и определяемых равенствами (2), а значение параметра Ь1 определяет максимальные значения переменной у в тех случаях, когда ее изменение ограничено.

В случае, если д =1, можно выделить три области значений параметров а0 и а1, при которых реализуются различные сценарии. Первая область задается неравенствами 0 < а0 < 1 и 0 < а1 < 0, вторая область — неравенствами 0 < а0 < 1 и а1 > 1, а третья — неравенством

а0 > 1.

В первой области с течением времени численность обоих популяций стремится к нулю. Во второй области численность нормальных клеток стремится к нулю, а численность аномальных — возрастает. В третьей области растет численность обеих популяций.

Вернемся к рассмотрению дискретной модели (1) и положим т0 = т1. Перейдем к нормализованным переменным х = Х/Х0, у = У/Х0, где Х0 — число посаженных клеток. Тогда х(0) = 1,

y(0) = 0 и решение соответствуюшей задачи Коши примет вид

x(k) = (а0)';

b1((a0)k - (а1)') (а0

b1(a0)kk 0 i ——, а0 = а1, а0

У(к)= ,1, 0а0 - а1)

, а0 = а1; (6)

Сравнивая (5) и (6), можно видеть, что поведение решений задачи Коши для непрерывной и дискретной моделей похоже.

Решение задачи параметрической идентификации математической модели (1). Решение указанной задачи с использованием теорем об эквивалентности [5] и байесовского подхода [6] удобно представить как последовательную реализацию следующих трех этапов.

Этап 1. Нахождение точечных оценок а0, а1, Ь1 параметров а0, а1, Ь1 математической модели (1).

Пусть далее {(2п0 - 1)т0; 2пт0} — п-й двухточечный шаблон первого уравнения состояния в (1), в узлах которого известны экспериментальные значения X2n-1,X2n суммарной численности популяций нормальных клеток в моменты фиксации £ = (2п0 - 1)т0 и £ = 2п0т0 соответственно. Формируем массивы экспериментальных данных

Xч 4 ^^...,X2N] е М^(Л), ^ч = , , . . . , X2N-1 ] е MlxN(Л).

Согласно (1) и (7)

^ = Лнч + ^ (8)

где матрица-строка е = [е1, е2,... бN] е M1xN (Л) — реализация случайного возмущения первого из уравнений состояния в исходной математической модели (1), обусловленного как случайными возмущениями самой математической модели, так и погрешностями экспериментальных данных. Таким образом, воспользовавшись теорией псевдообратных матриц [7] и равенством (8), можно утверждать, что оценка МНК параметра а0 может быть представлена в явном виде

а0 = xчxн+ч, (9)

где XнЧ — матрица, псевдообратная по отношению к матрице Xнч, определенной в (7).

Пусть далее {(2га1 - 1)т0; 2пт1} — п-й двухточечный шаблон второго уравнения состояния в (1), в узлах которого известны экспериментальные значения У2п-1, У2п суммарной численности популяций аномальных клеток в моменты фиксации £ = (2га1 - 1)т1 и £ = 2п1т1

' Y1,Y3,.. . , Y2M-1 \Y 1 1 нч

. , X2M-1 Хнч.

соответственно. Формируем массивы экспериментальных данных

(Д),

Унч — , , . . . , У2М-1] е М1хМ(10)

е М2хм(Л)

и согласно (1), (10) получаем

Уч = [а1,61]^ + п, (11)

где матрица-строка

п = [п1, П2,..., пм] е ^хи(Л)

является реализацией случайного возмущения второго из уравнений состояния в математической модели (1). Таким образом, воспользовавшись теорией псевдообратных матриц [7] и равенством (11), найдем оценки МНК параметров а1, Ь1:

[а1,Ь1] = к %+. (12)

Этап 2. Проверка статистической гипотезы о виде закона распределения случайных возмущений исходной модели состояния.

После нахождения точечных оценок для неизвестных параметров исходной математической модели (1) и непосредственно перед построением для них интервальных оценок согласно общей теории и в соответствии с (8), (9), (11), (12), необходимо проверить статистическую гипотезу о том, что невязки

б = [б1, ¿2,..., бN ] — хч — а0Хнч,

п = [П1, П2,..., пм ] — кч — [а1, Ь1^,

где матрицы экспериментальных данных Хч, Хнч, Уч, % определены равенствами (7) и (10), представляют собой реализации нормальных случайных величин при наличии стандартной альтернативной гипотезы. Для решения этой задачи при проведении настоящих исследований использовали один из наиболее "жестких" критериев — критерий Колмогорова-Смирнова в модификации Ю.Н. Тюрина [8], адаптированной к выборкам малого объема, и стандартный алгоритм, проиллюстрированный ниже на примере матрицы строки невязок б, определенной в (13):

Шаг 1. Вычисление выборочного среднего

N

А v г*1 ^ (13)

- А

N

1

к-1

е = ^ (14)

и исправленной выборочной дисперсии

i N

52 = - Ч" (15)

к-1

матрицы-строки невязок е, определенной равенствами (13), (9), (7);

Шаг 2. Центрирование и нормирование элементов матрицы-строки е, сортировка чисел

N

V}« ('6)

в порядке возрастания.

Шаг 3. Определение по таблицам или с использованием вычислительных алгоритмов значения {Е(уп)}^=1 функции распределения стандартного нормального закона, где через уп, п = , обозначен п-й элемент упорядоченного множества, полученного на шаге 2; Шаг 4. Вычисление значения

Л = sup F(yn) - NN

N - o,oi + °,86

N

Шаг 5. Применение решающего правила. При уровне значимости а = 0,05 и пороговом значении квантиля Аа = 0,895 воспользоваться следующим решающим правилом:

А < Аа (с уровнем значимости а нет оснований для отклонения основной гипотезы — приступить к реализации этапа 3);

А > Аа (отклонить основную гипотезу в пользу конкурирующей и вернуться к анализу экспериментальных данных и исходной математической модели для установления их соответствия изучаемому процессу).

Этап 3. Апостериорные законы распределения параметров а0, а1, Ь1 и их интервальное оценивание.

Если с заданным уровнем значимости нет оснований для отклонения основной гипотезы о том, что матрицы-строки невязок е, п, определенные равенствами (13), являются матрицами реализаций нормальных случайных величин, то можно доказать [5], что байесовские апостериорные законы распределения параметров а0, а1, б1, построенные с учетом теории инвариантности Джеффриса [6], являются распределением Стьюдента. Их плотности распределения вероятности могут быть представлены в следующем виде:

N

/(а0^^) - {52 + ^н^н/)(а0 - а0)2}--2, (17) 52 = (Хч - а^Х^ - а0^)Т = ббт;

м

/(а1,61|Уч,^) - {£2 + [а1 - а1,^1 - б1]^[а1 - а1 ,Ь1 - Ь1]т}-М,

52 = (уч - [а1,^)(уч - [а1 )т = .

(18)

При этом в (17) матрицы экспериментальных данных Xч,Xнч определены равенствами (7), а оценка а0 параметра а0 — равенством (9). В (18) матрицы экспериментальных данных У, ^ определены равенствами (10), а оценка [а1, Ь1] параметров а1, Ь1 — равенством (12). Согласно (17) число степеней свободы

= N - 1 (19)

и с используемым уровнем значимости а полуразмах симметричного доверительного интервала (а0 - /(а0) а0 + /(а0)) для параметра а0 определяется равенством [6]

/ £2

/(а0) = Да0^^нч^а) = Л N(^^г) М-а/гК),

где определена равенством (19), а ^1_а/2(^а) — квантиль распределения Стьюдента с числом степеней свободы ра уровня 1 - а/2, где определена равенством (19), а а = 0,05.

Число степеней свободы в двумерном распределении Стьюдента (18)

V = М - 2 (20)

и с уровнем значимости а полуразмахи /(а1) и /(б1) симметричных доверительных интервалов для параметров а1 и Ь1 определяются соответственно равенствами [6]

/(а1) = /(а1, Уч, V, а) = ^М-1£2^и );

/(Ь1) = /(Ь1, Уч, Z, V, а) = уМ_1£я22 ^1_а/2^),

где Zkk — к-й элемент матрицы (ZZT)-1, а число степеней свободы определяется равенством (20).

Определение вероятностей реализаций возможных сценариев изменения численностей взаимодействующих популяций нормальных и аномальных клеток. С использованием данных наблюдений за их состояниями эту задачу решали посредством:

а) проведения К-кратной имитации значений {а^а^^ параметров а0 и а1 с использованием байесовских апостериорных законов распределения (17)-(20);

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

в) использования комбинаторного определения вероятности (определения Лапласа).

Имитация значений {а°,а]}к=1 проводилась в два этапа.

На первом этапе с использованием машинного генератора псевдослучайных чисел имитировались реализации {у^у}}^ для двух независимых равномерно распределенных на отрезке [0,1] случайных величин.

На втором этапе, с использованием уже известных байесовских апостериорных законов распределения параметров а0 и а1 , а также реализаций {у^у1}К=1 имитировались реализации {а^а]}^ по стандартной схеме [3]. Для этого в соответствии со свойствами многомерного распределения Стьюдента [6] и равенствами (18) были определены байесовские апостериорные маргинальные плотности распределения вероятностей параметров а1 и Ь1 :

f (a11Y4, YH4,Хнч)

M1

4S2 + [a1-a1](YH4YHT4-Y^^^^r^Y^;)[a1-a1]T}- 2

(21)

M1

f (61|Уч,Унч,ХнЧ) -

~{^2+[б1—Ь1](ХнЧХнТч—ХнчУТ (УнчУнч )-1УнчХч )[Ь1 — Ь1]т}-^

(22)

где £2 — (уч — [а1,64% )(уч — [а1,64% )т = ппт.

Для проверки предложенного метода получения оценок вероятностей реализации вариантов развития клеточных популяций был проведен вычислительный эксперимент. Параметры математической модели (1) имели следующие значения: а0 = 1,0500, а1 = 1,1000, Ь1 = 0,0500.

Для нахождение точечных оценок а0, а1, Ь1 параметров а0, а1, Ь1 на первом этапе моделирования с использованием модели (1) были получены два набора измерений нормированных численностей популяций нормальных и аномальных клеток ж (к), у (к), к = 1,..., 20. При моделировании измерений к детерминированным значениям ж (к), у (к) были добавлены случайные возмущения. Вектор случайных возмущений был получен из нормально распределенной генеральной совокупности случайных чисел с математическим ожиданием, равным нулю, и СКО, равным 0,2, реализованной с помощью генератора случайных чисел системы МаНаЪ.

С использованием (11) и (13) были получены точечные оценки параметров а0, а1 и Ь1: а0 = 0,9847, а1 = 1,0705, Ь1 = 0,0493.

а

1.4 -

1.3 -1.2 -1.1 -

0.9 -

Рис. 5. Случайные реализации оценок а0 и а1

Далее с использованием генератора случайных чисел были получены случайные реализации оценок а0 и а1, распределенные по закону Стьюдента. Параметры распределений заданы соотношениями (17) и (21).

Результаты моделирования представлены на рис. 5, и по ним получены следующие оценки вероятностей реализации описанных выше трех сценариев: для первой области р1 = 0,076, для второй — р2 = 0,670, а для третьей — р3 = 0,254.

Заключение. Проведенные исследования показали возможность использования разработанной модели для качественного анализа динамики селективного размножения популяции аномальных клеток.

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

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, проект №09-04-00948а, и проекта № 2.1.1/227 аналитической ведомственной целевой программы "Развитие научного потенциала высшей школы (2009-2010 годы)".

СПИСОК ЛИТЕРАТУРЫ

1. Динамика скорости роста, иммунофенотипа и генетическая стабильность мезенхимальных стволовых клеток костного мозга человека на ранних и поздних пассажах при культивировании ex vivo / Е.Ю.Осипова, В.А.Никитина, Т.А. Астрелина и др. // Онкогематология. - 2009, № 1. - C. 44-50.

2. Математическая модель суммарных числено стей взаимодействующих клеточных популяций / Н.П.Бочков, М.С.Виноградова, И.К.Волков и др. // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2011, № 1. -С. 18-24.

3. Волков И. К., Загоруйко Е. А. Исследование операций. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2004. - 440 с.

4. ЭльсгольцЛ. Э. Дифференциальные уравнения и вариационное исчисление. - М.: Наука, 1969. - 424 с.

5. Волков И. К. Условия идентифицируемости математических моделей эволюционных процессов по результатам дискретных косвенных измерений вектора состояния // Известия РАН. Техническая кибернетика. - 1994. № 6. - C. 55-72.

6. Зельнер А. Байесовские методы в эконометрии. - М.: Статистика, 1980. -440 с.

7. А л б е р т А. Регрессия, псевдоинверсия и рекуррентное оценивание. - М.: Наука, 1977. - 224 с.

8. Т ю р и н Ю. Н. О предельном распределении Колмогорова-Смирнова для сложной гипотезы. // Известия АН СССР. Математика. - 1984. - Т. 48, № 6. -C.1314-1343.

Статья поступила в редакцию 31.03.2011

Николай Павлович Бочков родился в 1931г., окончил Второй Московский медицинский институт им. Н.И. Пирогова в 1955 г. Д-р медицинских наук, профессор, академик РАМН, заведующий лабораторией мутагенеза Медико-генетического научного центра РАМН, лауреат Государственных премий СССР и РФ в области науки и техники, лауреат премии Правительства РФ. Автор более 500 научных публикаций в области генетики человека и медицинской генетики.

N.P. Bochkov (b. 1931) graduated from the 2nd Moscow Medical Institute n.a. N.I. Pirogov in 1955. D. Sc. (Med.), professor, academician of Russian Academy of Medical Sciences, head of laboratory of mutagenesis of Medical and Genetic Scientific Center of RAMS, winner of the USSR and RF State Prizes in Science and Technology, winner of the RF Government Prize. Author of more than 500 publications in the field of human genetics and medical genetics.

Марина Станиславовна Виноградова родилась в 1961 г., окончила Московский институт электронного машиностроения в 1984 г. Старший преподаватель кафедры "Математическое моделирование" МГТУ им. Н.Э. Баумана. Автор шести научных работ в области математического моделирования и программирования.

M.S. Vinogradova (b. 1961) graduated from the Moscow Institute of Electronic Machine Building in 1984. Senior teacher of "Mathematical Simulation" department of the Bauman Moscow State Technical University. Author of 6 publications in the field of mathematical simulation and programming.

Игорь Куприянович Волков родился в 1946 г., окончил Казанский государственный университет в 1970 г. Д-р физ.-мат. наук, профессор, лауреат премии Правительства РФ в области науки и техники. Автор более 140 научных работ в области математической физики, математического моделирования и математической генетики.

I.K. Volkov (b. 1946) graduated from the Kazan State University in 1970. D. Sc. (Phys.-Math.), professor, winner of the RF Government Prize in Science and Technology. Author of more than 140 publications in the field of mathematical simulation and mathematical genetics.

Уважаемые коллеги!

Научно-исследовательский институт гиперкомплексных систем в геометрии и физике и Трансильванский университет (г. Брашов, Румыния) приглашают Вас принять участие в Международной конференции "Финслеровы обобщения теории относительности" (FERT-2011), которая будет проходить в период с 28 августа по 4 сентября 2011 г. в г. Брашов (Румыния). Подробную информацию о конференции можно получить на сайте http://www.hyper-complex.ru/index.php?lang=rus

Оргкомитет

Dear Colleagues,

Research Institute of Hypercomplex Systems in Geometry and Physics and Transilvania University of Brasov (Romania) invite you to attend the International Conference VI International Conference "Finsler Extensions of Relativity Theory" (FERT 2011), August 28th and September 4th, 2011, which will be held in Brasov, Romania. http://www.hyper-complex.ru/index.php?lang=eng

Organizing Committee