function notch_ring

%% simulate a simple Delta-Notch signalling model in a ring of cells

%% model parameters

n=50; %number of cells

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
init=[namp*rand(n,1);damp*rand(n,1)];

%% run the simulation -- do not change

sig=0.5*(circshift(eye(n),1)+circshift(eye(n),-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

 for j=1:nrep
     ringplotter(n,stout(:,j),tout(j))
 end

%% space-time plots of solution

figure(11), clf

subplot(2,1,1)
imagesc(stout(1:n,:)'), axis xy
set(gca,'FontSize',14)
colormap hot, colorbar
xlabel('cell number')
ylabel('time')
title('Notch activity')
nt=get(gca,'YTickLabel')
set(gca,'YTickLabel',linspace(tstart,tend,length(nt)))


subplot(2,1,2)
imagesc(stout(n+1:2*n,:)'), axis xy
set(gca,'FontSize',14)
colormap hot, colorbar
xlabel('cell number')
ylabel('time')
title('Delta activity')
nt=get(gca,'YTickLabel')
set(gca,'YTickLabel',linspace(tstart,tend,length(nt)))

%% plot nullclines for 2 cell quasi steady state system

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

figure(12), 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);
        delta=y(n+1:2*n);
        dndt=f(sig*delta)-notch;
        dddt=g(notch)-delta;
        dydt=[dndt;dddt];
    end

    function dydt=model_ddes(t,y,z)
        notch=y(1:n);
        delta=y(n+1:2*n);
        dnotch=z(1:n);
        dndt=f(sig*delta)-notch;
        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

%% function for plotting solutions on a ring of cells -- do not change

function ringplotter(n,d,t)

%given data d on a ring of n cells
%plot value of d in each cell

notch=d(1:n);
delta=d(n+1:2*n);

a=2*pi/n;
r1=1;
r2=1.2;
figure(10)
clf
subplot(1,2,1)
for i=1:n
    patch([r1*cos((i-1)*a),r2*cos((i-1)*a),r2*cos(i*a),r1*cos(i*a)],[r1*sin((i-1)*a),r2*sin((i-1)*a),r2*sin(i*a),r1*sin(i*a)],notch(i)*[1,1,1])
end
axis tight, axis off
set(gca,'FontSize',16)
set(gca,'DataAspectRatio',[1 1 1])
title(['Notch activity. Time = ',num2str(t)])
colorbar('SouthOutside')
colormap gray

subplot(1,2,2)
for i=1:n
    patch([r1*cos((i-1)*a),r2*cos((i-1)*a),r2*cos(i*a),r1*cos(i*a)],[r1*sin((i-1)*a),r2*sin((i-1)*a),r2*sin(i*a),r1*sin(i*a)],delta(i)*[1,1,1])
end

axis tight, axis off
set(gca,'FontSize',16)
set(gca,'DataAspectRatio',[1 1 1])
title(['Delta activity. Time = ',num2str(t)])
colorbar('SouthOutside')
colormap gray

end

end