Модел на Кимура

Двупараметричен модел на Кимура

Два типа бази, всеки с по два вида:

Пурини

Аденин (A)

Гуанин (G)

Пирамидини

Тимин (T)

Цитозин (C)

Мутации.

Транзиция (междувидово в един тип)

A ↔ G

T ↔ C

Трансверсия (междутипово): при човека три пъти по-малък шанс

A ↔ T

A ↔ C

G ↔ T

G ↔ C

Основни връзки. Получаване на матрицата

Нека въведем n – вектор на състоянието; показва броят нуклеотидни бази в даден момент (и в даден отрязък ДНК евентуално)

n(t)=(((({{(Subscript[n, A])(t), (Subscript[n, G])(t), (Subscript[n, C])(t), (Subscript[n, T])(t)}})))^T)

Имаме връзките:

(Subscript[n, A])(t)+(Subscript[n, G])(t)+(Subscript[n, C])(t)+(Subscript[n, T])(t)=const, ∀t

Т.е. броят на базите е константен във всеки един момент. Това е естествено следствие от предположението, че при мутация базите не изчезват, а само минават от един вид в друг.

Да въведем още:
α – вероятност за транзиция;
β – вероятност за трансверсия;

От казаното до тук следва, че:

0<β≪α≪1

Ако приемем, че моделът трябва да е линеен, то законът за еволюция на състоянията n, ще е от вида (имаме хомогенност):

n(t+1)=M n(t)

Тук M е 4×4 матрица и се нарича матрица на Кимура (виж. Motoo Kimura & James F. Crow).

За да покажем вида на матрицата, нека сформираме зависимостите (1) и (2) в математична форма. Прилагайки (2) имаме системата:

(Subscript[n, A])(t+1)=(Subscript[x, A])(Subscript[n, A])(t)+α (Subscript[n, C])(t)+β (Subscript[n, G])(t)+β (Subscript[n, T])(t)
(Subscript[n, C])(t+1)=α (Subscript[n, A])(t)+(Subscript[x, C]) (Subscript[n, C])(t)+β (Subscript[n, G])(t)+β (Subscript[n, T])(t)
(Subscript[n, G])(t+1)=β (Subscript[n, A])(t)+β (Subscript[n, C])(t)+(Subscript[x, C]) (Subscript[n, G])(t)+β (Subscript[n, T])(t)
(Subscript[n, T])(t+1)=β (Subscript[n, A])(t)+β (Subscript[n, C])(t)+β (Subscript[n, G])(t)+(Subscript[x, T]) (Subscript[n, T])(t)

Например първият ред означава: броят на базите Аденин в следващият момент е сума от мутиралите в аденин: цитозин (с вероятност α понеже е от същия тип), гуанин и тимин (с вероятност β понеже са от другия тип) плюс количеството (Subscript[x, A]) непроменили се бази аденин. По същия начин разчитаме и другите редове.

За да намерим хиксовете, трябва да приложим и връзка (1). Т.е. имаме:

(Subscript[n, A])(t+1)+(Subscript[n, G])(t+1)+(Subscript[n, C])(t+1)+(Subscript[n, T])(t+1)=(Subscript[n, A])(t)+(Subscript[n, G])(t)+(Subscript[n, C])(t)+(Subscript[n, T])(t)

Заместваме с полученото в (4) и приравняваме:

((Subscript[x, A])+α+2β)(Subscript[n, A])(t)+((Subscript[x, C])+α+2β)(Subscript[n, C])(t)+((Subscript[x, G])+α+2β)(Subscript[n, G])(t)+((Subscript[x, T])+α+2β)(Subscript[n, T])(t)=(Subscript[n, A])(t)+(Subscript[n, G])(t)+(Subscript[n, C])(t)+(Subscript[n, T])(t)

Тогава:

(Subscript[x, A])=(Subscript[x, C])=(Subscript[x, G])=(Subscript[x, T])=1-α-2β

Окончателно получаваме системата зависимости:

(Subscript[n, A])(t+1)=(1-α-2β)(Subscript[n, A])(t)+α (Subscript[n, C])(t)+β (Subscript[n, G])(t)+β (Subscript[n, T])(t)
(Subscript[n, C])(t+1)=α (Subscript[n, A])(t)+(1-α-2β) (Subscript[n, C])(t)+β (Subscript[n, G])(t)+β (Subscript[n, T])(t)
(Subscript[n, G])(t+1)=β (Subscript[n, A])(t)+β (Subscript[n, C])(t)+(1-α-2β) (Subscript[n, G])(t)+β (Subscript[n, T])(t)
(Subscript[n, T])(t+1)=β (Subscript[n, A])(t)+β (Subscript[n, C])(t)+β (Subscript[n, G])(t)+(1-α-2β) (Subscript[n, T])(t)

Непосредствено се проверява, че тогава M ще има вида:

(({{1-α-2β, α, β, β}, {α, 1-α-2β, β, β}, {β, β, 1-α-2β, α}, {β, β, α, 1-α-2β}}))

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

Нека въведем модела за последователни манипулации:

kimura:=(({{1-α-2β, α, β, β}, {α, 1-α-2β, β, β}, {β, β, 1-α-2β, α}, {β, β, α, 1-α-2β}}))

kN[t_]:={nA[t],nC[t],nG[t],nT[t]}

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

Кратко припомняне на смятане на собств. вектори и стойности

Характеристичен полином

CharacteristicPolynomial[kimura,λ]//Simplify

Нека

kB[x_]:=kimura-x IdentityMatrix[4]

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

Det@kB@λ//Simplify

За собств. стойности според Mathematica

kval=Eigenvalues[kimura]

Или по нашите сметки.

Solve[%%==0,λ]

Въпрос: собствената стойност (Subscript[λ, 3,4])=1-2α-2β двукратна ли е?

По теоремата за ранга имаме:

kB@kval[[3]]//MatrixRank

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

Собствени стойности според Mathematica

kH=Transpose@Eigenvectors@kimura

%//MatrixForm

Да пробваме ръчно да намерим (Subscript[x, 2]):

x2={x21,x22,x23,x24};

Solve[kB@kval[[2]].x2==0,x2]

За (H^(-1)): според теорията имаме:

H=(1/(det(H)))(Subscript[A, H]^T)

където (Subscript[A, H]) е адюнгираната на H, т.е.

(Subscript[A, H])={(Subscript[a, i j]):(Subscript[a, i j])=(((-1))^(i+j))(Subscript[Δ, i j])}

където (Subscript[Δ, i j]) е поддетерминантата на H, в която i-тия ред и j-тия стълб не участват

myinv[x_]:=
Transpose@
Table[
Table[
(1/(Det@x))(((-1))^(i+j))Det@x[[Range[i-1]~Join~Range[i+1,Length@x]]][[All,Range[j-1]~Join~Range[j+1,Length@x]]]
,{j,1,Length@x}
]
,{i,1,Length@x}
]

Според Mathematica

kHinv=Inverse@kH

%//MatrixForm

Според нас:

myinv@kH//MatrixForm

Нека за по-красиво вземем:

kHinv2=4kHinv;

Имаме жордановата форма:

Print["M="MatrixForm@kH,".",MatrixForm@DiagonalMatrix@kval,".(1/4)",MatrixForm@kHinv2]

Асимптотично решение

Както знаем за каноничната форма (в нашия случай с диагонална матрица) имаме:

(M^t)=H (D^t) (H^(-1))

Записвайки горното развитие като се абстрахираме от точния вид на собств. стойности за D имаме:

kDt=DiagonalMatrix[{1,((λ2)^t),((λ34)^t),((λ34)^t)}]

Или получаваме (като 1/4 е изведена предварително пред скоби):

kM1=kH.kDt.kHinv2

%//MatrixForm

Но имаме:

0<(Subscript[λ, 2])=1-4β<1
0<(Subscript[λ, 3,4])=1-2α-2β<1 (при условие, че 0<α+β<1)

Тогава:

(Subscript[λ, 2]^t)(Overscript[⟶, t->∞])0∧(Subscript[λ, 3,4]^t)(Overscript[⟶, t->∞])0

Т.е.

kM2=kM1/.{((λ2)^t)->0,((λ34)^t)->0}

%//MatrixForm

Умножавайки по n и по (да не забравяме) (1/4) получаваме:

(1/4)kM2.kN[t]//MatrixForm

Т.е. с времето количествата бази ще бъдат равни.

Визуално представяне

kModel[t_,α1_,β1_,n0_]:=MatrixPower[kimura/.{α->α1,β->β1},t].n0

kn0={400,700,0,100};

Manipulate[
BarChart[
kModel[t,α,β,kn0],BarSpacing->1,PlotRange->{{0,7},{-20,Max@kn0}},ChartLabels->{"A","G","C","T"}
]
,{t,0,512,1}
,{{α,0.1},0,1,0.05}
,{{β,0.005},0,1,0.005}
]

Маркови процеси

Въведение

Нека имаме матрица  (SuperStar[P])=(Subscript[{(Subscript[((SuperStar[p])), i j])}, m×m]). Елементът (Subscript[p, i j]^( * )) дава вероятността от състояние i да се премине в състояние j. Всъщност говорим за матрица на съседство, тъй като процесът се описва тъкмо с граф. За да може обаче матрицата да обрисува процес на Марков, необходимо е сумите на елементите от редовете да са единици, тъй като сборът от вероятностите за излизане от дадено състояние трябва да е единица. С други думи:

(∑_(j=1)^m)(Subscript[p, i j]^( * ))=1,j=1,...,m

Вероятността за преминаване от едно състояние в друго в нашия случай представлява именно матрицата на Кимура. Условие (9) е спазено.

Разглеждаме материална система от n на брой елемента S={(Subscript[e, 1]),...,(Subscript[e, n])}. За горния пример това означаваше система от 1200 елемента (1200 бази). Всеки нейн елемент може да заема едно дадено състояние σ={(Subscript[σ, 1]),...,(Subscript[σ, m])}. В нашия случай σ={A,C,G,T}. Системата S се състои от m на брой подсистеми (Subscript[S, i]), всяка от които включва елементите на S в едно от възможните състояния. Тогава нека (Subscript[x, i]) да е количествената мярка на i-тата подсистема, т.е. (Subscript[x, i]) показва броя елементи на (Subscript[S, i]). Нека x е векторът от тези мерки. По принцип (Subscript[x, i]) показва връка между число и съответната подсистема и може да има и други интерпретации.

Например в нашия случай системата изглежда така:

"S"[
"(Subscript[S, A])"[
"(Subscript[e, 1])","(Subscript[e, 4])","(Subscript[e, 13])","..."
],
"(Subscript[S, C])"[
"(Subscript[e, 7])","(Subscript[e, 8])","(Subscript[e, 23])","..."
],
"(Subscript[S, G])"[
"(Subscript[e, 2])","(Subscript[e, 5])","..."
],
"(Subscript[S, T])"[
"(Subscript[e, 6])","(Subscript[e, 10])","..."
]
]//TreeForm

Тогава, ако за една материална система имаме процес на Марков, то количествените й мерки x(t+1) в даден момент t ще са както следва:

x(t+1)=(SuperStar[P])x(t)

И следователно:

x(t)=((((SuperStar[P])))^t)(Subscript[x, 0]),(Subscript[x, 0])=x(0)=const

Нека положим:

P(t)=((((SuperStar[P])))^t)(Subscript[x, 0])

Оттук нататък ще изследваме поведението на оператора P(t). Необходимо (и естествено) е да въведем следните изисквания:

1.P(0)=I
2.P(t+s)=P(t).P(s),∀t,s≥0
3.P(t+s)=P(s+t),т.е. P(t).P(s)=P(s).P(t)

Коментар:

x(0)=P(0)x(0)⇔(P(0)-I)x(0)=0⇔P(0)=I

(2) Нарича се равенство на Чапман-Колмогоров, което е и основен принцип в динамичните системи (каквато е и нашата). То гласи, че системата ще еволюира еднакво за време t + s започвайки от 0 и при интервала s с начално t. (Тъй като уравнението се чете: прилагаме P(t) към P(s).
(3) Условието придава единственост (линейност) на еволюционния процес – независимо от пътя на еволюцията (т.е. състоянията през които минава системата) за еднакъв интервал от време и еднакви начални количествени мерки системата еволюира до едно и също свое състояние.

Фамилия от еднопараметрични линейни оператори {P(t), t≥0}, за които са изпълнени (13) се нарича еднопараметрична полугрупа.

Нека проверим дали условията са в сила за нашия пример. Имаме:

n(t)=(M^t)n(0)

Нека:

P(t)=(M^t)

Т.е.

n(t)=P(t)n(0)

Имахме:

(1/4)kM1//MatrixForm

Първо условие:

(1/4)kM1/.{t->0}//MatrixForm

Второто условие:

ReplaceAll[(1/4)kM1,t->t+s]==((1/4)kM1).((1/4)kM1/.{t->s})//Reduce

Третото условие очевидно е в сила.

Диференцируеми операторни функции

Прочетете материала за Кимура, стр. 12+

(Overscript[P, •])(t)(Overscript[=, def])(Overscript[P, •])(0)P(t)

Нека

A:=(Overscript[P, •])(0)

Тогава:

(Overscript[P, •])(t)=A P(t)

Теорема на Поанкаре

Нека имаме еднопараметричната група от линейни оператори:

{P(t):(ℝ^m)->(ℝ^m), такива че: (Overscript[P, •])(t)=A P(t), P(0)=I,(Overscript[P, •])(0)=A}

Когато (20) е диференцуема, A е диагонализуем оператор, то P(t) е от експоненциален тип.

P(t)=(e^tA),A=(Overscript[P, •])(0)

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

Първо нека вземем жордановото разлагане за P. Имахме:

(1/4)MatrixForm@kM1

Второ: нека намерим 4.A (за по-удобно с (1/4)). За тази цел диференцираме поелементно по t в t=0

(∂_t) kM1/.t->0//MatrixForm

kA=%;

Да построим фундаменталната за kA. За целта намираме собств. вектори и стойности на kA.

Eigensystem@((1/4)kA)

Тогава нейната фундаментална е:

kFA=(({{1, -(E^(Log@λ2 t)), 0, -(E^(Log@λ34 t))}, {1, -(E^(Log@λ2 t)), 0, (E^(Log@λ34 t))}, {1, (E^(Log@λ2 t)), -(E^(Log@λ34 t)), 0}, {1, (E^(Log@λ2 t)), (E^(Log@λ34 t)), 0}}))

%//MatrixForm

Тогава за P(t)=(e^(t A))=Φ(t).((Φ(0))^(-1)) имаме:

kFA.(Inverse@kFA/.t->0)//MatrixForm

Което е и тъкмо P(t).

Асиметрични апроксимации

Нека забележим, че за собствените стойности (Subscript[μ, i]) на kA получихме:

(Subscript[μ, i])=ln (Subscript[λ, i])

Развивайки тези собствени стойности в ред на Маклорен (Тейлър около нулата), получаваме приближенията:

(Subscript[μ, i])=ln (Subscript[λ, i])=(∑_(k=1)^∞)(((-1))^(k-1))((((1-(Subscript[λ, i])))^k)/k)

Тук условието е:

|1-(Subscript[λ, i])|<1

Което макар и от математическа гледна точка не винаги да е изпълнено, можем да смятаме, че практически е в сила.