Mathematica-File  for the experiments of Subsection 3.7 of the paper “An effective model of endogeneuos clocks and external stimuli” by Tim Breitenbach, Charlotte Helfrich-Förster and Thomas Dandekar

(*Generate unperturbed time curves of the agents of the model*)
d=120; (*time Horizon for the simulation*)
s=NDSolve[{MD'[t]==0.76*(1^4)/(1^4+(PN[t])^4)-0.65*MD[t]/(0.5+MD[t]),
P0'[t]==0.38*MD[t]-3.2*P0[t]/(2+P0[t])+1.58*P1[t]/(2+P1[t]),
P1'[t]==3.2*P0[t]/(2+P0[t])-1.58*P1[t]/(2+P1[t])-5P1[t]/(2+P1[t])+2.5*P2[t]/(2+P2[t]),P2'[t]==5*P1[t]/(2+P1[t])-2.5*P2[t]/(2+P2[t])-1.9*P2[t]+1.3*PN[t]-0.95*P2[t]/(0.2+P2[t]),
PN'[t]==1.9*P2[t]-1.3PN[t],MD[0]==0.5,P0[0]==0.5,P1[0]==0.5,P2[0]==0.6,PN[0]==1.5},{MD[t],P0[t],P1[t],P2[t],PN[t]},{t,0,d}];
Plot[Evaluate[{MD[t],P0[t],P1[t],P2[t],PN[t]}/.s],{t,0,d},PlotRange->All,PlotLegends->{M,(Subscript[P, 0]),(Subscript[P, 1]),(Subscript[P, 2]),(Subscript[P, N])}];
ZwM=Evaluate[MD[t]/.s];
ZwP0=Evaluate[P0[t]/.s];
ZwP1=Evaluate[P1[t]/.s];
ZwP2=Evaluate[P2[t]/.s];
ZwPN=Evaluate[PN[t]/.s];

(*Model for Subsection 3.7, Anti-phase coupling*)
u[t_]:=(Cos[2*Pi*t/24+Pi]+1)*0.25
s=NDSolve[{M'[t]==0.76*(1^4)/(1^4+(PN[t])^4)-0.65*M[t]/(0.5+M[t]),
P0'[t]==0.38*M[t]-3.2*P0[t]/(2+P0[t])+1.58*P1[t]/(2+P1[t])-u[t]*P0[t],
P1'[t]==3.2*P0[t]/(2+P0[t])-1.58*P1[t]/(2+P1[t])-5P1[t]/(2+P1[t])+2.5*P2[t]/(2+P2[t]),P2'[t]==5*P1[t]/(2+P1[t])-2.5*P2[t]/(2+P2[t])-1.9*P2[t]+1.3*PN[t]-0.95*P2[t]/(0.2+P2[t]),
PN'[t]==1.9*P2[t]-1.3PN[t],
MB'[t]==0.76*(1^4)/(1^4+(PBN[t])^4)-0.65*MB[t]/(0.5+MB[t]),
PB0'[t]==0.38*MB[t]-3.2*PB0[t]/(2+PB0[t])+1.58*PB1[t]/(2+PB1[t])-0.5*P0[t]*PB0[t],
PB1'[t]==3.2*PB0[t]/(2+PB0[t])-1.58*PB1[t]/(2+PB1[t])-5PB1[t]/(2+PB1[t])+2.5*PB2[t]/(2+PB2[t]),PB2'[t]==5*PB1[t]/(2+PB1[t])-2.5*PB2[t]/(2+PB2[t])-1.9*PB2[t]+1.3*PBN[t]-0.95*PB2[t]/(0.2+PB2[t]),
PBN'[t]==1.9*PB2[t]-1.3PBN[t],M[0]==0.5,P0[0]==0.5,P1[0]==0.5,P2[0]==0.6,PN[0]==1.5,MB[0]==0.5,PB0[0]==0.5,PB1[0]==0.5,PB2[0]==0.6,PBN[0]==1.5},{M[t],P0[t],P1[t],P2[t],PN[t],MB[t],PB0[t],PB1[t],PB2[t],PBN[t]},{t,0,d}];
Plot[Evaluate[{M[t],P0[t],P1[t],P2[t],PN[t]}/.s],{t,0,d},PlotRange->All,PlotLegends->{(M^u),((Subscript[P, 0])^u),((Subscript[P, 1])^u),((Subscript[P, 2])^u),((Subscript[P, N])^u)}]
Plot[Evaluate[{MB[t],PB0[t],PB1[t],PB2[t],PBN[t]}/.s],{t,0,d},PlotRange->All,PlotLegends->{OverTilde[M],(Subscript[OverTilde[P], 0]),(Subscript[OverTilde[P], 1]),(Subscript[OverTilde[P], 2]),(Subscript[OverTilde[P], N])}]
Plot[{Evaluate[{M[t],u[t],MB[t]}/.s],ZwM},{t,0,d},PlotRange->All,PlotLegends->{(M^( u)),u,OverTilde[M],M}]