# File symmetric_network.ode
#Model of symmetric network from Figure 4.6. 
#This code can be used to generate Figures 4.7, 4.8, 4.9, and 4.19

par k1=20, k2=5, k3=5, k4=5, k5=2, n1=4, n2=1
#set n2=4 for Figures 4.8, 4.9 and 4.19A

#model equations
s1' = k1/(1+s2^n1) - k3*s1
s2' = k2/(1+s1^n2) - k4*s2

# initial condition for Figure 4.7A
init s1=3, s2=1

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

# 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.19B can be produced by accessing Auto,
#choosing k1 as the main parameter, setting the window size and 
#parameter range, and running the continuation.

done