# File oscillatory_network.ode
#Model of oscillatory network from Figure 4.14. 
#This code can be used to generate Figures 4.15, 4.16, 4.17, 4.21

par k0=8, k1=1, k2=5, n=2
#set n=2.5 for Figures 4.16 and 4.17

#model equations
s1' = k0 - k1*s1*(1+s2^n)
s2' = k1*s1*(1+s2^n) - k2*s2

# initial condition for Figure 4.15A
init s1=1.5, s2=1

# final time
@ total=4
# set the plotting window size:
@ xlo=0, xhi=8, ylo=0, yhi=3.5

# run the simulation in the GUI by selecting (I)nitialconds|(G)o
# plot species s2 by selecting (G)raphic stuff|(A)dd curve and 
# assigning the y-axis

#multiple trajectories can be generated by (G)raphic stuff|(F)reeze and
#running simulations from different initial conditions

#generate a phase plane by choosing (V)iewaxis|(2)D and assigning the axis
#the nullclines can be generated with (N)ullclies|(N)ew
#a direction field can be generated with (D)ir.field/flow|(D)irect field

#the bifurcation diagram in Figure 4.21 can be produced by accessing Auto,
#choosing n as the main parameter, setting the window size and 
#parameter range, and running the continuation.

done