#file methionine_model.ode
#from Martinov et al. (2000) Journal of Theoretical Biology vol. 204 pp. 521-532
#generates figures 5.10, 5.11

#declare parameters
pars v_MATI_max=561
pars K_MATI_m=41, K_MATI_i=50
pars v_MAT3_max=22870, K_MIII_m2=21.1
pars v_MET_max=4544, AbyK_METm2=0.1
pars v_GNMT_max=10600, K_GNMT_m=4500, K_GNMT_i=20
pars alpha_d=1333
pars K_AHC=0.1, Adenosine=1
pars Met=48.5

#declare dynamics
AdoMet' = (v_MATI_max * (1/(1+ (K_MATI_m/Met)*(1+AdoMet/K_MATI_i)))) + (v_MAT3_max * (1/(1+ ((20000/(1+ 5.7*(AdoMet/(AdoMet+600))^2))*21.1)/(Met^2+Met*21.1)))) - (v_GNMT_max * (1/(1+(4500/AdoMet)^2.3)) * (1/(1+AdoHcy/20)))-(v_MET_max * (1/(1+ (10 * (1 + AdoHcy/4))/AdoMet + 1/0.1 + (1/0.1)*((10 * (1 + AdoHcy/4))/AdoMet))))

AdoHcy' = (((v_GNMT_max * (1/(1+(4500/AdoMet)^2.3)) * (1/(1+AdoHcy/20)))+(v_MET_max * (1/(1+ (10 * (1 + AdoHcy/4))/AdoMet + 1/0.1 + (1/0.1)*((10 * (1 + AdoHcy/4))/AdoMet))))) - alpha_d * AdoHcy*K_AHC/Adenosine)/(1+K_AHC/Adenosine)

#declare initial conditions
init AdoMet=10, AdoHcy=10

#reduce time step to resolve fast dynamics, and increase allows variable size
@ dt=0.005
@ bound=1000

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

#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




done