function [ t, s ] = spatial_model( initial_cell, time, grid_size, HSL_pos )
% function [ t, s ] = spatial_model( initial_cell, time, grid_size, HSL_pos )
%
% Runs a spatial simulation of the interval switch model with specific
% starting concentrations for each species in each cell an a user-defined
% grid. A 'sender cell' can also be added so a ring of GFP can be observed
% where the concentration of HSL is in between a certain margin.
%
% INPUTS:
%
% initial_cell: must be an array of length 4, representing the concentration
%               (in uM) of each species for each cell. The indexes 
%               correspond to:
%     GFP = 1;
%     HSL = 2;
%     cI = 3;
%     HSLluxR = 4;
%
% time: is the (virtual) simulation time, in minutes
%
% grid_size: a two-dimensional array representing the width and height of
%            the grid on which the simulation will run
%
% HSL_pos: an array representing the positions of the 'sender cells', where
%          the cell positions are (HSL_pos(1),HSL_pos(2)),
%          (HSL_pos(3),HSL_pos(4)), etc.
%
% OUTPUTS:
%
% t: an array with t(I) the time at time-step I
%
% s: the concentration of a certain species Sp of cell (I,J) at time-step T is
%    stored in S(T,I,J,Sp)
%

%----- Parameters -----
V_GFP = 2; % maximal GFP expression
V_cI = 1;  % maximal cI expression
V_HSL = 1; % maximal HSL production sender cell

Km = 0.01;   % HSLluxR concentration for half-maximal expression of GFP
Kmi = 0.008; % cI concentration for half-maximal expression of GFP
Km_ = 0.08;  % HSLluxR concentration for half-maximal expression of cI, assumed to be 8 times higher then the unmodified hybrid promoter

n = 1;  % HSLluxR cooperativity constant
ni = 2; % cI cooperativity constant

k_dec_cI = 0.0692;      % decay cI
k_dec_HSL = 0.01;       % decay HSL
k_dec_GFP = 0.0692;     % decay GFP
k_dec_HSLluxR = 0.0231; % decay HSLluxR

k_diff_HSL = 0.001; % diffusion speed HSL

k_dim_HSLluxR = 0.5; % rate at which HSL binds to luxR

luxR = 0.5; % luxR concentration is kept constant at 0.5

%----- Species (indices of species array) -----
n_species = 4; % total number of species in the model

gfp = 1;
HSL = 2;
cI = 3;
HSLluxR = 4;

%----- Prepare and run simulation ----- 

% increase grid size to add padding
x = grid_size(1)+2;
y = grid_size(2)+2;

z = n_species;

% set initial values
initial_values = zeros(x,y,z);
for index = 1:n_species
  initial_values(:,:,index) = initial_cell(index);
end
initial_values = reshape(initial_values,1,x*y*z);
% start simulation
[t,s] = ode15s(@cell_sim,[0:time],initial_values);

% reshape the returned array
s = reshape(s,time+1,grid_size(1)+2,grid_size(2)+2,n_species);
% remove the padding
%s(:,1,:,:) = [];
%s(:,end,:,:) = [];
%s(:,:,1,:) = [];
%s(:,:,end,:) = [];

  function d_cells = cell_sim(t,cells)
    cells(find(cells < 0)) = 0;
    cells = reshape(cells,x,y,z);
    d_cells = zeros(x,y,z);

    for I = 2:(x-1)
      for J = 2:(y-1)
        % calculation production rates (derivatives of species concentrations)
        d_cells(I,J,gfp)  = V_GFP*cells(I,J,HSLluxR)^n/(Km^n+cells(I,J,HSLluxR)^n)*(1/(1+(cells(I,J,cI)/Kmi)^ni)) - k_dec_GFP*cells(I,J,gfp);
        d_cells(I,J,cI)   = V_cI*cells(I,J,HSLluxR)^n/(Km_^n+cells(I,J,HSLluxR)^n) - k_dec_cI*cells(I,J,cI);
        d_cells(I,J,HSLluxR) = k_dim_HSLluxR*cells(I,J,HSL)^2*luxR^2 - k_dec_HSLluxR*cells(I,J,HSLluxR);
        % calculate HSL diffusion to and from the 4 orthogonal neighbors
        d_cells(I,J,HSL) = k_diff_HSL*(cells(I-1,J,HSL)+cells(I,J-1,HSL)+cells(I+1,J,HSL)+cells(I,J+1,HSL)) - k_diff_HSL*4*cells(I,J,HSL) - k_dec_HSL*cells(I,J,HSL);

        % production of HSL by sender populations
        for S = 1:2:length(HSL_pos)
          % if this is a sender cell, make it produce the specified amount of HSL
          if(I == (HSL_pos(S) + 1) && J == (HSL_pos(S+1) + 1))
            d_cells(I,J,HSL) = d_cells(I,J,HSL) + V_HSL;
            break; % quit for-loop
          end
        end
        
      end
    end

    d_cells = reshape(d_cells,x*y*z,1,1);
  end

end