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

Математика к Математическое

моделирование

УДК 51-76; 57.087.1

Получение законов распределения оценок параметров модели

популяционной системы численными методами Штраус Е. Ю.1, Ткачев С. Б.1'2'*, Волков И. К.1

Ссылка на статью:

// Математика и математическое моделирование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. №4. С. 81-92.

Б01: 10.7463/шаШш.0415.0812686

Представлена в редакцию: 28.08.2015 © МГТУ им. Н.Э. Баумана

mathmod@bmstu.ru

1МГТУ им. Н.Э. Баумана, Москва, Россия 2ФИЦ «Информатика и управление» РАН, Москва, Россия

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

Ключевые слова: динамическая система; оценка параметров; численный метод; маргинальный закон распределения; критерий Колмогорова; популяционная система

Введение

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

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

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

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

Отметим, что особенностью предложенной в [6] модели является ее нелинейность по переменным и линейность по параметрам. Методы идентификации параметров на экспериментальных выборках ограниченного объема для таких моделей известны [8, 9] и позволяют получить оценки параметров модели, а также проверить адекватность полученной модели экспериментальным данным.

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

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

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

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

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

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

Рассмотрим задачу параметрической идентификации математической модели популяци-онной системы, предложенной в [5, 6]:

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

Перейдем к решению задачи идентификации параметров модели (1). Будем предполагать, что:

а) по результатам эксперимента известны численности нормальных и аномальных клеток, измеренные через равные промежутки времени;

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

Для решения задачи параметрической идентификации будем использовать ее дискретный аналог [12]. Пусть Хк, Ук — нормированные численности нормальных и аномальных клеток в момент времени Ьк, а Хк+1, Ук+1 — нормированные численности клеток в момент времени

1. Идентификация параметров модели

¿к+ъ измерения количества клеток проводятся через равные промежутки времени ДЬ = = tk+l — ¿к.

Дискретный аналог модели (1) имеет вид

Хк+1 = "о1Хк + "о2ХА; + "о3Хк Ук,

п+1 = + ь02Х2 + "ОЛ2 + ^14Хк Ук + "О^2,

где параметры дискретной модели определяются соотношениями

й01Д£ +1, "02 = "02 ДЬ, "03 = йозД£, ^ = "цДЬ, "12 Дг + 1, "03 = &1зДг, "14 = "МД£, "15 = ^15Дг.

Отметим, что система (2) является нелинейной, однако относительно оцениваемых параметров она линейна, что удобно для решения задачи идентификации.

Для решения задачи параметрической идентификации математической модели (2) воспользуемся методом, предложенным в [8, 9, 12], и будем решать задачу идентификации параметров модели (2) по ограниченной выборке результатов измерений в два этапа. На первом этапе найдем точечные оценки параметров модели, на втором — проведем проверку статистической гипотезы о том, что закон распределения случайных ошибок измерений является нормальным. Если эта статистическая гипотеза будет принята, то модель с найденными оценками значений параметров признается адекватной.

"01 "12

2. Точечные оценки параметров нелинейной модели

Следуя [12], приведем основные соотношения для решения задачи идентификации. Обозначим наблюдаемые значения переменных через Х, У и запишем модель (2) в матричном

виде:

/

"01

"02 "03

1,(1. 1,(1. 1,(1. 1,(1. 1,(1. "11 "12 "13 "14 "15

\

Х2к-1

У2к-1 Х 2

Х2к-1У2к-1

У 2 У 2к— 1

\

+

/

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

е =

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

Пусть ¿п = пДЬ, п = 1, ..., , — моменты времени, а {¿2п-1; ¿2га) — п-й двухточечный шаблон. В узлах этого шаблона известны экспериментальные значения Х2п-1, У2п-1 и Х2п, У2п численностей популяции нормальных и аномальных клеток. Согласно

Хч = Н0 +

0

0

е

п

где

Хч = Х2, Х4, . . , X2N^J , Хнч = (Xb X3, . . . , X2N-1) H0d =

Х1 Хэ . . . X2N-1 ^ ( X Хнч

ZХ = Х12 Хэ2 X 2 . . . X2N-1 =

X3Y3 . . . X2N-1Y2N-1 X V у Х нч1 нч

hd hd h h01, '02, h

ьоз

Аналогично

Уч = HdZy + n,

где

Y4

= Y2

Y2,

1

2N

Zy —

V

Y нч = (1, 1з, . . . , Y2N-1) , H1d = (hdn,

X1 Хз ... X2N -1 \ X

11 1 з . . . Y2N -1 Y нч

X12 X ... X 2 X2N-1 = ^ч

X111 X3Y3 ... X2N -1Y2N - 1 ^нч^нч

V2 Y 2 Y3 ... Y 2 12N -1 / v 2 \ 1нч

h11, h12, h13, h14, h

Оценки векторов параметров Hd, Hd имеют вид

Я" d _ Z Z+ Hd _

0 = X4Zx , H1 —

где — псевдообратные матрицы к матрицам ¿х и '¿у соответственно.

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

е — Хч — n = Y4 — H dZy,

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

Поскольку оценки получены с использованием выборок малого объема, предлагается воспользоваться критерием Колмогорова в модификации Ю.Н. Тюрина [13]. Напомним, что в критерии Колмогорова в качестве меры расхождения между теоретическим и эмпирическим распределением рассматривают максимальное значение абсолютной разности между теоретической функцией распределения F(х) и эмпирической функцией распределения Еп (х):

Dn = sup |F(x) — Fn(x) |.

Х

Для малых выборок используют следующую форму статистики:

Лг = Dn  (VN — 0.01 + Vf).

Критическое значение Аа при заданном уровне значимости а, определяют по таблице, которая приведена, например, в [13].

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

3. Нахождение маргинальных законов распределения оценок параметров

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

При найденных значениях параметров модели определим траекторию системы (2), которая согласуется с результатами наблюдений. Для этого необходимо решить задачу идентификации начальных данных. Здесь также будем использовать метод наименьших квадратов. Найдем такие начальные значения Х(0), У^(0), при которых расчетная траектория, реализующая решение задачи Коши для системы (2) с векторами параметров Н°, Н°, имеет наименьшее (в смысле метода наименьших квадратов) отклонение от результатов наблюдений.

Полученную траекторию будем называть опорной. Вычислим значения численностей популяций в моменты времени Ьк = кДЬ, к = 1, ..., 2Ж. Найдем отклонения расчетных значений от экспериментальных, определим дисперсии соответствующих случайных величин — отклонений по Х и У. Эти дисперсии будем использовать в качестве параметров для генерации нормально распределенных случайных ошибок с помощью датчика нормально распределенных случайных чисел.

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

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

4. Численный эксперимент

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

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

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

Фактически ошибки измерений возникают уже в процессе развития популяционной системы, поскольку оценка количества производится на основе анализа контрольного объема. С учетом сделанного предположения примем начальные условия для модельной траектории Х(0) = 1, У (0) = 0. Среднее квадратичное отклонение а случайной величины с нулевым начальным ожиданием, используемой для моделирования шума измерений, была выбрана равной 0.05. Ошибки измерения были приняты одинаковыми для нормальных и аномальных клеток.

На рис. 1 приведена траектория, полученная из детерминированной модели (2) при указанных выше значениях параметров и стандартных начальных данных Х(0) = 1, У (0) = 0. Траектория, приведенная на рис. 2, получена путем добавления к значениям переменных Хк и Ук в узлах сетки значений нормально распределенных случайных величин с нулевым начальным ожиданием, с помощью которых моделировались ошибки измерений.

Рис. 1. Исходная траектория системы Рис. 2. Траектория системы

при наложении шума

На рис. 3 представлены две траектории: траектория системы с наложенными на нее шумами, которые моделируют ошибки измерений, и траектория системы с найденными оценками параметров.

Стартовая точка второй траектории найдена в результате решения задачи идентификации начальных данных.

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

1.4 1.2

1

0.8 0.6 0.4 0.2 0

1 2 3 4 5 6 7 8 X

Рис. 3. Траектория системы при найденной оценке начальных условий

Дальнейшие расчеты проводились в соответствии с предложенной в разд. 3 схемой — к значениям нормированных численностей популяций в узлах прибавлялись сгенерированные реализации случайных ошибок. Добавляемые случайные величины имели нулевые математические ожидания и дисперсии, найденные при анализе невязок (8).

В результате проведенного вычислительного эксперимента при уровне значимости 0.1 оказалось возможным принять гипотезы о том, что все маргинальныезаконы распределения параметров — нормальные.

Заключение

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

Согласно известным теоретическим оценкам [11,12], полученным на основе баесовского подхода, маргинальные распределения оценок параметров являются распределениями Стъ-юдена. Оказалось, что с приемлемой точностью в качестве маргинальных распределений можно использовать нормальные распределения.

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

Работа выполнена при финансовой поддержке РФФИ, проект № 13-07-00720.

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

1. Бочков Н.П., Никитина В.А., Воронина Е.С., Кулешов Н.П. Методическое пособие по тестированию клеточных трансплантатов на генетическую безопасность // Клеточные технологии в биологии и медицине. 2009. № 4. C. 183-189.

2. Duesberg P., Mandrioli D., McCormack A., Nicholson J.M. Is carcinogenesis a form of speciation? // Cell Cycle. 2011. no. 10. P. 2100-2114.

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. supp. 01. P. 871-892. DOI: 10.1142/S0218202511005404

4. 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

5. Виноградова М.С. Динамическая модель клеточной популяционной системы // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2013. №12. С. 175-192. DOI: 10.7463/1213.0646463

6. Виноградова М.С. Исследование нелинейной модели развития клеточной популяционной системы // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2014. № 8. С. 123-138. DOI: 10.7463/0814.0720269

7. Виноградова М.С. Анализ сценариев развития клеточной популяционной системы // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журнал. 2014. № 11. DOI: 10.7463/1114.0735732

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

9. Волков И.К. Идентифицируемость математических моделей эволюционных процессов // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2005. № 3. C. 64-73.

10. Бочков Н.П., Виноградова М.С., Волков И.К. Оценка вероятности реализации вариантов развития взаимодействующих клеточных популяций // Вестник Моск. гос. техн. ун-та им. Н.Э. Баумана. Сер. Естественные науки. 2011. №3. С. 31-43.

11. Виноградова М.С. Параметрическая идентификация модели взаимодействующих клеточных популяций на основе байесовского подхода // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2012. №11. С. 155-182. Режим доступа: http://technomag.edu.ru/doc/490900.html, дата обращения: 20.08.2015.

12. Виноградова М.С. Математическое моделирование динамики развития изолированной клеточной популяционной системы: Дисс. . .. канд. физ.-мат. наук. М., 2012. 128 с.

13. Тюрин Ю Н. Непараметрические методы статистики. М.: Знание, 1978. 64 с.

Mathematics i Mathematical Modelling

Electronic journal of the Bauman MSTU

Mathematics and Mathematical Modelling of the Bauman MSTU, 2015, no. 4, pp. 81-92.

DOI: 10.7463/mathm.0415.0812686

Received: 28.08.2015

© Bauman Moscow State Technical University

http://mathmjournal.ru

Investigating Distribution of the Model Parameters Estimates for the Population System Using Numerical Methods

Strauss E. Yu.,1,Tkachev S. B.12 *, Volkov I. K.1

mathmod@bmstu.ru

1 Bauman Moscow State Technical University, Russia 2FRC "Computer Science and Control" of RAS, Russia

Keywords: dynamical system, parameter estimating, numerical method, the marginal distribution law, Kolmogorov test, population system

The paper is aimed at determining the marginal distributions of the estimates of model parameters for a nonlinear dynamic model of the cell population dynamics. The cell population system evolves in the laboratory (in vitro) and comprises two types of cells.

Basic estimates of model parameters are obtained on a single limited sample of experimental data. The value of the number of populations of each cell type is obtained from the experiment at equal intervals of time. This paper proposes a technique for determining the marginal distributions of the parameter estimates using numerical modeling.

This technique includes the identification of parameters of nonlinear models and test of the obtained model adequacy with base parameter estimates, the identification of the initial data and finding the reference trajectory. The initial data for the trajectory are found using the least squares method, while minimizing the deviation from the experimental trajectory. Data sets on measured densities of the populations at specific points in time are generated using the reference trajectory and the normally distributed random numbers generator. The problem of obtaining estimates of system parameters is solved for each data set. Tests of hypotheses about the marginal distribution of each parameter are carried out based on the calculated set of estimated parameters. To prove hypothesis, the Kolmogorov test is used. The description of a numerical example is included. The obtained marginal distributions of the parameter estimates can be further used to evaluate the probabilities of different scenarios of the population system development.

References

1. 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).

2. Duesberg P., Mandrioli D., McCormack A., Nicholson J.M. Is carcinogenesis a form of speciation? Cell Cycle, 2011, no. 10, pp. 2100-2114.

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. supp. 01, pp. 871-892. DOI: 10.1142/S0218202511005404

4. 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, pp. 913920. DOI: 10.1039/c2mb05439f

5. 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

6. 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

7. 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, p. 607-622. DOI: 10.7463/1114.0735732

8. Volkov I.K. Conditions for the Identifiability of Mathematical Models of Evolutionary Processes from the Results of Discrete Indirect Measurements of the State Vector. Izvestija RAN. Teorija i sistemy upravlenija, 1994, №6, pp. 55-72.

9. Volkov I.K. Identifiability of Mathematical Models of Evolutionary Processes. VestnikMGTU im. N.E. Baumana. Ser. Estestvennye nauki = Herald of the Bauman MSTU. Ser. Natural science, 2005, no. 3, pp. 64-73.

10. 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.

11. Vinogradova M.S. Parametrical identification of cooperating cellular populations model on the basis of Bayesian inference. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2012. №11. pp. 155-182. Available at: http://technomag.edu.ru/doc/490900.html, accessed: 20.08.2015.

12. Vinogradova M.S. Matematicheskoe modelirovanie dinamiki razvitija izolirovannojkletochnoj populjatsionnoj sistemy. Kand. diss. [Mathematical modeling of the dynamics of an isolated cell population system. Cand. diss.]. Moscow, 2012. 128 p.

13. Tyurin Yu.N. Nepametricheskie metody statistiki [Nonparametric methods of statistics]. Moscow, Znanie publ., 1978. 64 p.