% composomes.m (commented version of find_envi_22_6_dan.m) 
% Program used for generating the data in
% Segre' D., Ben-Eli D. and Lancet D. 
% Proc. Natl. Acad. Sci. USA (2000), 97(8), 4112-4117: 
% "Compositional genomes: prebiotic information transfer in mutually
% catalytic noncovalent assemblies."
% See http://ool.weizmann.ac.il
% for further info
%
% Written by Daniel Segre' and Dafna Ben-Eli, June 1999
% at the dept. of Molecular Genetics of the Weizmann Institute, Israel.
% For questions, please contact daniel.segre@weizmann.ac.il
%
% The program is written for Matlab 5
%
% This is a matlab version of fluxes3, with the following assumptions:
% rho is constant ( very large food suply), and equal to 1/NG.
% Constant Population

%----------------------------------------------------------------------------------

clear;

cycles=4000;       % n. of time steps
NG=100;            % n. of different molecular types 
Kf=10^-2;          % forward (joining) rate
Kb=10^-5;          % backward rate
rho=1/NG;          % concentration in solution
Nmin=40;           % initial size of assemblies 

do_split=1         % if 1 do splitting
limited_food=1     % if 1 use limited food supply
n_TOT=1000;        % total available number of monomers of each kind
W=1;               % number of assemblies (this version works only with W=1)
G_split=2;         % ratio to initial volume for splitting
                   % for normal splitting put =2
dt = 0.05;         % time scale
Poisson_threshold=20; % threshold below which reactions are stochastic

% The beta matrix
% element row i column j: species j catalyses species i by that.

%load Beta_dan_100_22_6_99; 
load Beta_mu5_23_6_99

%----------------------------------------------------------------------------------

environment=zeros(cycles,NG);   % initializes compositions of assemblies
concentration=zeros(cycles,NG); % initializes concentrations

%----------------------------------------------------------------------------------
% CREATE INITIAL CONDITIONS:
% Note the two parts below are two possible ways of entering initial conditions.
% The unwanted choice should be commented.

% OPTION 1
% Create initial random assemblies; matrix n is NG*W; sum of each column Nmin.

n_rand_index = fix(rand(Nmin,W)*NG) + 1;
n_col_mat = zeros(NG,W);

for i=1:Nmin,
   for j=1:W,
      n_col_mat(n_rand_index(i,j),j) =  n_col_mat(n_rand_index(i,j),j) + 1; 
   end;
end;


% OPTION 2
% READS INITIAL CONDITION from FILE
load initial_condition_22_6_stepwise_growth
n_col_mat=initial_condition';

%----------------------------------------------------------------------------------
% Here the growth - splitting loop starts

% The kinetics equation: [Kf*rho*(sum(n(j)))-Kbn(i)] * (1+sum(beta(ij)n(j) / sum(n(j)))
for loop=1:cycles,
   sum_n = sum(n_col_mat);
   if (limited_food==1)
     rho=(n_TOT-n_col_mat)/(NG*n_TOT); % limited food supply
   end
   sum_n_rows = sum_n(ones(size(n_col_mat,1),1),:);
   delta_n = ((Kf*rho.*(sum_n_rows) - Kb*n_col_mat) .* ( 1 + (Beta*n_col_mat)./sum_n_rows)) * dt;

   poissonize=find(delta_n<Poisson_threshold);
   disc_delta_n=round(delta_n);
   disc_delta_n(poissonize) = Poiss(delta_n(poissonize));
   n_col_mat = n_col_mat + disc_delta_n;
   if (limited_food==1)
      n_col_mat(find(n_col_mat>n_TOT))=n_TOT;
   end

   if do_split==1  % do this part when we want splitting
     sum_n = sum(n_col_mat);
     n_to_split = find (sum_n >= G_split*Nmin);
     if n_to_split
        % calculate the progeny using binom distribution
        progeny1 = binom(n_col_mat(:,n_to_split),0.5);
        %progeny2 = n_col_mat(:,n_to_split) - progeny1; 
	
        % change the matrix of live compositions	
        n_col_mat(:,n_to_split) = progeny1;
        %dead_n = fix(rand(length(n_to_split),1)*W) + 1;
        %n_col_mat(:,dead_n) = progeny2;
		
     end;
   end; % end if do split
   environment(loop,:) = n_col_mat';
   concentration(loop,:) = n_col_mat' ./ sum(n_col_mat);

end;

%----------------------------------------------------------------------------------

% HOMEOSTASIS calculations:
delay=40
clear H_i
clear H_m
clear H_f
load vec_fittest_dan_100_22_6_99 ;  % composition of fittest
%vec_fittest=n_col_mat/norm( n_col_mat);
for y=1:(cycles-delay)
  vec1=environment(1,:)/norm(environment(1,:)); % initial composition 
  vec2=environment(y,:)/norm(environment(y,:)); % current composition
  vec3=environment(y+delay,:)/norm(environment(y+delay,:)); % comp. in delay steps
  H_i(y)=dot(vec1,vec2);        % similarity to initial
  H_m(y)=dot(vec2,vec3);        % similarity to delayed
  H_f(y)=dot(vec2,vec_fittest); % similarity to fittest
end


%figure;plot(n_col_mat);title('best composition no splitting');
%figure;plot(1:cycles,environment);title('environment');
%figure;plot(1:cycles,concentration);title('concentration');


%figure;subplot(2,1,1);plot(1:cycles,concentration);title('concentration');
%subplot(2,1,2);plot(1:cycles,environment);title('environment');

%figure
%plot(dt*(1:cycles),environment)
figure
plot(dt*(1:cycles),concentration)
set(gca,'fontsize',[18])
figure
plot(dt*(1:(cycles-delay)),H_m,'k--');
hold on; 
%plot(dt*(1:(cycles-delay)),H_i,'k--','linew',[2])
plot(dt*(1:(cycles-delay)),H_f,'r');
set(gca,'fontsize',[18])

drawnow


cont=input('Do correlation plot? (1=yes)');
if cont==1
sampling=input('Every how many steps do sampling?');

n_result=environment;
counterx=0;
for y1=1:sampling:cycles
  counterx=counterx+1;
  countery=0;
  for y2=1:sampling:cycles
    countery=countery+1;
    vec1=n_result(y1,:)/norm(n_result(y1,:));  
    vec2=n_result(y2,:)/norm(n_result(y2,:));
    H(counterx,countery)=dot(vec1,vec2);
  end
end
figure; image((1:sampling:cycles)*dt,(1:sampling:cycles)*dt,H*64)
set(gca,'fontsize',[18],'dataaspectratio',[1 1 1])


end