# 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)