Модел на Sally Otto

Феноменология

Записването на връзките за модела има следния смисъл.

(Subscript[r, n]) е количеството налични кръвни телца в n-тия ден (а не броят унищожени такива);

(Subscript[μ, n]) е количеството новопроизведени (от костния мозък) телца;

процесът по унищожаването на стари кръвни телца ще се отчита в (Subscript[r, n]);

процесите протичат в реда (Subscript[r, n]) (унищожаване), (Subscript[μ, n]) (произвеждане), а не едновременно

Тогава след протичането на двата процеса ще имаме общо налични (Subscript[r, n])+(Subscript[μ, n]) кръвни телца за деня. Нека означим с w тази величина.

(Subscript[w, n])=(Subscript[r, n])+(Subscript[μ, n])

Първата връзка е ясна и интуитивна: от кръвните телца за миналия ден искаме да унижожим част f, но само от старите r (а не от новопроизведените μ)

(Subscript[r, n+1])=(1-f)(Subscript[r, n])+(Subscript[μ, n])

Или с други думи

(Subscript[r, n+1])=(Subscript[w, n])-f (Subscript[r, n])

Втората връзка: нека част γ от унищожените кръвни телца (f (Subscript[r, n])) да бъде възобновена от костния мозък. Интуитивно е ясно защо само за γ=1 броят на клетките няма да клони към нула или безкрайност.

(Subscript[μ, n+1])=γ f (Subscript[r, n])

Тогава параметрите се разглеждат така:

f – коефициент на унищожение на клетки от далака

γ – коефициент на производство на клетки от костния мозък

Дискретен линеен еволюционен модел

Обичайното представяне на дискретни линейни модели с еволюционен характер е:

Рекурентно задаване

v(n+1)=M v(n)

или по-общо (придвид възможното съществуване на константи във връзките):

v(n+1)=M v(n)+(Subscript[c, 0])

Общо решение

v(n)=(M^n)(Subscript[v, 0])

или по-общо (придвид възможното съществуване на константи във връзките):

v(n)=(M^n)(Subscript[v, 0])+(∑_(k=0)^(n-1))(M^k)(Subscript[c, 0])

Въвеждане на модела

За всеки случай първо изтриваме всички глобални символи (тези, които ние ползваме).

Remove["Global`*"]

Въвеждаме матрицата.

msally=(({{1-f, 1}, {f γ, 0}}))

Търсене на общото решение

Каноничната форма на една матрица M имаше вида

M=H J (H^(-1))

където H се нарича матрица на прехода, чиито стълбове са собств. вектори на M.
J е матрица, съставена от жорданови клетки за съответните собств. стойности на M.

Също така знаем, че:

(M^n)=H (J^n) (H^(-1))

Собствени вектори и стойности

Принципно прилагането на  матрица върху вектор (т.е. линейна трансформация) води до променяне на посоката и големината на вектора. Възможно е обаче матрицата да изменя големината без да променя посоката на даден вектор (или евентуално да обръща посоката). Тези вектори се наричат собствени вектори, а множителя с които се променя дължината се нарича собствена стойност.

При горната трансформация синият вектор не е собствен, тъй като променя посоката си. Червеният обаче запазва посоката си и затова е собствен. Нещо повече, тъй като запазва и големината си, то съответстващата му собствена стойност е единица.

Формално, дефиницията за собствен вектор x (имащ собствена стойност λ) на матрица A тогава е:

A x = λ x

От дефиницията за собствен вектор извеждаме и собствените стойности:

A x=λ x
A x-λ x=0
A x-λ I x=0
(A-λ I)x=0

Ако матрицата  A-λ I има обратна, то чрез умножение на двете страни по нея получаваме тривиалното решение x=0. Тогава искаме матрицата да няма обратна, т.е. детерминантата й да е нула, т.е.:

det(A-λ I)=0

Собствени стойности

Пишем матрицата M-λI.

msally-λ IdentityMatrix[2]

Използваме MatrixForm, за да видим резултата в матричен вид. Използваме постфиксен запис (//) за по-удобно.

a+b+c//f е същото като f[a+b+c]

msally-λ IdentityMatrix[2]//MatrixForm

Намираме детерминантата й. Процентът (%) замества последния изход. Ще ползваме и префиксен запис.

f@a+b е същото като f[a]+b

Det@%

Намираме решенията на характеристичния полином.

Solve[%==0,λ]

Да сравним дали Mathematica ще даде същия резултат.

Eigenvalues@msally

Да дадем имена на собств. стойности.

{λ1,λ2}=%;

За употребата на ; може да се мисли като за начин да не се изписва резултата (както в Matlab). В действителност чрез ; се образуват CompoundExpreesion-и, а в конкретния случай втория израз просто е Null и не се принтира, т.е:

1;//FullForm//HoldForm

Собствени вектори

Да зададем собствените вектори.

x1={x11,x12};
x2={x21,x22};

Намираме x1 от системата (A-(Subscript[λ, 1])I)x=0. В случая ползваме ".", за да умножим матрици и Solve за намиране на решение.

Solve[(msally-λ1 IdentityMatrix[2]).x1=={0,0},x1]

Както може да се очаква не получаваме единствено решение. При това решението е от вида „списък от правила“, т.е. A->B. Нека го заместим в x1, като изполваме ReplaceAll (/.). Тази команда позволява да се заместват елементи по списък с правила.

x1/.%

Нека нашият вектор да има за втора координата единица.

x1sol=%/.x12->1

Правим същото и за втория собствен вектор.

Solve[(msally-λ2 IdentityMatrix[2]).x2=={0,0},x2]

x2/.%

x2sol=%/.x22->1

Да сравним резултата с Mathematica.

Eigenvectors@msally

x1sol~Join~x2sol

Използвахме инфиксен запис. a~Join~b+c е същото като Join[a,b]+c

Join се използва за конкатениране на списъци (и въобще изрази с еднаква глава).

Наистина са еднакви:

%==%%

Канонична форма

Пишем матрицата H от каноничната форма. Използваме инфиксна форма (~) за по-удобно слепване на векторите.

hsally=x1sol~Join~x2sol

Да видим дали сме получили матрицата с колони – собств. вектори.

hsally//MatrixForm

Не. Трябва да я транспонираме. Можехме да получим матрицата директно (без транспониране), ако бяхме ползвали Join[...,...,2]. Тогава за H имаме

hsally=Transpose@%

Вече нещата са наред:

hsally//MatrixForm

Ще намерим обратната на H, но да видим дали H не е особена.

Det@hsally

В общия случай не е. Да намерим и обратната тогава.

hsallyinv=Inverse@hsally;

Да запишем жордановата матрица (в случая диагонална).

jsally=DiagonalMatrix@{λ1,λ2};

Проверка

Да проверим дали сме смятали правилно. Очакваме да получим msally.

hsally.jsally.hsallyinv//Simplify//MatrixForm

Анализ на модела

Аналитично тълкувание

Да припомним смисъла на параметрите. f е частта от произведените от костния мозък клетки, които далакът унищожава, т.е. 0<f<1.
γ в случая е скоростта на костния мозък – каква част от унищожените от далака клекти костния мозък ще възстанови. Тогава γ>0.

Тъй като по спектралната теорема имаме:

(M^t)=H(({{((λ1)^t), 0}, {0, ((λ2)^t)}}))(H^(-1))

видно е, че:
1.  Ако и λ1, и λ2 са отрицателни моделът няма да работи
2.  Ако |λ1|<1 и |λ2|<1, броят на телцата ще отива в нула при t->∞
3.  Ако поне едната собствена стойност е по-голяма от единица, броят на телцата ще отива в безкрайност при t->∞
4.  Очевидно двете стойности не могат да са еднакви, както вече знаем.
5.  Остава едната стойност да е единица, другата по-малка от едно. Ще видим, че има значение за коя от дрете говорим, тъй като само едната може да приеме такава стойност.

Анализ от такъв род провеждаме, използвайки функцията Reduce, която при дадени зависимости (във формата на равенстнва и неравенства), връща зависимостите изразени чрез желани променливи.

Нека проверим 1: дали собствените стойности могат да бъдат едновременно неположителни.

λ1≤0&&λ2≤0&&0<f<1&&γ>0//Reduce

Т.е. при тези предположения за параметрите и собствените стойности системата неравенства е несъвместима.

Всъщност, (Subscript[λ, 2]) е положително число, а (Subscript[λ, 1]) – отрицателно:

(λ2≤0||λ1≥0)&&0<f<1&&γ>0//Reduce

С други думи получаваме, че (Subscript[λ, 1])<(Subscript[λ, 2]):

λ2≤λ1&&0<f<γ&&γ>0//Reduce

А сега 4.: дали собстевните стойности могат да съвпадат. Предвид последната забележка, очевидно не:

λ1==λ2&&0<f<1&&γ>0//Reduce

Остава да се захванем с анализиране на петата възможност.

Какви са условията за f и γ, за да имаме

(1) λ1<λ2=1
или
(2) λ2<λ1=1

Предвид казаното малко по-горе, очевидно само (1) е възможно. Ще зададем {f,γ} като втори аргумент, тъй като искаме да получим зависимост изразена чрез тях.

Reduce[λ2==1&&λ1<1&&0<f<1&&γ>0,{f,γ}]

Т.е. условие (1) да е в сила е γ=1

Да опростим израза за λ1 при условие, че γ=1

Simplify[λ1,Assumptions->0<f<1&&γ==1]

Т.е. получаваме, че при γ=1:

λ1=-f
λ2=1

Да видим кога броят на клетките ще клони към нула.

Reduce[Abs@λ1<1&&Abs@λ2<1&&0<f<1&&γ>0,{f,γ}]

А кога ще расте в безкрайност:

Reduce[(Abs@λ1>1||Abs@λ2>1)&&0<f<1&&γ>0,{f,γ}]

Извод за собствените стойности и параметрите на аналитичното решение

Собствени числа

|λ1|<1 и |λ2|<1 – клони към 0

0<λ1=-f<1 и λ2=γ=1 – имаме решение различно от 0 и ∞

|λ1|>1 и |λ2|>1 – клони към ∞

Параметри

0<f<1

γ>0

γ<1 – клони към 0

γ=1 – имаме решение различно от 0 и ∞

γ>1 – клони към ∞

Хомеостатичност

От разглежданите модели се очаква броят на клетките да не зависи от началните данни (а по-скоро да зависи от параметрите).

За да намерим общото решение при t->∞ трябва да направим граничен преход:

Limit[hsally.(jsally^n).hsallyinv.{r0,m0}/.{λ2->1,γ->1},n->∞,Assumptions->0<f<1]

Ако съберем двата скалара ще получим граничната стойност за дневното количество кръвни телца. Използваме Apply (@@), като по този начин заменяме главата List със Plus и ефективно събираме скаларите.

Plus@@%//Simplify

Резултатът показва, че общото решение зависи от началните условия, т.е. моделът не е хомеостатичен.

Графична интерпретация

Нека зададем модела във вид на функция, за да можем да го рисуваме. За параметрите задаваме стойности по подразбиране. Използваме ReplaceAll, за да можем директно да приложим резултатите до момента. Добра практика е преди всяка дефиниция да ползваме ClearAll, за да сме сигурни, че ще се изпълни тъкмо тази дефиниция, а не някоя минала, т.е. въведена от нас преди това. Това естествено има смисъл, ако се очаква, че е възможно дефиницията да бъде променяна.

ClearAll@rmsally
rmsally[n_,r0_,m0_,myf_:0.8,myγ_:1]:=hsally.DiagonalMatrix[{((λ1)^n),((λ2)^n)}].hsallyinv.(({{r0}, {m0}}))/.{f->myf,γ->myγ}//Flatten

Графика при дадени параметри и начални условия

Да зададем конкретни стойности за f и γ за да можем нарисуваме графика.

Plot[rmsally[n,10,10],{n,0,100},PlotRange->{{0,100},{8,12}}]

Не чертае. Да видим дали стойностите, които функцията връща не са проблемът.

За дискретна стойност.

rmsally[1,10,10]

А сега за недискретна стойност.

rmsally[1.1,10,10]

Обичаен отговор за повдигане на нецелочислена степен. Тогава ще ползваме само реалната част чрез Re

Plot[Re[rmsally[n,10,10]],{n,0,100}]

Тъй като работим с дискретен модел, употребата на командата за рисуване на недискретни функции не е коректна. За тази цел се ползва DiscretePlot (изисква Mathematica 7). Подобен ефект бихме могли да постигнем и с ListPlot.

DiscretePlot[rmsally[n,10,10],{n,0,100}]

Да приближим, като зададем граници на чертежа чрез PlotRange.

DiscretePlot[rmsally[n,10,10],{n,0,100},
PlotRange->{{0,100},{8,12}}]

Графика при фиксирани начални стойности и различни параметри

Функционалният начин за създаване на цикли е чрез употребата на команди от рода на Table и Array. В долния пример показваме стойностите на модела за някои f и γ=1

Table[rmsally[1,10,10,i,1],{i,0.1,1,0.1}]

Да начертаем полученото.

DiscretePlot[
Table[rmsally[n,10,10,i,1],{i,0.1,1,0.1}],
{n,0,10}]

Не стана. Вероятно проблемът е в подадената информация. Да видим размерността й.

Dimensions[Table[rmsally[n,10,10]/.{f->i,γ->1},{i,0.1,1,0.1}]]

Има десет елемента от първо ниво и по два от второ. Нека използваме Flatten, за да бъдат резултатите само от едно ниво.

Dimensions[Table[rmsally[n,10,10,i1],{i,0.1,1,0.1}]//Flatten]

Получихме 20 елемента от едно ниво. Да пробваме да начертаем.

DiscretePlot[
Table[rmsally[n,10,10,i,1],{i,0.1,1,0.1}]//Flatten,
{n,0,10}]

Поради многото елементи чертежът стана нечетлив.

Ще ползваме Joined за да свържем точките.

DiscretePlot[
Table[rmsally[n,10,10,i,1],{i,0.1,1,0.1}]//Flatten,
{n,0,10},Joined->True]

По-добре е, но се получава нежелан фонов градиент (поради наслагване на линийките за различните функции). Ползваме Filling, за да премахнем ефекта.

DiscretePlot[
Table[rmsally[n,10,10,i,1],{i,0.1,1,0.1}]//Flatten,
{n,0,10},Joined->True,Filling->None]

Да направим същото този път за фиксирано f и да меним γ. Свързващата линия между двата клона е заради еднаквата втора стъпка за r.

DiscretePlot[
Table[rmsally[n,10,10,0.5,i],{i,0.1,1,0.2}]//Flatten,
{n,0,10},Joined->True,Filling->None]

Нека разгледаме само клона за r (първия).

DiscretePlot[
Table[rmsally[n,10,10,0.5,i][[1]],{i,0.1,1,0.2}]//Flatten,
{n,0,10},Joined->True,Filling->None]

Сега да го видим по-близо до критичната стойност на γ.

DiscretePlot[
Table[rmsally[n,10,10,0.5,i][[1]],{i,0.8,1.2,0.1}]//Flatten,
{n,0,10},Joined->True,Filling->None]

И последно за малко по-голям интервал на t.

DiscretePlot[
Table[rmsally[n,10,10,0.5,i][[1]],{i,0.8,1.2,0.1}]//Flatten,
{n,0,30},Joined->True,Filling->None]

Вижда се как при γ<1 (синята и виолетовата крива) броят на клетките клони към нула. При γ=1 (кафяво) се стабилизира (условно казано), а при γ>1 (зелено и горното синьо) отива в безкрайност.

Интерактивна графика

Използвайки Manipulate лесно можем интерактивно да меним параметрите.

Manipulate[
DiscretePlot[
Table[
Plus@@rmsally[n,r0,0,f,γ],{r0,0,10,1}
]//Flatten,
{n,0,10},Joined->True,Filling->None],
{{f,0.8},0.1,0.9,0.1},{{γ,1},0.2,2,0.2}
]

Изводи

Да резюмираме изводите на Sally Otto за модела.

Ако γ<1 клони към 0

Ако γ>1 клони към ∞

Ако γ=1

Тъй като (Subscript[λ, 1]) е отрицателно (-f), стойностите ще осцилират към равновесната стойност: (OverHat[r])=((r0+m0)/(1+f)).

Равновесие (но не в смисъла на хомеостатичност) ще се достига по-бързо за малко f и по-бавно за голямо. Тогава когато далакът е по-малко активен по-бързо ще се постига равновесие, това не е интуитивно.