# fitness model # assign fitness according to a table # define a mutation matrix # 32 bit genome # directions: # click on InitialCond Range. Then click OK. After a short bit, click Restore # shows x0,x15,x23,x30,x31 as function of u # see below for more infor par u=0 # some bit manipulation functions # shift 1 sft(x)=flr(x/2) # shift n sftn(x,n)=if(n<1)then(x)else(sftn(sft(x),n-1)) # nth bit of x bit(x,n)=mod(sftn(x,n),2) # hamming distance up to n bits hd(x,y,n)=sum(0,n-1)of(abs(bit(x,i')-bit(y,i'))) # here is the HW table table f fbook.tab # this makes a table of the hamming distances between two sequences table h % 1024 0 1023 hd(flr(t/32),mod(t,32),5) # the matrix q table q % 1024 0 1023 (1-u)^(5-h(t))*u^h(t) # fx[0..31]=f([j])*x[j] # average fitness phi=sum(0,31)of(shift(fx0,i')) # matrix multiplication special qf=mmult(32,32,q,fx0) x[0..31]'=qf([j])-phi*x[j] # everyone starts out equal! init x[0..31]=.0325 aux fbar=phi aux mut=u aux xtot=sum(0,31)of(shift(x0,i')) # control stuff # integrate to 100 but just keep last value - probably steady state @ total=100,dt=.1,meth=qualrk,tol=1e-5,atol=1e-8 @ trans=99.99 # plot stuff 5 plots @ nplot=5,xp=mut,xlo=0,xhi=.5,ylo=0 @ xp5=mut,xp2=mut,xp3=mut,xp4=mut @ yp=x0,yp5=x15,yp2=x23,yp3=x30,yp4=x31 # range stuff @ rangeover=u @ rangelow=0,rangehigh=0.5,rangestep=50,rangereset=no done This file solves the quasi species equations it is set up now to only look at the steady states since all transients are integrated away. plots are xi vs u - the mutation you can change the fitness profile by altering fbook.tab this is a list of fitnsses for 0-31 (first line starts with 10)