% % boadraw08.m Tony Kellems and Steve Cox % % Draw a basin of attraction of a gene net that looks more % like Wuensche's (radially outward). % % Usage: boadraw08(s,wire,rule) % % where s = initial state vector [1 by n] % wire = wiring table matrix [n by m] % rule = rule vector [n by 1] % % s = [1 0 0 1 0 1]; % wire = [4 2 3;5 3 2;3 6 1;5 4 6;6 1 2;3 5 6] % rule = [231 64 5 108 61 62]'; % % or n=6; boadraw08(round(rand(1,n)),gring(n),round(255*rand(n,1))); % if you have saved gring to separate file % % helper functions: d2b, getatt, pre, bpre, gring % function boadraw08(s,wire,rule) color = {'b' 'r' 'g' 'k' 'm' 'c'}; disp('Getting attractor...') att = getatt(s,wire,rule); % get the attractor [latt,n] = size(att); ang = 0:2*pi/latt:2*pi; % plot the attractor disp(['Found ' num2str(latt) ' states in attractor']) plot(cos(ang),sin(ang),'linewidth',3) hold on plot(cos(ang),sin(ang),'yo') for k=1:n, ruletab(k,:) = d2b(rule(k),8); % build the rule table end % prepare for brute force pre func s1 = zeros(2^n,n); s2 = s1; for j=0:2^n-1 s1(j+1,:) = d2b(j,n); % each possible state s2(j+1,:) = suc(s1(j+1,:),wire,rule); % its successor end for j=1:latt, cnum = mod(j,6); if cnum == 0 cnum = 6; end clab = char(color(cnum)); disp(['Computing predecessors for state ' num2str(j) '...']) p = bpre(att(j,:),s1,s2); np = size(p,1); % strip attractor preimage [extra, attind, pind] = intersect(att,p,'rows'); p = p([1:pind-1 pind+1:end],:); np = size(p,1); fang = linspace(ang(j)-pi/8,ang(j)+pi/8,np + 1); % fan of angles for k=1:np % plot the preimages XX(k,:) = [cos(ang(j)) 1*cos(ang(j))+cos(fang(k))]; YY(k,:) = [sin(ang(j)) 1*sin(ang(j))+sin(fang(k))]; plot(XX(k,:),YY(k,:),clab) plot(XX(k,:),YY(k,:),'yo') end if np > 0, X = XX; Y = YY; end c1 = 1; % start of unprocessed preimages c2 = np; % end of unprocessed preimages gen = 1; % generation counter for plotting while c2 >= c1 % compute the preimages of the preimages gen = gen + 1; disp(['Computing preimages for generation ' num2str(gen) '...']) cnt = c2; pc = 0; % preimage counter fang = linspace(ang(j)-pi/8,ang(j)+pi/8,c2 - c1 + 1); % fan of angles for k=c1:c2 q = bpre(p(k,:),s1,s2); cq = size(q,1); if cq > 0, p(cnt+1:cnt+cq,:) = q; cnt = cnt + cq; fxang = linspace(fang(k-c1+1)-pi/8,fang(k-c1+1)+pi/8,cq); % new fan of angles for next set of preimages for kk=1:cq % plot the new preimages x(pc+kk,:) = [X(k-c1+1,2) gen*cos(ang(j))+cos(fxang(kk))]; y(pc+kk,:) = [Y(k-c1+1,2) gen*sin(ang(j))+sin(fxang(kk))]; plot(x(pc+kk,:),y(pc+kk,:),clab) plot(x(pc+kk,:),y(pc+kk,:),'yo') end % kk pc = pc + cq; end % if end % k c1 = c2 +1; c2 = cnt; if pc > 0, X = x; Y = y; clear x y end end % while clear p end % j disp('Finished') ax = axis; xt = ax(2)+1; yt = ax(4); ys = ax(4)-ax(3); for j=1:n text(xt,yt-j*ys/(n+1),[num2str(wire(j,:)) ' ' num2str(rule(j))]) end hold off axis(ax+[0 3 0 0]) axis equal title(['seed = ' num2str(s)]) axis off % % getatt.m % % follow a given state into its attractor % % usage: att = getatt(s,wire,rule) % % where s is the initial state % wire is the wire diagram % rule is the vector of gene rules % and % att is the attractor % function att = getatt(s,wire,rule); n = length(s); iter = 1; % state counter mi = 0; % match place while mi == 0, ns = suc(s(iter,:),wire,rule); % get successor of s [jnk,mi] = ismember(ns,s,'rows'); % mi = row of s that matches ns iter = iter + 1; s(iter,:) = ns; % add ns to stack end att = s(mi:iter-1,:); % strip attractor from stack % % suc.m % % find the immediate successor of a given state % % usage ns = suc(s,wire,rule) % % where s is the given state % wire is the wiring diagram % rule is the rule vector % and % ns is the successor of s % % functions called : d2b % function ns = suc(s,wire,rule) n = length(s); for k=1:n, ruletab(k,:) = d2b(rule(k),8); % build the rule table end for k=1:n, pat = s(wire(k,:)); % input pattern to gene k cc = 1 + pat*[4 2 1]'; % base 10 of pat ns(k) = ruletab(k,cc); % new state of gene k end % % perm = pre(s,wire,ruletab) % % find all immediate predecessors of s % % this is the EXTRA CREDIT nonbruteforce method % function perm = pre(s,wire,ruletab) n = length(s); itab = [0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 1]; ind = cell(1,n); % holding cell for permissible precursor % states of each gene for j=1:n ind{j} = find(s(j)==ruletab(j,:)); end perm = []; % initialize full precursor list for i = 1:length(ind{1}), ts = NaN(1,n); ts(wire(1,:)) = itab(ind{1}(i),:); % ith trial candidate from s(1) list p = ts; % initialize trial at gene 1 pc = 1; % trial counter for j=2:n % step through genes 2 through n and turn trial members % into predecessors when appropriate p2 = []; % true precursors consistent with ind{1}(i) pc2 = 0; % counter for k=1:length(ind{j}), % check all possibilities at gene j ts2 = NaN(1,n); ts2(wire(j,:)) = itab(ind{j}(k),:); % gene j trial for m=1:pc % check for conflicts in p-list act = 1 - isnan(p(m,:)); % non-NaN indices in p(m,:) tact = intersect(find(act),wire(j,:)); % common indices if ts2(tact) == p(m,tact) tmp = p(m,:); tmp(wire(j,:)) = ts2(wire(j,:)); pc2 = pc2 + 1; p2 = [p2; tmp]; % add candidate to p2-list end end % m end % k p = p2; pc = pc2; end % j perm = [perm; p2]; % append p2 to full precursor list end % i % % d2b.m % % convert base 10 to base 2 % % usage bin = d2b(dec,n) % % where dec is the base 10 number % n is the number of bianry digits % and % bin is the binary translation of dec % function bin = d2b(dec,n) bin = zeros(1,n); % initialize if dec == 0 % if dec = 0 return return end h = floor(log2(dec)); % build sequential powers of two until no remainder bin(h+1) = 1; dec = dec - 2^h; while dec > 0 h = floor(log2(dec)); bin(h+1) = 1; dec = dec - 2^h; end % % optional, this is here for convenience, % create wire for a ring of n genes. % function wire = gring(n) wire = zeros(n,3); wire(1,:) = [n 1 2]; for j=2:n-1, wire(j,:) = [j-1 j j+1]; end wire(n,:) = [n-1 n 1]; % % the brute force pre function % function perm = bpre(s,s1,s2) match = find(ismember(s2,s,'rows')); perm = s1(match,:); if isempty(perm) return end m = ismember(perm,s,'rows'); if m>0 perm = perm([1:m-1 m+1:end],:); end