toggle.bngl

## Title: A model of a genetic toggle switch
## Description:   In this model, a transcriptional repressor, X, inhibits the synthesis
#                 of a second repressor, Y, which inhibits the synthesis of X.
#                 Repression is captured using the following function: 
#                 1/(1+([repressor]/K)^n).
## Author: William S. Hlavacek
## Date: 2015
## Reference: Gardner TS, Cantor CR, Collins JJ (2000) Nature 403: 339-342.


begin model

begin parameters
# Note that parameter values are different from those chosen in Ref 1.\
# maximum rate of synthesis of repressor X
bx 50 # [=] molecules per cell per unit time
# maximum rate of synthesis of repressor Y
by 50 # [=] molecules per cell per unit time
# rate constant for clearance of X (through degradation and dilution)
ax 1 # [=] per unit time
# rate constant for clearance of Y (through degradation and dilution)
ay 1 # [=] per unit time
# parameters in the function characterizing Y's effect on synthesis of X
nyx 3 # [=] dimensionless
Kyx 20 # [=] same units as Y()
# parameters in the function characterizing X's effect on synthesis of Y
nxy 3 # [=] dimensionless
Kxy 20 # [=] same units as X()
# initial abundance of X 
Xinit 0 # [=] copies per cell
# initial abundance of Y
Yinit 0 # [=] copies per cell
end parameters

begin species
X() Xinit
Y() Yinit
end species

begin observables
Molecules X_tot X()
Molecules Y_tot Y()
end observables

begin reaction rules
# synthesis of X according to a user-defined rate law
#   0 is a source.
0->X() bx/(1+(Y_tot/Kyx)^nyx) 
# clearance of X (through degradation and dilution)
#   0 is a sink.
X()->0 ax
# synthesis of Y according to a user-defined rate law
0->Y() by/(1+(X_tot/Kxy)^nxy) 
# clearance of Y
Y()->0 ay
end reaction rules

end model

generate_network({overwrite=>1})
# Write ODEs for chemical kinetics of the generated network in a MATLAB M-file 
writeMfile()

# set initial system state close to the stable steady state where
#   expression of repressor X is off.
setConcentration("X()",0)
setConcentration("Y()",50)

# find steady states as a function of Kxy, which reflects the strength
#   of X's repression of Y.
#   Values of Kxy are scanned from the minimum specified value to the maximum.
#   Values are also scanned in the opposite direction.
#   If bistability exists, this procedure will find both steady states at each
#   tested value of the bifurcation parameter in the bistable region.
bifurcate({parameter=>"Kxy",par_min=>1.0,par_max=>30,n_scan_pts=>100,\ 
           log_scale=>1,method=>"ode",t_start=>0,t_end=>1000,n_steps=>10,\
           steady_state=>1})

# If the parameters of the system are such that the system is bistable, running a long 
# stochastic trajectory may reveal transitions between the two steady states.
#simulate({method=>"ssa",t_end=>100000,n_steps=>1000})