function [mCherry,T,tnew,ynew] = v6_ssa_Ecoli_mCherry_full(PqsR_number,PpqsA_number,j,mCherry)

% Gillian Forsyth, Aleksandra Plochocka, Robyn Shuttleworth 
%
% Implementation of the Stochastic Simulation Algorithm
% (or Gillespie's algorithm) on the PQS system which
% was developed by Dundee iGEM Modelling team.
% D J Higham, 2006.
% 
% The function uses Gillespie's algorithm to calculate numerically the Stochastic change 
% in mCherry numbers, for the jth realisation and the given
% PqsR_conc and PpqsA_conc input configuration. The change is 
% based on our system over the selected time interval and the values are 
% stored in the jth column of mCherry matrix. 
% The function passes additional plot data and the modified mCherry matrix 
% into programs calling it. 
% Inputs: 
%      PqsR_number: amount of PqsR, receptor, molecules
%      PpqsA_number: amount of PpqsA, promoter, molecules
%      j: realisation number
%      mCherry: matrix to store mCherry numbers for each time and realisation
% 
% Outputs:
%      mCherry: matrix with mCherry numbers for each time and realisation
%      T: uniform time vector for the time interval
%      ynew: is the non-interpolated mCherry numbers
%      tnew: non-uniform time corresponding to reaction time for ynew
%
% Usage:
% [mCherry,T,tnew,ynew] = v6_ssa_Ecoli_mCherry_full(PqsR_number,PpqsA_number,j,mCherry)


V = [ 1 -2 2 0 0 -1 0;
      0 -1 1 0 0 0 0;
      0 1 -1 -1 1 0 0;
      0 0 0 -1 1 0 0;
      0 0 0 1 -1 0 0;
      0 0 0 0 0 0 1];
  
T = 1:1200;                                     
T=T'; 

nA = 6.023e23;                                  % Avagadro's number
vol = 1e-15;                                    % E.coli volume

Se_conc = 0.2e-6;                                % concentration of PQS in sputum
S = (Se_conc * nA * vol);                       % corresponding number of molecules 

%%%%% Deterministic Rates

k = 0.001/60;                                   % PQSe movement into the cell (/s)
k1 = 4.76/60;                                   % PQSi and PQSR association rate (/M/s)
k2 = 1/60;                                      % PQSi and PQSR dissociation rate (/s)
k3 = 1/60;                                      % PQSi2.PQSR and PpqsA association rate (/M/s)
k4 = 2/60;                                      % PQSi2.PQSR and PpqsA dissociation rate (/s)
d = 0.001/60;                                   % PQSi movement out of the cell (/s)
K = 1/60;                                       % Maximal rate of mCherry expression per promoter (/s)

%%%% corresponding Stochastic rates

c(1) = k;
c(2) = k1*(nA*vol);
c(3) = k2;
c(4) = k3*(nA*vol);
c(5) = k4;
c(6) = d;
c(7) = K;

X = zeros(6,1); 

X(2) = PqsR_number;
X(4) = PpqsA_number; 

t = 0;                                           
tfinal = 1200;

count = 1;
tvals(1) = 0;
Xvals(:,1) = X;

while t < tfinal
    % Propensity functions correspond to picking a discrete random variable
    % 
         a(1) = c(1)*S;
         a(2) = c(2)*X(1)*X(1)*X(2);
         a(3) = c(3)*X(3);
         a(4) = c(4)*X(3)*X(4);
         a(5) = c(5)*X(5);
         a(6) = c(6)*X(1);
         a(7) = c(7)*X(5);
         
         asum = sum(a);
         b = min(find(rand<cumsum(a/asum)));
         % selects reaction to take place
         tau = log(1/rand)/asum;
         % time for reaction
         X = X + V(:,b);
         
         count = count + 1;
         t = t + tau;
         tvals(count) = t;
         Xvals(:,count) = X;
        
    
end

%%%%%%%%%% Plot data %%%%%%%%%%

% Creates a non-uniform time matrix to store reaction time for
% corresponding mCherry matrix, ynew.
Length = length(tvals);
tnew = zeros(1,2*(Length-1));
tnew(1:2:end-1) = tvals(2:end);
tnew(2:2:end) = tvals(2:end);
tnew = [tvals(1),tnew];

tnew=tnew';

% Creates a non-interpolated matrix of mCherry numbers for the number of j
% realisations. 
mCherry_vals = Xvals(6,:);
ynew = zeros(1,2*Length-1);
ynew(1:2:end) = mCherry_vals;
ynew(2:2:end-1) = mCherry_vals(1:end-1);

ynew=ynew';


%%%%%%%%%% Interpolation %%%%%%%%%%

mCherry(:,j) = interp1q(tnew,ynew,T);

end

