function [err, FId, FCd] = SAkin29proofM(Gmc,Gse,rs,N,k) % [err, FId, FCd] = SAkin29proofM(Gmc,Gse,rs,N,k) % 29 state kinetic trap model for 2x2 proofreading tiles % forward rate f = k*exp(-Gmc) k = 20*10^6 % reverse rate r_b = k*exp(-b*Gse) % freezing rate rs for complete 2x2 blocks % assumes a BCA tileset with N+1 symbols % % error rate is the number of bad BLOCKS. % the per/tile error rate is 2x lower. % % uses matrix inversion to do the calculation, rather than algebra. f = k*exp(-Gmc); r4 = k*exp(-4*Gse); r3 = k*exp(-3*Gse); r2 = k*exp(-2*Gse); r1 = k*exp(-Gse); % state diagram is awful -- see my notebook. M = zeros(29); % kf M(j,i) += kf; M(i,i) -= kf; % i <---> j ==> % kr M(i,j) += kr; M(j,j) -= kr; % i j kf kr reacts = [ 1 2 2*N*f r1 ; ... 1 3 2*N*f r1 ; ... 1 4 f r2 ; ... 1 5 2*f r1 ; ... 2 6 f r1 ; ... 2 7 f r2 ; ... 3 7 f r2 ; ... 3 8 f r1 ; ... 3 9 f r2 ; ... 4 9 2*N*f r1 ; ... 4 10 2*f r2 ; ... 5 10 f r3 ; ... 5 11 f r1 ; ... 6 12 f r2 ; ... 6 18 f r1 ; ... 7 12 f r1 ; ... 7 13 f r1 ; ... 7 19 N*f r1 ; ... 7 20 f r1 ; ... 8 14 f r1 ; ... 8 20 f r2 ; ... 8 21 f r2 ; ... 9 15 f r2 ; ... 9 21 f r1 ; ... 10 15 N*f r1 ; ... 10 16 f 2*r2; ... 10 22 f r1 ; ... 11 17 f 2*r2; ... 11 22 f r3 ; ... 12 23 f r2 ; ... 13 25 f r1 ; ... 14 23 f r3 ; ... 14 26 f r2 ; ... 15 27 f r1 ; ... 16 29 f r2 ; ... 17 28 2*N*f r1 ; ... 17 29 f r4 ; ... 18 23 f r3 ; ... 19 24 f r1 ; ... 20 23 f r2 ; ... 21 26 f r1 ; ... 21 27 f r2 ; ... 22 29 f 2*r3; ... ]; for n=1:size(reacts,1) i=reacts(n,1); j=reacts(n,2); kf=reacts(n,3); kr=reacts(n,4); M(j,i) = M(j,i)+kf; M(i,i) = M(i,i)-kf; M(i,j) = M(i,j)+kr; M(j,j) = M(j,j)-kr; end %%%%% freezing reactions % innermost states only -- right shape, but error rate 2x too low badguys = [27]; goodguys = [29]; % completely wrong state -- matches SAkin9proofM, but not data badguys = [23]; goodguys = [29]; % all states with both outer tiles on both sides -- best agreement with data badguys = [14 23 24 25 26 27]; goodguys = [17 28 29]; for i=[goodguys badguys] M(i,i)=M(i,i)-rs; end xss = inv(M) * [-1 zeros(1,28)]'; % steady-state, dumping 1 unit/sec into E FId = rs*sum(xss(badguys)); FCd = rs*sum(xss(goodguys)); err = FId/(FId+FCd); % this should be exactly equal to FId