% Simulates reversible edge tile binding % B B % B B % B B % C simulations = 10; edge_tiles = 10000; corner_tiles = 500; on_rate = 1; off_rate = 100; corner_assembly_sizes = zeros(corner_tiles*2,simulations); for i=1:simulations corner_tile_left_sizes = zeros(1,corner_tiles); number_of_left_additions = zeros(1,corner_tiles); corner_tile_right_sizes = zeros(1,corner_tiles); number_of_right_additions = zeros(1,corner_tiles); edge_assemblies = ones(1,edge_tiles); number_of_edge_assemblies = edge_tiles; % edge_bond = zeros(1,edge_tiles); steps = 0; step_limit = 50000; circular_assemblies = 0; on_sum = 1; off_sum = 0; while (steps < step_limit) % Find a random possible bond if number_of_edge_assemblies > 0 off_sum = sum(corner_tile_left_sizes) + sum(corner_tile_right_sizes) ... + sum(edge_assemblies(1: number_of_edge_assemblies)) ... - sum(edge_assemblies(1:number_of_edge_assemblies) ... == 1); else off_sum = sum(corner_tile_left_sizes) + ... sum(corner_tile_right_sizes); end; on_sum = number_of_edge_assemblies*(2*corner_tiles + number_of_edge_assemblies); corner_bond_ratio = number_of_edge_assemblies *2*corner_tiles * ... on_rate/(off_sum*off_rate + on_sum*on_rate); edge_bond_ratio = number_of_edge_assemblies^2 * on_rate/(off_sum*off_rate + ... on_sum*on_rate); % off_ratio = 1 - corner_bond_ratio - edge_bond_ratio bond = random ('Uniform', 0, 1); % And a random possible other bond: lucky_other_assembly = random('Discrete Uniform', ... number_of_edge_assemblies); if bond < corner_bond_ratio % Bond is from a left corner tile: pick the assembly lucky_assembly = ceil(bond / corner_bond_ratio * corner_tiles); corner_tile_left_sizes(lucky_assembly) = ... corner_tile_left_sizes(lucky_assembly) + ... edge_assemblies(lucky_other_assembly); number_of_left_additions(lucky_assembly) = ... number_of_left_additions(lucky_assembly)+1; edge_assemblies(lucky_other_assembly) = ... edge_assemblies(number_of_edge_assemblies); number_of_edge_assemblies = number_of_edge_assemblies - 1; % edge_bond(number_of_edge_assemblies) = 1; elseif bond < 2 * corner_bond_ratio % Bond is from a right corner tile lucky_assembly = ceil((bond - corner_bond_ratio)/corner_bond_ratio ... * corner_tiles); corner_tile_right_sizes(lucky_assembly) = ... corner_tile_right_sizes(lucky_assembly) + ... edge_assemblies(lucky_other_assembly); number_of_right_additions(lucky_assembly) = ... number_of_right_additions(lucky_assembly)+1; edge_assemblies(lucky_other_assembly) = ... edge_assemblies(number_of_edge_assemblies); number_of_edge_assemblies = number_of_edge_assemblies - 1; % edge_bond(number_of_edge_assemblies) = 1; elseif bond < (2*corner_bond_ratio + edge_bond_ratio) % Bond is between two edge assemblies lucky_assembly = ceil((bond - 2 * corner_bond_ratio) / ... (1 - 2 * corner_bond_ratio) * ... number_of_edge_assemblies); if lucky_assembly == lucky_other_assembly 'Circular assembly; skipping'; circular_assemblies = circular_assemblies + 1; else edge_assemblies(lucky_assembly) = edge_assemblies(lucky_assembly) ... + edge_assemblies(lucky_other_assembly); edge_assemblies(lucky_other_assembly) = ... edge_assemblies(number_of_edge_assemblies); number_of_edge_assemblies = number_of_edge_assemblies - 1; end; else % Off reaction off_reaction = random('Discrete Uniform',off_sum); if off_reaction < sum(corner_tile_left_sizes) for x=1:corner_tiles if off_reaction <= sum(corner_tile_left_sizes(1:x)) corner_tile_left_sizes(x) = sum(corner_tile_left_sizes(1:x)) - ... off_reaction; number_of_edge_assemblies = number_of_edge_assemblies + ... 1; edge_assemblies(number_of_edge_assemblies) = ... off_reaction - sum(corner_tile_left_sizes(1:(x-1))); break; end; end; elseif off_reaction < (sum(corner_tile_left_sizes)+ ... sum(corner_tile_right_sizes)) off_reaction = off_reaction - sum(corner_tile_left_sizes); for x=1:corner_tiles if off_reaction <= sum(corner_tile_right_sizes(1:x)) corner_tile_right_sizes(x) = sum(corner_tile_right_sizes(1:x)) - ... off_reaction; number_of_edge_assemblies = number_of_edge_assemblies + ... 1; edge_assemblies(number_of_edge_assemblies) = ... off_reaction - sum(corner_tile_right_sizes(1:(x-1))); break; end; end; else off_reaction = off_reaction - sum(corner_tile_left_sizes) ... - sum(corner_tile_right_sizes); for x=1:number_of_edge_assemblies if off_reaction <= (sum(edge_assemblies(1:x)) - ... sum(edge_assemblies(1:x) == 1)) edge_assemblies(x) = sum(edge_assemblies(1:x)) - off_reaction; number_of_edge_assemblies = number_of_edge_assemblies + 1; edge_assemblies(number_of_edge_assemblies) = ... off_reaction - sum(edge_assemblies(1:(x-1))); break; end; end; end; end; steps = steps+1; if (mod(steps,500) == 0) steps end; end; corner_assembly_sizes(:,i) = [corner_tile_left_sizes, ... corner_tile_right_sizes(1:corner_tiles)]'; end; % Display the current assemblies %max_corner_size = max(max(corner_assembly_sizes)); max_corner_size = 60; min_corner_size = 0; bins = 21; size(histc(corner_assembly_sizes(:,1),linspace(min_corner_size, ... max_corner_size,bins))) for i=1:simulations corner_assembly_sizes_hist(:,i) = histc(corner_assembly_sizes(:,i),linspace(min_corner_size, ... max_corner_size,bins)); end; corner_assembly_sizes_mean = mean(corner_assembly_sizes_hist,2); corner_assembly_sizes_std = std(corner_assembly_sizes_hist,0,2); figure(5); hold off; bar(linspace(min_corner_size,max_corner_size,bins),corner_assembly_sizes_mean,'y'); hold on; errorbar(linspace(min_corner_size,max_corner_size,bins), ... corner_assembly_sizes_mean,corner_assembly_sizes_std,'.'); axis([-2,62,0,800]); xlabel('Corner Assembly Length, in Number of Tiles','FontSize',20); ylabel('Number of Assemblies','FontSize',20); title('Reversible Aggregation 10:1','FontSize',20); set(gca,'FontSize',20); hold off; saveas(gcf, 'reversible_aggregation_length_distribution', 'fig'); print -depsc2 reversible_aggregation_length_distribution.eps;