classdef WriFish < matlab.mixin.SetGet & dynamicprops % The software models the change in frequency of the B allele subject to initial f(B) and genotype fitness rules. %----------------------------------------------- % Author: Adam Gogacz $ Steven M. Carr % Version: 1.0 % Last Revision: 2020-10-30 % Copyright: Adam Gogacz %----------------------------------------------- %% --- parameters and constants properties (SetAccess=private,SetObservable=true,Hidden=true) %---- modelling parameters and variables q_init =0.5; % f(a) = q; w0 =1; % fitness of AA genotype; must be between 0 and Inf w1 =1; % fitness of AB genotype; must be between 0 and Inf w2 =1; % fitness of BB genotype; must be between 0 and Inf Nindiv =100; % number of individuals in the population Ngen =200; % number of generations Npop =100; % number of populations N_BB double; N_AB double; qT double; %--- GUI controls axes matlab.ui.control.UIAxes; end %--- GUI constants properties (Constant=true,Hidden=true) uiHt=20; uiGap=10; uiHelp=15; Name='Wright-Fisher Model'; end methods function this=WriFish(varargin) matlabScreenrSize=get(0,'screensize'); scW=matlabScreenrSize(1,3); scH=matlabScreenrSize(1,4); figTag='wrightfisher'; delete(findall(0,'Type','figure','-and','Tag',figTag)); hf=uifigure('Name',this.Name,'Tag',figTag,'Position',[1,0.25*scH,scW,0.75*scH-36],'Visible',0); mainGridLayout=matlab.ui.container.GridLayout('Parent',hf,'RowHeight',{'1x','1x'},'ColumnWidth',{220,'1x','1x'},... 'RowSpacing',48,'ColumnSpacing',10,'Scrollable','off'); this.axes=matlab.ui.control.UIAxes('Parent',mainGridLayout); this.axes.Layout.Row=[1,2]; this.axes.Layout.Column=[2,3]; %--- labels this.axes.XAxis.Label.String='Generations'; this.axes.XAxis.Label.FontWeight='bold'; this.axes.YAxis.Label.String = 'f(B)'; this.axes.YAxis.Label.FontWeight='bold'; this.axes.XLim=[1,this.Ngen]; this.axes.YLim=[-0.01,1.01]; this.axes.XGrid='on'; this.axes.YGrid='on'; % --- add input controls glTP=matlab.ui.container.GridLayout('Parent',mainGridLayout,'RowSpacing',6,'ColumnSpacing',6,... 'ColumnWidth',{'2x','1x'},'Scrollable','on','Padding',[0,0,0,0],... 'RowHeight',{... this.uiHt,... % EVOLUTION CONTROLS header this.uiHt,this.uiHt,this.uiHt,... % w0, w1, and w2 uicontrols this.uiHt,... % initial f(B) frequency this.uiHt,this.uiHt,this.uiHt... % population size, number of generations, number of populations this.uiHt,... % spacer this.uiHt,... % STATS header 2.7*this.uiHt,... % stats table '1x',... 1.5*this.uiHt,... % RUN button });glTP.Layout.Row=[1,2];glTP.Layout.Column=1; % --- modelling parameters ir=1; hdummy=matlab.ui.container.Panel('Parent',glTP,'BorderType','none','Title','Model Parameters',... 'FontWeight','bold','TitlePosition','centertop','ForegroundColor','#0072BD'); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=[1,2]; hdummy=matlab.ui.control.Label('Parent',glTP,'Text','AA fitness W0','FontWeight','bold',... 'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1; hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%g','Value',this.w0,'Limits',[0,Inf],'UpperLimitInclusive',0,... 'ValueChangedFcn',@(h,~)set(this,'W0',h.Value)); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2; hdummy=matlab.ui.control.Label('Parent',glTP,'Text','AB fitness W1','FontWeight','bold',... 'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1; hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%g','Value',this.w1,'Limits',[0,Inf],'UpperLimitInclusive',0,... 'ValueChangedFcn',@(h,~)set(this,'W1',h.Value)); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2; hdummy=matlab.ui.control.Label('Parent',glTP,'Text','BB fitness W2','FontWeight','bold',... 'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1; hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%g','Value',this.w2,'Limits',[0,Inf],'UpperLimitInclusive',0,... 'ValueChangedFcn',@(h,~)set(this,'W2',h.Value)); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2; hdummy=matlab.ui.control.Label('Parent',glTP,'Text','Initial q=f(B)','FontWeight','bold',... 'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1; hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%g','Value',this.q_init,'Limits',[0,Inf],'UpperLimitInclusive',0,... 'ValueChangedFcn',@(h,~)set(this,'q_init',h.Value)); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2; hdummy=matlab.ui.control.Label('Parent',glTP,'Text','Population size','FontWeight','bold',... 'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1; hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%d','Value',this.Nindiv,'Limits',[1,Inf],'LowerLimitInclusive',0,... 'UpperLimitInclusive',0,'ValueChangedFcn',@(h,~)set(this,'Nindiv',h.Value),'RoundFractionalValues',1); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2; hdummy=matlab.ui.control.Label('Parent',glTP,'Text','No. generations','FontWeight','bold',... 'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1; hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%d','Value',this.Ngen,'Limits',[1,Inf],'LowerLimitInclusive',0,... 'UpperLimitInclusive',0,'ValueChangedFcn',@(h,~)set(this,'Ngen',h.Value),'RoundFractionalValues',1); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2; hdummy=matlab.ui.control.Label('Parent',glTP,'Text','No. populations','FontWeight','bold',... 'HorizontalAlignment','left','VerticalAlignment','bottom'); hdummy.Layout.Row=ir;hdummy.Layout.Column=1; hdummy=matlab.ui.control.NumericEditField('Parent',glTP,'ValueDisplayFormat','%d','Value',this.Npop,'Limits',[1,Inf],'LowerLimitInclusive',0,... 'UpperLimitInclusive',0,'ValueChangedFcn',@(h,~)set(this,'Npop',h.Value),'RoundFractionalValues',1); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=2; % --- output stats ir=ir+1; hdummy=matlab.ui.container.Panel('Parent',glTP,'BorderType','none','Title','Results',... 'FontWeight','bold','TitlePosition','centertop','ForegroundColor','#0072BD'); hdummy.Layout.Row=ir;ir=ir+1; hdummy.Layout.Column=[1,2]; htable=matlab.ui.control.Table('Parent',glTP,'ColumnEditable',[false,false],'ColumnFormat',{'char','char'},'ColumnWidth',{'2x','1x'},'Data',cell.empty,... 'RowName',{},'ColumnName',{}); addStyle(htable,uistyle('FontWeight','bold'),'column',1); htable.Layout.Row=ir; htable.Layout.Column=[1,2]; hdummy=matlab.ui.control.Button('Parent',glTP,'Text','Run Evolution','FontWeight','bold','BackgroundColor','#EDB120',... 'FontSize',this.uiHt,'Interruptible',false,'ButtonPushedFcn',@(~,~)run_buttonpushedfcn()); hdummy.Layout.Row=numel(glTP.RowHeight);hdummy.Layout.Column=[1,2]; function run_buttonpushedfcn() try run(); plotData(); catch ex hdummy=uiconfirm(hf,{'CRITICAL ERROR,';' ';['Identifier: ',ex.identifier];['Message: ',ex.message]},... this.Name,'Options',{'OK'},'Icon','error'); end end %--- turn figure visibility on hf.Visible=1; %% --- supporting functions function run() try hbar=matlab.ui.dialog.ProgressDialog(hf,'Indeterminate',0,'ShowPercentage',1,'Title',this.Name,'Message','Evolution in progress...'); wmax=max([this.w0,this.w1,this.w2]); if(0==wmax); uialert(hf,'At least one fitness value must be non-zero.',this.Name); return; end %--- calculated values this.N_BB = zeros(this.Npop,this.Ngen); this.N_AB = zeros(this.Npop,this.Ngen); w0loc=this.w0/wmax; w1loc=this.w1/wmax; w2loc=this.w2/wmax; Nindiv=this.Nindiv; Ngen=this.Ngen; Npop=this.Npop; for s = 1:Npop hbar.Value=s/Npop; q = this.q_init; %--- initialize variables this.N_BB(s,1) = round(Nindiv*q^2); this.N_AB(s,1) = round(Nindiv*2*q*(1-q)); %--- proceed with looping over generations for m = 2:Ngen %--- reset variables g=1; ct=1; genGenotype=nan(1,Nindiv); %--- begin looping until maximum capacity is reached (Npop) while(g <= Nindiv) %--- generate a genotype donAlpha = rand(1) >= q; % first allele; 1 represents presence of "A" allele, and 0 for "B" allele donBeta = rand(1) >= q; % second allele; 1 represents presence of "A" allele, and 0 for "B" allele genotype = donAlpha+donBeta; % genotype of the individual before fitness application if(2==genotype) % AA genotype if(1==w0loc || w0loc >= rand(1)); genGenotype(g)=2; g=g+1; end elseif(1==genotype) % AB genotype if(1==w1loc || w1loc >= rand(1)); genGenotype(g)=1; g=g+1; end else % BB genotype if(1==w2loc || w2loc >= rand(1)); genGenotype(g)=0; g=g+1; end end %--- update loop counter ct=ct+1; end this.N_BB(s,m)=nnz(0==genGenotype); this.N_AB(s,m)=nnz(1==genGenotype); %--- update q q = (this.N_BB(s,m) + this.N_AB(s,m)/2)/Nindiv; % frequency of "B" allele %--- test if q went to 0 or 1, if so, break out of the loop and update the N_AB and N_BB matrices if(0==q || 1==q) this.N_BB(s,m+1:end)=this.N_BB(s,m); this.N_AB(s,m+1:end)=this.N_AB(s,m); break; end end end %--- update variables qT=transpose((this.N_BB + this.N_AB/2)/Nindiv); %--- update table htable.Data={... 'f(B)=1 fixation',nnz(1==max(qT,[],1));... 'f(B)=0 loss',nnz(0==min(qT,[],1))}; %--- clean up delete(hbar); catch ex delete(hbar); rethrow(ex); end end function plotData() g = 1:this.Ngen; plot(this.axes,g,transpose((this.N_BB + this.N_AB/2)/this.Nindiv)); this.axes.XLim=[1,this.Ngen]; end end end %--- ICONS properties(Constant=true,Hidden=true) helpIcon=[255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255; 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255; 255,255,255,255,255,255,255,239,180,146,121,155,154,129,146,180,239,255,255,255,255,255,255,255; 255,255,255,255,255,255,178,126,186,238,255,254,250,246,227,186,132,178,255,255,255,255,255,255; 255,255,255,255,247,140,158,246,255,255,255,254,250,246,242,239,227,158,140,247,255,255,255,255; 255,255,255,255,137,164,255,255,255,255,242,205,202,234,242,238,234,230,183,137,255,255,255,255; 255,255,255,171,143,255,255,255,255,155,57,60,64,78,177,237,233,230,226,148,171,255,255,255; 255,255,236,102,235,255,255,255,155,57,60,87,90,70,74,186,232,229,225,208,108,236,255,255; 255,255,167,154,251,253,254,254,70,60,144,247,244,126,79,121,231,228,222,211,143,167,255,255; 255,255,124,199,248,249,250,250,107,109,223,244,242,190,85,125,229,226,218,206,176,124,255,255; 255,255,89,214,244,245,246,246,246,245,243,241,239,124,90,145,227,222,212,201,189,89,255,255; 255,255,86,239,241,242,242,242,242,241,240,120,104,90,104,203,224,216,205,195,182,86,255,255; 255,255,81,235,237,238,238,238,238,237,236,85,90,153,211,224,218,208,198,186,175,81,255,255; 255,255,78,193,233,234,234,234,234,233,232,90,96,227,224,218,209,200,190,179,162,78,255,255; 255,255,108,161,229,230,231,231,230,230,229,127,132,222,216,208,200,191,180,170,133,108,255,255; 255,255,152,100,225,226,227,227,227,226,224,163,165,200,205,198,190,180,171,160,95,152,255,255; 255,255,233,69,171,219,220,220,219,217,214,107,112,178,194,186,179,170,159,126,69,233,255,255; 255,255,255,149,75,189,207,207,207,205,202,132,135,175,182,175,168,157,137,79,149,255,255,255; 255,255,255,255,101,89,178,194,194,192,189,185,181,176,170,163,151,132,76,101,255,255,255,255; 255,255,255,255,244,98,69,137,179,178,175,172,168,163,154,142,110,69,98,244,255,255,255,255; 255,255,255,255,255,255,145,59,73,106,137,134,129,121,97,68,59,145,255,255,255,255,255,255; 255,255,255,255,255,255,255,231,144,94,57,57,57,57,94,144,231,255,255,255,255,255,255,255; 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255; 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255] end end