function notch_square

%% simulate a simple Delta-Notch signalling model on an array of square cells

%% model parameters

n=15; %number of cells in x-direction
m=15; %number of cells in y-direction

thresh1 = 0.2; %N activation threshold
thresh2 = 0.1; %D repression threshold
n1 = 2; %N activation sensitivity
n2 = 2; %D repression sensitivity
delay = 0; %time delay

%% simulation settings

tstart=0;  %start time
tend=50;   %end time
nrep=100;  %number of solution reports

%% initial conditions

namp = 0;   %Notch
damp = 0.1; %Delta
ninit=namp*rand(n,m);
dinit=damp*rand(n,m);
init=[ninit(:);dinit(:)];

%% run the simulation -- do not change

nsig=circshift(eye(n),1)+circshift(eye(n),-1);
msig=circshift(eye(m),1)+circshift(eye(m),-1);

if delay==0
    sol=ode45(@model_odes,[tstart,tend],init);
else
    sol=dde23(@model_ddes,delay,init,[tstart,tend]);
end
tout=linspace(tstart,tend,nrep);
stout=deval(sol,tout);

%% plot the solution at nrep regular time intervals

figure(31),clf
for j=1:nrep
    nout=reshape(stout(1:n*m,j),n,m);
    dout=reshape(stout(n*m+1:2*n*m,j),n,m);
    imagesc(nout'),axis xy,colorbar,colormap hot
    caxis([min(min(stout(1:n*m,:))),max(max(stout(1:n*m,:)))]);
    title(['Notch activity. Time = ',num2str(tout(j))])
    pause(0.1)
end

%% plot nullclines for 2 cell quasi steady state system

x=0:0.01:1.1;
gf=g(f(x));

figure(32), clf
set(gca,'FontSize',14)
plot(x,gf,'b',gf,x,'r','LineWidth',2)
axis tight
xlabel('D_1'), ylabel('D_2')
legend('D_2 nullcline','D_1 nullcline')

%% specification of the ODEs or DDEs

    function dydt=model_odes(t,y)
        notch=y(1:n*m);
        delta=y(n*m+1:2*n*m);
        nm=reshape(notch,n,m);
        dm=reshape(delta,n,m);
        dsig=0.25*(nsig*dm + dm*msig);
        dndtm=f(dsig)-nm;
        dndt=dndtm(:);
        dddt=g(notch)-delta;
        dydt=[dndt;dddt];
    end

    function dydt=model_ddes(t,y,z)
        notch=y(1:n*m);
        delta=y(n*m+1:2*n*m);
        dnotch=z(1:n*m);
        nm=reshape(notch,n,m);
        dm=reshape(delta,n,m);
        dsig=0.25*(nsig*dm + dm*msig);
        dndtm=f(dsig)-nm;
        dndt=dndtm(:);
        dddt=g(dnotch)-delta;
        dydt=[dndt;dddt];
    end

%% regulation functions -- do not change

    function ff=f(x)
        ff=(x.^n1)./(thresh1.^n1 + x.^n1);
        %increasing Hill function
    end

    function gg=g(x)
        gg=1./(1 + (x./thresh2).^n2);
        % decreasing Hill function
    end

end