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

НАУКА и ОБРАЗОВАНИЕ

Эл № ФС77 - 30569. Государственная регистрация №0421100025. ГББМ 1994-0408

Качественный анализ модели функционирования взаимодействующих клеточных популяций 77-30569/251409

# 10, октябрь 2011 М. С. Виноградова

УДК 51.76:517.9:57.085.23

МГТУ им. Н.Э. Баумана mathmod@bmstu.ru

1. Введение

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

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

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

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

Объектом исследований является размножающаяся in vitro культура клеток, в которой различают два типа клеток: нормальные (”здоровые”) и аномальные (”анеуплоидные”).

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

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

2. Учет влияния плотностного фактора

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

dx , А ,

— = (а — 1)х, t > 0 d t v 7

= д((а1 — 1 )y + b1x), t> 0 (1)

dt

x(0) = 1, y (0) = 0

где t — безразмерное время; x — нормированная суммарная численности популяции нормальных клеток; у — нормированная суммарная численности популяции аномальных клеток; д = т0/т1 — коэффициент, характеризующий отношение скоростей размножения аномальных и нормальных клеток, т0 и т1

— средняя длительности продолжительность клеточных циклов нормальных и аномальных клеток, а комплексные параметры а0, а1, b1 определены соотношениями

а0 = (1 — Mi0)(1 — AE0) + 2Mi0(1 — AE 0)(1 — 70),

а1 = (1 — Mi1)(1 — AE1) + 2Mi1(1 — AE1), (2)

b1 = 2Mi0 (1 — AE°)7°.

Здесь AEj — доля клеток в j-ой популяции, погибающих на временном интервале длительности т?; Mij — доля клеток в j-ой популяции, разделившихся на временном интервале длительности т?; Y0 — доля нормальных клеток, каждая из которых в процессе деления на временном интервале длительности т0 становится аномальной.

При дальнейших исследованиях будем предполагать, что AEj £ (0,1),

Vj £ {0, 1} и Y0 £ (0,1), поскольку при значениях этих параметров, равных нулю или единице, модель (1) не интересна с биологической точки зрения.

В модели (1) безразмерное время t связано с текущим временем t соотношением t = т0t, а сама модель получена из исходной дискретной модели [2] стандартными методами [4] с использованием предельного перехода при т0 ^ 0.

Известно [5-7], что борьба за ресурсы оказывает доминирующее влияние на процессы гибели и размножения. Применительно к изучаемой популяционной системе этот факт находит свое отражение(в первую очередь) в наличие функциональной зависимости доли разделившихся клеток (параметры {Mij }]=0) от суммарной численности популяций нормальных (j = 0) и аномальных (j = 1) клеток.

Развитие популяции ”in vitro” происходит в условиях ограниченного пространства и, следовательно, ограниченных питательных ресурсов, поэтому численность популяций нормальных и аномальных клеток зависит от количества клеток в единице объема, то есть от плотностного фактора.

Ограничимся линейным приближением [7,8] и с учетом естественных ограничений Mij £ (0,1), Vj £ {0,1}, примем

!0; aj — в'0Х — в 1У < 0

a — в'0Х — в j 1У ; aj — в'0Х — в j 1У £ (0,1) (3)

1; aj — ej0X — в? 1У > 1

Содержательная интерпретация параметров, входящих в правую часть ра-

венства (3), аналогична стандартной используемой для вольтерровских моделей [6-8], но с учетом специфики рассматриваемой популяционной системы:

вз — уменьшение доли делящихся клеток на временном интервале длительности т3 как за счет внутри популяционной ] = г, так и за счет межпопуляци-онной ] = г конкуренции. В силу указанной интерпретации примем вз > 0,

V г, ^ Е {0,1}. Интерпретация параметра а3- идентична интерпретации параметра Мг3, и 0 < аз < 1, V ] Е {0,1}.

Согласно классификации типов взаимодействия ”по результату” [7], для популяции нормальных клеток взаимодействие с популяцией аномальных клеток является вредным, так как в процессе деления часть нормальных клеток трансформируется в аномальные, а для популяции аномальных клеток это взаимодействие по указанной выше причине является полезным. Таким образом в изучаемой популяционной системе тип взаимодействия между образующими ее клеточными популяциями — ”хищник - жертва” [7], где в роли ’’хищника.” выступает популяция аномальных клеток = 1), а в роли ’жертвы” — популяция нормальных клеток (^ = 0).

В дальнейших рассуждениях при проведении качественного анализа математической модели (1)—(3) ограничимся областью О фазовой плоскости х0у, определяемой соотношенинями

О = О0 П О1 = 0,

00 = {(х, у) | х Е К, у Е К; 0 < ао - воох - во1У < 1; х > 0; у > 0}, (4)

01 = {(х, у) | х Е К, у Е К; 0 < а1 — в10х — впУ < 1; х > 0; у > 0}.

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

^ = НщЖ + ^02х2 + ^озху, ш

^у = ^цх + 12 у + ^1зх2 + ^14ху + ^15у2,

а£

(5)

где

Ло1 = (1 — АЕ0)(1 — 27>о — АЕ0, /02 = —воо(1 — АЕ0)(1 — 27°), Ноз = —в01 (1 — АЕ0)(1 — 270), Ни = 2Мао(1 — АЕ0)70,

/12 = м((1 — АЕ 1)а1 — АЕ1), Н13 = —2двоо(1 — АЕ0)70,

Н14 = —М(в10(1 — АЕ 1) + 2в01(1 — АЕ0)70), Н15 = —Мв11(1 — АЕ1).

(6)

Заметим, что параметры, входящие в коэффициенты Ну, меняются в следующих диапазонах: АЕ3 Е (0,1), 70 Е (0,1), ву > 0, ау Е (0,1) согласно содержательной интерпретации параметров АЕ0, АЕ1, 70, в*з, ау.

3. Точки покоя

Согласно [9] точки покоя изучаемой популяционной системы определяются неотрицательными решениями системы алгебраических уравнений, задаваемых правыми частями уравнений состояния из (5):

Н01Ж + Но2х2 + Нозху = 0, (7)

Нцх + Н12у + Н13 х2 + Н14ху + Н15у2 = 0,

Для удобства дальнейших рассуждений будем использовать следующие понятия [10]:

(1) вырожденная точка покоя — точка покоя О&(х&, у&), хоть одна из координат которой х^ или у^ является нулевой;

(2) внутренняя точка покоя — точка покоя О& (х&, у&), обе координаты которой являются положительными числами.

Вырожденные точки покоя. Из анализа системы (7) вытекает, что в общем случае изучаемая популяционная система не может иметь более двух вырожденных точек покоя в области О, определяемой соотношениями (4). При этом точка покоя О1 (0, 0) соответствует полному вымиранию популяционной системы, а точка покоя 02(0, — Н12/Н15) соответствует ее вырождению в популяцию аномальных клеток предельно возможной численности.

Согласно (6)

± ( _ ае 1 \

Н15 вп V1 1 — АЕV .

А так как ви > 0, то для существования вырожденной точки покоя О2 должно выполняться условие положительности ее второй координаты

АЕ1

а1 > 1 — АЕ1' (8)

При этом реализация ограничения (8) с учетом требования а1 Е (0,1), АЕ1 Е (0,1) возможна лишь при условии АЕ1 Е (0,1/2), поскольку в интервале (0,1) функция г/(1 — г) монотонно возрастает и равна 1 при г * = 1/2.

При AE1 £ (0,1) неравенство (8) эквивалентно неравенству

(1 — AE 1)a1 > AE1.

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

При выполнении условия

2«=^. (10)

h11 h13

реализуется особый случай — имеется третья вырожденная точка покоя

O3(—h01/h02, 0). Однако, непосредственно подставляя соответствующие выражения из (6) в (10), можно показать, что равенство (10) имеет место только при AE0 = 0. Последнее условие означает, что нормальные клетки не умирают в процессе развития популяции. При этом значении первой координаты точки O3 при AE0 = 0 равно а0/(в00), что проводит к равенству Mi0 = 0, то есть численность ”вечно живущих” нормальных клеток достигает своего предела и их размножение прекращается. При этом аномальные клетки полностью вымирают.

Заметим, что с содержательной точки зрения рассмотренная ситуация нереализуема, и в модели (1) предполагается, что AE0 > 0.

Внутренние точки покоя. При y0 = 1/2 система алгебраических уравнений (7) определяет координаты лишь вырожденных точек покоя изучаемой клеточной популяционной системы, поскольку в этом случае правая часть первого уравнения системы (5) становится тождественно равной нулю. Поэтому для определения ее внутренних точек покоя полагаем y0 £ (0,1/2) U (1/2,1).

При x = 0 первое уравнение системы (7) упрощается и принимает следующий вид:

h01 + h02x + h03y = 0. (11)

При указанных выше диапазонах изменения параметров модели и дополнительном предположении в01 = 0 согласно (6) выполняется неравенство h03 = 0, и

это уравнение можно преобразовать к виду

y = — kx + b,

где

k = h02 = в00

h03 /$01 ’ /1 ^

b = h01 = 1 \ AE0 ] (13)

h>3 &Л“0 (1 — AE0)(1 — 2y0) J .

Подставив (12) во второе уравнение системы (7) и выполнив преобразования

с учетом (6), получим квадратное уравнение

Ax2 + Bx + C = 0, (14)

где

A = h13 + h^k + h^k2 = k^e

^01

2

^01

B — hn + /112k + ^14& + 2h^bk — —a в + ¿1k — ¿0 — вюЬ, (15)

с = н12ь + н15ь2 = в11ь2 — ^16.

Здесь Ав = det(вj),

АЕ1

^ = а1 — 1 — АЕ1, (16)

270АЕ0

^ = (1 — АЕ1)(1 — 270). ( 7)

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

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

Условия существования внутренних точек покоя. Исследуем условия, при которых существуют внутренние точки покоя. Одним из необходимых

условий существования внутренней точки покоя изучаемой популяционной системы является выполнение требования Ь > 0, где комплексный параметр Ь определен в (13).

Утверждение 1. Если 7о Е (1/2,1), то Ь > 0.

Доказательство. Поскольку в00 > 0, то условие Ь > 0 эквивалентно условию

АЕ 0

> (1 — АЕ0)(1 — 270)' (18)

При этом в силу ограничения АЕ0 Е (0,1) при 70 Е (1/2,1) неравенство (18) будет удовлетворяться, поскольку правая часть неравенства будет отрицательной.

Утверждение 2. Если Ь > 0 и 70 Е (0,1/2), то АЕ0 < 1

2(1 —70)'

Доказательство. Из условия Ь > 0 с учетом ограничений на а0 следует, что

АЕ0

< а0 < 1.

(1 — АЕ0)(1 — 270)

При 70 Е (0,1/2)

откуда вытекает

АЕ0 < (1 — АЕ0)(1 — 270):

0 < 1 — 270 — 2АЕ0(1 — 70).

Из последнего неравенства с учетом диапазонов изменения входящих в него параметров можно получить условие

01

АЕ0 < 1

2(1 — 70)'

4. Исследование устойчивости точек покоя

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

Для достижения поставленной цели проведем исследование устойчивости внутренних точек покоя 01(0,0) и 02(0, — ^12/^15) по первому приближению. Корректность использованния этого подхода непосредственно следует из известных результатов [9,11].

Матрица Якоби правых частей системы обыкновенных дифференциальных уравнений (5) имеет вид

Р(х,у)=( ^ + 2Л“Х +ЛозУ ^ V (19)

^ Нц + 2^13Ж + /¿14У ^12 + ^14X + 2^15У у

где соответствующие параметры заданы соотношениями (6).

Для вырожденной точки покоя 01(0,0) матрица Якоби (19) упрощается:

P(0, 0) =

h01 0 h11 h12

и её собственные числа равны

Л1 = h01 = (1 — AE 0)(1 — 2y0)«0 — AE0, (20)

Л, = h12 = р((1 — AE')«! —AE1). ( )

Утверждение 3. Если одновременно выполняются неравенства

AE0

(1 — 27>0 < ITA1S, (21)

AE1

< 1 — AE1, (22)

то точка покоя O1 является устойчивым узлом.

Действительно, если имеют место неравенства (21), (22) то согласно( 20) Л1 = h01 < 0 и Л2 = h12 < 0, что означает справедливость утверждения.

Заметим, что

1) при y0 £ [0.5,1) неравенство (21) выполняется автоматически, и в проверке нуждается лишь неравенство (22);

2) при AE1 Е [1/2,1) условие (22) выполняется автоматически, поскольку в этом случае AE1/(1 — AE1) > 1, а а1 < 1;

3) выполнене условия (21) при Y0 Е (0, 1 /2) эквивалентно выполнению условия b < 0, а при y0 Е (1/2,1) — условию b > 0;

4) с учетом (16) условие (22) можно записать в виде d1 < 0.

Неравенство (21) эквивалентно неравенству

«0(1 — AE0)(1 — Y0) < AE0 + a0Y0(1 — AE0). (23)

С учетом (23) реализация неравенств (21), (22) обусловлена совокупностью следующих факторов:

1) доля нормальных клеток, ’’трансформирующихся” в аномальные в процессе деления, составляет менее 50 % (70 Е (0,1/2)), но в силу (23) доля нормальных клеток, которые на временном интервале разделились (без учета плотностного фактора), выжили и остались нормальными, меньше суммарной доли ’вымерших” нормальных клеток и нормальных клеток, которые в процессе деления стали аномальными, то есть популяция нормальных клеток вымирает, и при этом в силу справедливости (22) доля аномальных клеток, разделившихся и выживших, должна быть меньше доли вымерших аномальных клеток, то есть популяция аномальных клеток также является вымирающей;

2) доля нормальных клеток, ’трансформирующихся” в аномальные в процессе деления, составляет не менее 50 % (70 Е [1/2,1)), и доля аномальных клеток, разделившихся и выживших, должна быть меньше доли ’вымерших” аномальных клеток, то есть обе популяции вымирают.

Утверждение 4. Если одновременно выполняются неравенства

АЕ0

(1 — 2Т>0 > , (24)

АЕ1

> 1 — АЕ1 ’ ( )

то точка покоя 01 является неустойчивым узлом.

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

С учетом (16) условие (25) можно записать в виде ¿1 > 0.

Заметим, что неравенство (24) может выполняться только при 70 Е (0,1/2) и АЕ0 Е (0,1/2), а неравенство (25) может иметь место только при АЕ1 Е (0,1/2).

Неравенство (24) эквивалентно следующему:

ао(1 — АЕ0)(1 — 70) > АЕ0 + ао70(1 — АЕ0). (26)

Таким образом реализация неравенств (24), (25) обусловлена совокупностью следующих факторов: доля нормальных клеток, ”трансформирующихся” в аномальные в процессе деления, составляет менее 50 % (70 Е (0,1/2)),

доля ’вымерших” нормальных клеток менее 50 % (АЕ0 Е (0,1/2)), и в силу (26) доля нормальных клеток, которые на временном интервале разделились (без учета плотностного фактора), выжили и остались нормальными, больше суммарной доли ”вымерших” нормальных клеток и нормальных клеток, которые в процессе деления стали аномальными, то есть популяция нормальных клеток растет, при этом доля ”вымерших” аномальных клеток меньше 50 % (АЕ1 Е (0,1/2)) и в силу справедливости (25) доля аномальных клеток, разделившихся и выживших, должна быть больше доли ”вымерших” аномальных клеток, то есть популяция аномальных клеток также является растущей.

Утверждение 5. Если одновременно выполняются неравенства (21) и (25) или неравенства (22) и (24) то точка покоя 01 является седловой точкой.

Действительно, при выполнении условий (21), (25) имеет место А1 < 0, Л2 > 0, а при выполнении условий (22), (24) — А1 > 0, А2 < 0, откуда и вытекает справедливость утверждения.

Напомним, что седловая точка не является устойчивой.

Совокупность факторов, при которых справедливы неравенства (21), (25) или неравенства (22), (24) можно получить комбинированием соответствующих факторов, при которых выполняются утверждения 3 и 4.

Особые случаи возникают, если А1 =0 или А2 = 0 или оба собственных числа равны нулю одновременно.

Утверждение 6. Если

АЕ0

“0 = (1 — АЕ0)(1 — 2"/>) ’ ( 7)

и справедливо неравенство (25) или

АЕ1

«1 ! _ ае 1 ’ (28)

и справедливо неравенство (24), то точка покоя 01 является неустойчивой.

При выполнении условия (25) имеет место А2 > 0, а при выполнении условия (24) — А1 > 0. Следовательно [9], точка покоя 01 неустойчива.

Заметим, что равенство (27) может иметь место только при 70 Е (0,1/2) и АЕ0 Е (0,1/2), поскольку а0 Е (0,1).

Из аналогичных соображений вытекает, что равенство (28) может иметь место только при АЕ1 Е (0,1/2).

Рассмотрим случай, когда А1 = 0, а А2 < 0, то есть имеет место (27) и справедливо неравенство (22). В этом случае требуется дополнительный анализ, который можно провести согласно методике, приведенной в [9]. Система (5) в рассматриваемом случае имеет вид

(29)

где X(х, у) = ^02Ж2 + Нозху, У (х, у) = Н^х2 + НМху + Н^у2.

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

Н

11

х = х, у = -— х + у. Н12

(30)

Заметим, что Н11 > 0 из содержательных соображений, а Н12 < 0 в силу справедливости (22), поэтому данная замена определена корректно. В новых переменных система (29) имеет вид

где

¥ = Х (х,й,

^ = Н12У + У>(х,у),

Х(х,у) = ^х2 + Нозух У^(х,у) = ^2х + #зху + Н15у ,

(31)

= Н02 Н03

Н

11

Н

12

Р1

^2

$3 = Т— (Н03 — 2Н15) + Н14.

Н13 + 7— (Н02 — Н14) + Г 7—1 (Н15 — Н03),

Н12 Н12

Н

Н

12

При этом

X (х, 0) = ^1х2, У(х, 0) = д2х2.

Поскольку х(х, 0) и У(х, 0) содержат только квадраты х, то при £ = 0 и £2 = 0 согласно [9] точка покоя 01 будет неустойчивой в переменных х, у. В силу линейности замены эта точка покоя будет неустойчивой и в исходных переменных.

Случаи, когда £ = 0 и/или £2 = 0, требуют отдельного исследования, которое в силу достаточной громоздкости выражений, определяющих £ и £2, по-видимому возможно лишь при конкретных значениях параметров системы.

Рассмотрим случай, когда А1 < 0, а А2 = 0, то есть справедливо неравенство (21) и имеет место (28). Проведем анализ устойчивости точки покоя, также следуя [9].

Система (5) в рассматриваемом случае имеет вид

d x

— = hoix + X (x,y), ddy = hiix + Y (x,y),

(32)

где X(x, y) = ho2x2 + ho3xy, Y(x, y) = hi3x2 + hwxy + hi5y2.

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

hii I ~

y=- -—x + y, x = x. (33)

h0i

Поскольку h0i < 0, замена определена корректно.

В новых переменных система (32) примет вид

!=^x), ^ = hoix + X(y,x),

где

— . \ r\ r\

Y(y, x) = hi5y + giyx + g2x , X(y, x) = ho3yx + (ho2 + hoH x2,

(34)

h

oi

н

Р1 — Т— (2Н1Б — Н0з) + H14,

н

01

д2—

11

При этом

где Н1Б < 0, а

н13 + Т11 (Н14 — Н02) + Г "Г

Н01 V Н01

У (у, 0) — Н1БУ2

(Н1Б — Н03)

Х(у, 0) —0.

Согласно [9], сделать вывод о неустойчивости точки покоя (0, 0) системы (34) можно в том случае, когда Х(у, 0) = /1 ут + /2ут+1 + ..., где т > 2, однако в данном случае это условие формально нарушено. Тем не менее, следуя рекомендациям, приведенным в [9], для системы (34) можно построить знакопеременную функцию V(у, X), производная которой в силу системы (34) будет положительно определена в некоторой окрестности исследуемой точки покоя. Следовательно, по первой теореме Ляпунова о неустойчивости точка покоя (0,0) будет неустойчивой в переменных X, у. В силу линейности замены эта точка покоя будет неустойчивой и в исходных переменных.

Случай, когда оба собственных числа А1 = А2 = 0, требует отдельного исследования.

Перейдем к анализу устойчивости второй вырожденной точки — точки 02(0, —^12/^15). Матрица Якоби в этой точке определена следующим равенством

/ , \

Р(0, —Н12/Н1Б) —

Её собственные числа равны

Н01 — Н03

н

12

V

н11 — н

н1Б

н12

0

14

н

н

12

(35)

1Б

/

А1 — Н01 — н

н

12

¿031— н1Б

00

— (1 — АЕ0)(1 — 270) «0 —

АЕ0_______________в01

(1 — АЕ 0)(1 — 27°) 011

АЕ1

«1

1 АЕ1

А2 — —Н12 — — д((1 — АЕ )а1 — АЕ ).

C использованием (13) и (16) можно получить более удобное представление для первого из них:

Ai = (1 - AE°)(1 - 27°)во1 (b - вЦ .

Напомним, что из условия положительности второй координаты точки покоя O2 вытекает условие (8), которое может реализоваться только при AE1 G (0,1/2). Поэтому имеем А2 < 0.

Утверждение 7. Если d1 > 0 и справедливо одно из двух условий:

d

1) y ° G (0,1/2), b<^ ;

Pii

d

2) y ° G (1/2,1), b>^ ;

Pii

то точка покоя O2 является устойчивым узлом.

Действительно, при di > 0 вторая координата точки покоя положительна и

A2 < 0, а при справедливости первого или второго условия Ai < 0.

Утверждение 8. Если di > 0 и справедливо одно из двух условий:

d

1) y ° G (0,1/2), b>/ ;

Pii

d

2) y ° G (1/2,1), b<^ ;

Pii

то точка покоя является седлом.

При выполнении условий утверждения имеет место Ai > 0, А2 < 0, поэтому точка покоя O2 является седлом и неустойчива.

Заметим, что при выполнении условия b > —— в силу положительности di

Pii

и Pii имеет место b > 0, то есть наряду с вырожденными точками покоя Oi и

O2 могут существовать и внутренние точки покоя популяционной системы.

Трудности содержательной интерпретации условий реализуемости различных вариантов устойчивости вырожденной точки покоя 02(0, di /Pii) в основном связаны с трудностями содержательных интерпретаций комплексных параметров b и di /Pii.

Согласно (16)

d^_ /а°(1 - AEi) - AEi

Pii Pii \ 1 - AEi

где a°(1 - AEi) — доля аномальных клеток, которые на единичном временном интервале разделились и ”выжили”, AEi и (1 - AEi) — доли аномальных

клеток ’умерших” и ’выживших” на единичном временном интервале соответственно, а в11 — коэффициент, определяющий степень падения способности к делению аномальных клеток за счет внутрипопуляционной конкуренции с аномальными клетками.

Таким образом, ¿1 = а1 — АЕ1/(1 — ав1) определяет приращение аномальных клеток на единицу ’выживших” аномальных клеток (на единичном для популяции аномальных клеток временном интервале), а ^1/011 — приведенное значение этого приращения, то есть удельный коэффициент роста популяции аномальных клеток.

Выражение (13) можно представить в следующем виде

1 / ао(1 — АЕ0)(1 — 7°) — ао(1 — АЕ0)70 — АЕ0 = Л V (1 — АЕ0)(1 — 7°) — (1 — АЕ°)70

где для единичного временного интервала, характерного для популяции нормальных клеток, а0(1 — АЕ0)(1 — 70) — доля нормальных клеток, разделившихся, ’выживших” и оставшихся нормальными; а0(1 — АЕ0)70 — доля нормальных клеток, которые разделились, ’выжили” и стали аномальными; АЕ0

— доля ’’вымерших” нормальных клеток; (1 — АЕ0)(1 — 70) и (1 — АЕ0)70 — доли ’’выживших” нормальных клеток, которые остались нормальными и стали аномальными соответственно; 001 — коэффициент, определяющий степень падения способности к делению нормальных клеток за счет межпопуляционной конкуренции с аномальными клетками.

Таким образом, а0 — АЕ0/((1 — АЕ0)(1 — 270)) определяет приращение нормальных клеток на единицу ’’выживших” нормальных клеток (на единичном временном интервале), а Ь — удельный коэффициент роста популяции нормальных клеток.

С использованием введенных удельных коэффициентов можно дать содержательную интерпретацию условий реализуемости различных вариантов устойчивости точки покоя 02, задаваемых утверждениями (7) и (8).

Так, условие 2) утверждения (7) можно интерпретировать следующим образом: на любом единичном временном интервале популяция аномальных клеток имеет положительное приращение на каждую ’’выжившую” аномальную клетку (¿1 > 0), удельный коэффициент роста для популяции нормальных клеток превосходит аналогичный показатель для популяции аномальных клеток

Ь > ¿1 /в11 но в процессе деления более 50 % нормальных клеток становятся аномальными.

Рассмотрим особый случай, когда в точке покоя 02(0, — ^12/^15) имеет место А1 = 0. Этот случай реализуется при ^01 — ^03111 = 0. При этом, как

^15

указывалось выше, А2 < 0 при всех возможных значениях параметров модели популяционной системы. Сделаем замену переменных

X = X, у = у — У0, У0 = —^12/^15.

В новых переменных система (5) примет вид

= ^02Ж2 + ^0зХу = X(X, у),

^ = ^11 — ^ X — ^12У + ^1зХ2 + ^14Ху + ^15У2 = (36)

= дХ — ^12У + У(Х,у),

где

д = — ^12^14, I7(X, у) = ^1зХ2 + ^14Ху + ^15У2.

^15

Согласно [9] первое уравнение системы (36) имеет требуемый вид, однако для проведения анализа устойчивости положения равновесия необходимо преобразовать второе уравнение так, чтобы в нем не содержалось слагаемое вида дХ. Воспользуемся методом, предложенным в [9] для общего случая. Представим в виде суммы степенного ряда

у(Х) = В1X + В2Х2 + В3Х3 + ... (37)

решение уравнения

дХ — ^12 у(Х) + У(Х,у(Х)) = 0, (38)

где коэффициенты Bi, г = 1, 2,..., подлежат определению. Для их нахождения подставим (37) в (38), выпишем коэффициенты при различных степенях X и приравняем эти коэффициенты нулю. В результате получим

Г) д Г) ^13 + ^14В1 + ^15В1

Вт = ---, В2 =

^12

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

X(X, у(Х)) = (^02 + ^03Вх)Х2 + ^03В2X3 + ...

Поскольку наименьшая степень X равна 2, то при ^02 + ^03Вх = 0 точка покоя

02 будет неустойчивой.

5. Пример

Рассмотрим динамику популяции при следующих значениях параметров модели:АЕ0 = 0.200, АЕ1 = 0.275, 70 = 0.200, а0 = 0.900, а = 0.900, в00 = 0.050, 001 = 0.250, 0ю = 0.250, 0П = 0.150, ^ = 1.

Рис. 1

На рис. 1 приведен вид векторного поля системы (5). При указанных значениях параметров в системе (5) имеются три точки покоя (1, 2 и 3), отмеченные на рисунке.

При указанных значениях параметров точка 1 с координатами (0,0) является неустойчивым узлом (согласно утверждению 4), точка 2 с координатами (0,3.471) — устойчивым узлом (согласно утверждению 7), и точка 3 с координатами (2.532,1.427) — седлом. В точке 3 собственные числа имеют следующие значения: Ai ~ 0.021, А2 ~ -0.676.

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

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

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

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

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, проекты № 09-04-00948a и № 09-07-00468.

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

1. Анеуплоидия в стволовых клетках, выделенных из жировых тканей человека /Бочков Н.П., Никитина В.А., Буяновская О.А., Воронина Е.С., Гольд-

штейн Д.В., Кулешов Н.П., Ржанинова A.A., Чаушев И.Н. //Бюллетень экспериментальной биологии и медицины. 2008. Т. 146, № 9. C. 320-323.

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

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

4. Демидович В. П., Марон И. A. Основы вычислительной математики. М.: ГИФМЛ, 1960.-660 с.

5. Базыкин А.Д. Математическая биофизика взаимодействующих популяций. М.: Наука, 1985. 181 с.

6. Вольтерра В. Математическая теория борьбы за существование. М.: Наука, 1976. 286 с.

7. Ризнеченко Г. Ю. Рубин А.В. Математические модели биологических продукционных процессов. М.: Изд-во МГУ, 1993. 301 с.

8. Моран П. Статистические процессы эволюционной теории. М.: Наука, 1973.288 с.

9. Малкин И.Г. Теория устойчивости движения. М.: Наука, 1966. 530 с.

10. Волков И.К., Крищенко А.П. Качественный анализ модели развития популяции. // Дифференциальные уравнения 1996. Т. 32, № 11. С. 1457-1465.

11. Баутин Н.Н., Леонтович Е.А. Методы и приемы качественного исследования динамических систем на плоскости. М.: Наука, 1976. 468 с.

SCIENCE and EDUCATION

El № FS77 - 30569. №0421100025. ISSN 1994-0408

The qualitative analysis of model of functioning cooperating cellular populations 77-30569/251409

# 10, October 2011 M. S. Vinogradova

Bauman Moscow State Technical University

mathmod@bmstu.ru

The mathematical model of selective reproduction dynamics of abnormal cells in culture of human stem cells in vitro is designed considering influence of the density factor of sell population on cell reproduction. The model is represented by the system of the ordinary differential equations of the second order with polynomial right parts of a special kind. For this model all rest points with non-negative coordinates at which at least one coordinate is equal to zero are found. The stability analysis of these rest points is carried out at various values of parameters of model. It is shown that the system can’t have more than two rest points which both coordinates are positive. Results of numerical modeling are included.

References

1. Aneuploidija v stvolovyh kletkah, vydelennyh iz zhirovyh tkanej cheloveka /Bochkov N.P., Nikitina V.A., Bujanovskaja O.A., Voronina E.S., Gol’dshtejn D.V., Kuleshov N.P., Rzhaninova A.A., Chaushev I.N. //Bjulleten’ jeksperimen-tal’noj biologii i mediciny. 2008. T. 146, № 9. C. 320-323.

2. Matematicheskaja model’ summarnyh chislenostej vzaimodejstvujuwih kle-tochnyh populjacij. Bochkov N.P., Vinogradova M.S., Volkov I.K. i dr.// Vestnik MGTU im. N.Je. Baumana. Estestvennye nauki. 2011, № 1. S. 18-24.

3. Bochkov N.P., Vinogradova M.S., Volkov I.K. Ocenka verojatnosti realizacii variantov razvitija vzaimodejstvujuwih kletochnyh populjacij // Vestnik MGTU im. N.Je. Baumana. Estestvennye nauki. 2011, № 3. S. 31-43.

4. Demidovich V. P., Maron I. A. Osnovy vychislitel’noj matematiki. M.: GIFML, 1960.-660 s.

5. Bazykin A.D. Matematicheskaja biofizika vzaimodejstvujuwih populjacij. M.: Nauka, 19S5. 1S1 c.

6. Vol’terra V. Matematicheskaja teorija bor’by za sushhestvovanie. M.: Nauka, 1976. 2S6 s.

7. Riznechenko G. Ju. Rubin A.V. Matematicheskie modeli biologicheskih produk-cionnyh processov. M.: Izd-voMGU, 1993. 301 s.

S. Moran P. Statisticheskieprocessyjevoljucionnoj teorii. M.: Nauka, 1973. 2SS s.

9. Malkin I.G. Teorija ustojchivosti dvizhenija. M.: Nauka, 1966. 530 s.

10. Volkov I.K., Krishchenko A.P. Kachestvennyj analiz modeli razvitija populjacii. //Differencial’nye uravnenija 1996. T. 32, № 11. S. 1457-1465.

11. Bautin N.N., Leontovich E.A. Metody i priemy kachestvennogo issledovanija dinamicheskih sistem naploskosti. M.: Nauka, 1976. 46S s.