function out=FAMonomers5d(kox, kelong, bias, mut, draw)
 
% Will simulate replication (=fibers elongation-breakage) due to food bound at fibers, and how it changes with the cofactor.
% Species: Monomers=A, B; Hexamers=A6, A3B3, and B6 (assumed all are assembled); Trimers=A3, A1.5B1.5, and B3 (each can be bound to the different fibers (hexamers)).
% Oxidation is mon->bound trimers (kox parameter). To represent the selectivity of the cofactor, kox can be biased when near A6, A3B3, or B6 hexamers, using the 'bias' parameter.
% Elongation is the concerted step of converting trimers into assembled hexamers. Nucleation of new fibers is omitted.
% Mutations: how likely it is for food-attached to A6 to give rise to A3B3 and B6 (and likewise for the other hexamers)
% The default parameters giving rise to Figure S123 are: kox=1e-1, kelong=1e3, bias=[2 1 1], and mut=[0.9 0.1 0; 0.2 0.6 0.2; 0 0.1 0.9].

% Feel free to contact if you have any questions or comments: Omer Markovitch https://sites.google.com/site/omermar/

if exist('draw', 'var');
	if ~islogical(draw); draw=strcmpi(draw, 'draw'); end;
else
	draw=true;
end;

secs=tic;
out=[];
out.function=mfilename('fullpath');
out.machine=computer;
out.Date=datestr(now);
out.license=license;
out.matlabversion=version;

trange=[0:0.1:100]; %hours
out.c0=[zeros(1,2+3*3) 0.2/6/10 8*0.2/6/10 0.2/6/10]; %Initial concentrations in mM: A(free), B(free), {A3, mixed, B3}@A6, {A3, mixed, B3}@A3B3, {A3, mixed, B3}@B6, A6, A3B3, B6
out.nbb=[1 1 3*ones(1,9) 6*ones(1,3)]; %how many BB are in each species
out.inflow=[0.3 0.3]*10/500; %mM/hour
out.outflow=20/500; %1/hour
if ~exist('kox', 'var') || isempty(kox); out.kox=0.1; else out.kox=kox; end;
if ~exist('kelong', 'var') || isempty(kelong); out.kelong=1e3; else out.kelong=kelong; end;
if ~exist('bias', 'var') || isempty(bias); out.bias=[2 1 1]; else out.bias=bias; end; %how much is the relative bias of food ataching to fibers of certein type (unitless)
if ~exist('mut', 'var') || isempty(mut); out.mut=[0.9 0.1 0; 0.2 0.6 0.2; 0 0.1 0.9]; else out.mut=mut; end; %mut(i,j) is the probability that food-attached-to-hex_i will give rise to h_j

[out.fitN, out.trange]=myfun5c(out.kox, out.kelong, out.bias, out.mut, out.inflow, out.outflow, out.c0, trange);
out.inBB=out.fitN.*out.nbb;
out.inBB=out.inBB./sum(out.inBB,2);

if draw;
	colors={'green' 'cyan' 'black' 'blue' 'red'};
	
	close all;	
	figure;
	hold on;
	idx=[1:2 length(out.nbb)-2:length(out.nbb)]; %monomers and hexamers
	for i=1:length(idx); plot(out.trange, out.inBB(:,idx(i)), 'LineWidth',3, 'Color',colors{i}); end;
	xlabel('Hours');
	ylabel('Fraction of BB');
	box on;
	set(gca, 'LineWidth',1, 'FontName', 'Verdana','FontSize',12, 'TickLength',[0.02 0], 'TickDir', 'both');
	legend({'A' 'B' 'A_6' 'A_3B_3' 'B_6'}, 'Location', 'best');
	hold off;
end;

beep;
out.seconds=toc(secs);
return;

function [fitN, t]=myfun5c(kox, kass, bias, mut, inflow, outflow, c0, trange)
	odeopt=odeset('NonNegative', [1:numel(c0)]);
	r13=[1 1 3 1; %A+A->A3(trimer), factor=1;
		1 2 4 3;
		2 2 5 1];
	rr=[3 3 12; %A3+A3->A6 (index=12)
		3 5 13;
		4 4 13;
		5 5 14];
	r36=[rr mut(1,[1 2 2 3])'; rr mut(2,[1 2 2 3])'; rr mut(3,[1 2 2 3])'];
	
	[t, fitN]=ode15s(@odesimple5d, trange, c0, odeopt, kox*bias, kass, inflow, outflow, r13, r36);
return;

function dcdt=odesimple5d(t, c, kox, kass, inflow, outflow, r13, r36)
dcdt=zeros(size(c));

dcdt(1)=dcdt(1)+inflow(1)-outflow*c(1); %A monomer
dcdt(2)=dcdt(2)+inflow(2)-outflow*c(2); %B monomer
dcdt(end-2:end)=dcdt(end-2:end)-c(end-2:end)*outflow; %hexamers
dcdt(3:end-3)=dcdt(3:end-3)-c(3:end-3)*outflow; %trimers

hass=c(end-2:end)/sum(c(end-2:end));

for h=1:3;
	
	k=kox(h)*(1+hass(h)); %oxidation is catalyzed by the total amount of hexamers
	for i=1:size(r13,1);
		r1=r13(i,1);
		r2=r13(i,2);
		p1=r13(i,3)+(h-1)*3;
		r=k*sqrt(c(r1)*c(r2))*r13(i,4); %order 1 in monomers
		dcdt(r1)=dcdt(r1)-1.5*r;
		dcdt(r2)=dcdt(r2)-1.5*r;
		dcdt(p1)=dcdt(p1)+r;
	end; %oxidation

end;

for i=1:size(r36,1);
	h=ceil(i/4);
	k=kass*hass(h);
	r1=r36(i,1)+3*(h-1); %the 'food' that is assembled on fib/hex 'h'
	r2=r36(i,2)+3*(h-1); %the 'food' that is assembled on fib/hex 'h'
	p1=r36(i,3);
	r=k*(c(r1)*c(r2))*r36(i,4);
	dcdt(r1)=dcdt(r1)-r;
	dcdt(r2)=dcdt(r2)-r;
	dcdt(p1)=dcdt(p1)+r;
end; %elongation

return;

