Программа, реализующая стохастическую модель

alphabet = ['A', 'G', 'T', 'C'];

alphabet1 = [1 2 3 4];

Npop = 1000; %численность популяции

Ngene = 1800; %длина генома

Pmut = 0.2; %вероятность мутации

Prec = 0.5; %вероятность рекомбинации

Kiter = 1000; %количество поколений

Npairs = Npop*(Npop-1)/2; %количество пар геномов

HamGlobal = zeros(1, Kiter); % для хранения расстояния Хэмминга

MutSpeedA=[2.27743125852128, 0.318614435046603, 0.542883200161911];%A-G,A-T,A-C

MutSpeedG=[7.36088882597026, 0.655035702037978, 0.261999891600743];%G-A,G-T,G-C

MutSpeedT=[0.809503738020646, 0.58021101931253, 2.8093006895185];%T-A,T-G,T-C

MutSpeedC=[1.85795724212526, 0.228442777406494, 4.60051637685454];%C-A,C-G,C-T

%скоростьмутации

 

RouletteA=[MutSpeedA(1)/sum(MutSpeedA), (MutSpeedA(1)+MutSpeedA(2))/sum(MutSpeedA), 1];

%A->(G,T,C)

RouletteG=[MutSpeedG(1)/sum(MutSpeedG), (MutSpeedG(1)+MutSpeedG(2))/sum(MutSpeedG), 1];

%G->(A,T,C)

RouletteT=[MutSpeedT(1)/sum(MutSpeedT), (MutSpeedT(1)+MutSpeedT(2))/sum(MutSpeedT), 1];

%T->(A,G,C)

RouletteC=[MutSpeedC(1)/sum(MutSpeedC), (MutSpeedC(1)+MutSpeedC(2))/sum(MutSpeedC), 1];

%C->(A,G,T)

%расчет секторов рулетки для мутации

 

Type=['WT','M41L','T215N','T215S','T215Y','M41L+T215N','M41L+T215S','M41L+T215Y'];

Func0=[1, 0.60398, 0.20344, 0.25286, 0.70386, 0.29144, 0.49304, 0.78211];

Func03=[0.16226, 0.16627, 0.011927, 0.01554, 0.29251, 0.038069, 0.11565, 0.38495];

%функцииприспособленности

 

Roulette0 = zeros(1,8);

Roulette0(1) = Func0(1)/sum(Func0);

Roulette03 = zeros(1,8);

Roulette03(1) = Func03(1)/sum(Func03);

for i = 2:8

Roulette0(i) = Roulette0(i-1)+ Func0(i)/sum(Func0);

Roulette03(i) = Roulette03(i-1)+ Func03(i)/sum(Func03);

end

%рассчет рулетки для функции приспособленности

 

GlobType = zeros(Kiter, 8);

 

%% создание начальной популяции

pop1 = randsrc(Npop, Ngene, alphabet1); %функция создает матрицу заданных

%размеров, в которой pop(i,j) заполняется случайным числом от 1 до 4

%с равной вероятностью

pop = alphabet(pop1); % числа [1 2 3 4] -> буквы ['A', 'G', 'T', 'C']

clear pop1%удалим временно созданную матрицу чисел

pop(1, 121:123) = ['A''T''G'];

pop(1, 643:645) = ['A''C''C'];

for i = 2:Npop

pop(i, :) = pop(1, :);

end

newpop = pop; %матрица, в которую будет отбираться новое поколение

 

for k = 1:Kiter

 

%% рекомбинация

i = 1;

while i<=(Npop-1)

RandRec = rand; %(0, 1)

if RandRec<=Prec

RandGene = randi(Ngene); %[1, 1800]

Temp = pop(i, RandGene:Ngene);

pop(i, RandGene:Ngene) = pop(i+1, RandGene:Ngene);

%pop(i+1, :) = pop(i, :);

pop(i+1, RandGene:Ngene) = Temp;

%pop(i+3, :) = pop(i+2, :);

clear Temp;

i = i + 2; %если произошла рекомбинация со следующим, пропускаем его

else

i = i + 1; %если нет, то продолжаем со следующего

end

 

end

 

%% мутация

for i = 1:1:Npop

RandMut = rand;

if RandMut<=Pmut

RandGene = randi(Ngene); %номермутирующегогена

RandMutGene = rand;

switch pop(i, RandGene)

case'A'%мутирует A

if RandMutGene<=RouletteA(1)

pop(i, RandGene) = 'G' ;

elseif RandMutGene<=RouletteA(2)

pop(i, RandGene) = 'T';

else

pop(i, RandGene) = 'C';

end

case'G'%мутирует G

if RandMutGene<=RouletteG(1)

pop(i, RandGene) = 'A' ;

elseif RandMutGene<=RouletteG(2)

pop(i, RandGene) = 'T';

else

pop(i, RandGene) = 'C';

end

case'T'%мутирует T

if RandMutGene<=RouletteT(1)

pop(i, RandGene) = 'A' ;

elseif RandMutGene<=RouletteT(2)

pop(i, RandGene) = 'G';

else

pop(i, RandGene) = 'C';

end

case'C'%мутируетС

if RandMutGene<=RouletteC(1)

pop(i, RandGene) = 'A' ;

elseif RandMutGene<=RouletteC(2)

pop(i, RandGene) = 'G';

else

pop(i, RandGene) = 'T';

end

end

end

end

 

%% пересчет

NType=zeros(1,8); %количество генов каждого типа в данной популяции

for i = 1:1:Npop

if strcmp(pop(i, 121:123),'TTG')

switch pop(i, 643:645)

case'AAC'% (TTG AAC)- тип M41L+T215N

NType(6)=NType(6) + 1;

case'TCC'% (TTG TCC)- тип M41L+T215S

NType(7)=NType(7) + 1;

case'TAC'% (TTG TAC)- тип M41L+T215Y

NType(8)=NType(8) + 1;

otherwise% (TTG ---)- тип M41L

NType(2)=NType(2) + 1;

end

else

switch pop(i, 643:645)

case'AAC'% (--- AAC)- тип T215N

NType(3)=NType(3) + 1;

case'TCC'% (--- TCC)- тип T215S

NType(4)=NType(4) + 1;

case'TAC'% (--- TAC)- тип T215Y

NType(5)=NType(5) + 1;

otherwise% тип WT

NType(1)=NType(1) + 1;

end

end

end

 

%% отборч.1

MType = zeros(1,8); %количество генов каждого типа, отбираемых в новую популяцию из данной

Mpop = 0;

while Mpop<Npop

RandSelect = rand;

if RandSelect<=Roulette03(1)

if MType(1)<NType(1)*20

MType(1)=MType(1)+1;

Mpop = Mpop + 1;

end

elseif RandSelect<=Roulette03(2)

if MType(2)<NType(2)*20

MType(2)=MType(2)+1;

Mpop = Mpop + 1;

end

elseif RandSelect<=Roulette03(3)

if MType(3)<NType(3)*20

MType(3)=MType(3)+1;

Mpop = Mpop + 1;

end

elseif RandSelect<=Roulette03(4)

if MType(4)<NType(4)*20

MType(4)=MType(4)+1;

Mpop = Mpop + 1;

end

elseif RandSelect<=Roulette03(5)

if MType(5)<NType(5)*20

MType(5)=MType(5)+1;

Mpop = Mpop + 1;

end

elseif RandSelect<=Roulette03(6)

if MType(6)<NType(6)*20

MType(6)=MType(6)+1;

Mpop = Mpop + 1;

end

elseif RandSelect<=Roulette03(7)

if MType(7)<NType(7)*20

MType(7)=MType(7)+1;

Mpop = Mpop + 1;

end

else

if MType(8)<NType(8)*20

MType(8)=MType(8)+1;

Mpop = Mpop + 1;

end

end

end

 

%% отборч.2

Nnewpop=0; %количество уже отобранных в новую популяцию цепочек генов

NnewType=zeros(1, 8); %количество генов каждого типа, отобранных в новую популяцию из данной

i = 1; %индекс для старой популяции

j = 1; %индекс для новой популяции

while Nnewpop<Npop

if strcmp(pop(i, 121:123),'TTG') %проходим

switch pop(i, 643:645)

case'AAC'

if NnewType(6)<MType(6)

newpop(j, :) = pop (i, :);

NnewType(6)= NnewType(6) +1;

Nnewpop = Nnewpop + 1;

j = j + 1;

end

case'TCC'

if NnewType(7)<MType(7)

newpop(j, :) = pop (i, :);

NnewType(7)= NnewType(7) +1;

Nnewpop = Nnewpop + 1;

j = j + 1;

end

case'TAC'

if NnewType(8)<MType(8)

newpop(j, :) = pop (i, :);

NnewType(8)= NnewType(8) +1;

Nnewpop = Nnewpop + 1;

j = j + 1;

end

otherwise

if NnewType(2)<MType(2)

newpop(j, :) = pop (i, :);

NnewType(2)= NnewType(2) +1;

Nnewpop = Nnewpop + 1;

j = j + 1;

end

end

else

switch pop(i, 643:645)

case'AAC'

if NnewType(3)<MType(3)

newpop(j, :) = pop (i, :);

NnewType(3)= NnewType(3) +1;

Nnewpop = Nnewpop + 1;

j = j + 1;

end

case'TCC'

if NnewType(4)<MType(4)

newpop(j, :) = pop (i, :);

NnewType(4)= NnewType(4) +1;

Nnewpop = Nnewpop + 1;

j = j + 1;

end

case'TAC'

if NnewType(5)<MType(5)

newpop(j, :) = pop (i, :);

NnewType(5)= NnewType(5) +1;

Nnewpop = Nnewpop + 1;

j = j + 1;

end

otherwise

if NnewType(1)<MType(1)

newpop(j, :) = pop (i, :);

NnewType(1)= NnewType(1) +1;

Nnewpop = Nnewpop + 1;

j = j + 1;

end

end

end

if i>=Npop

i = 1;

else

i = i + 1;

end

end

 

%% подготовка к следующей итерации

pop = newpop;

GlobType(k, :) = MType;

 

%% подсчет расстояния Хэмминга

 

HamDifPairs = zeros(1, Ngene);

HamTotal = 0;

for j = 1:1:Ngene

HamNums = zeros(1, 4); % количестваразныхнуклеотидов

HamPairs = zeros(1, 4); % количество одинаковых пар для каждого нуклеотида

for i = 1:1:Npop %ищем количества разных нуклеотидов в столбце j

switch pop(i,j)

case'A'

HamNums(1) = HamNums(1) + 1;

case'G'

HamNums(2) = HamNums(2) + 1;

case'T'

HamNums(3) = HamNums(3) + 1;

case'C'

HamNums(4) = HamNums(4) + 1;

end

end

for l = 1:1:4 % считаем количество одинаковых пар для каждого нуклеотида

if (HamNums(l)>=2)

HamPairs(l) = HamNums(l)*(HamNums(l)-1)/2;

else

HamPairs(l) = 0;

end

end

HamDifPairs(1, j) = Npairs - sum(HamPairs);

% получим количество разных пар нуклеотидов в столбце j

end

HamTotal = sum(HamDifPairs)/Npairs;

% среднее расстояние Хэмминга для популяции

HamGlobal(k) = HamTotal;

 

end

 

%% построение графика квазивидов

plot(GlobType(:,:))

title('КонцентрацияAZT?M = 0,03')

xlabel('Поколение')

ylabel('Количествовирусныхцепочек')

legend('WT','M41L','T215N','T215S','T215Y','M41L+T215N','M41L+T215S','M41L+T215Y')

 

%% построение графика Расстояния Хэмминга

plot(HamGlobal(:))

title('Среднее расстояние Хэмминга')

xlabel('Поколение')

ylabel('Среднее расстояние Хэмминга')

 

 

 

ПРИЛОЖЕНИЕ В