# “Unlimited multistability and Boolean logic in microbial signaling”, Varun B. Kothamachu, Elisenda Feliu, Luca Cardelli, and Orkun S Soyer
# Supplementary file 2.
# This file contains a mathematical (ordinary differential equations) model for the yeast two-component system involved in osmosensing, as depicted in Figure 5A of the main paper.
# This model is executable in XPPAUT, a modeling tool capable of bifurcation analysis. XPPAUT is available online at http://www.math.pitt.edu/~bard/xpp/xpp.html

# In this osmosensing relay, there is one HK (known as Sln1), one Hpt (known as Ypd) and two RRs (known as Sln7 and Ssk1).
# The Sln1 is a hybrid HK (i.e. the HK and REC proteins of a regular phosphorelay are combined). Hence, we model Sln1 as a protein with two phosphorylation sites.
# In the model, we denote by:

# O: un-phosphorylated site. 
# P: phosphorylated site.
# H1: Sn1_00
# H2: Sln1_P0
# H3: Sln1_0P
# H4: Sln1_PP 
# T: Ypd
# T: Ypdp
# RR1: Ssk1
# RR1p: Ssk1p
# RR2: Sln7
# RR2p: Sln7p
# Htot: Sln(Total)
# Ttot: Ypd(Total)
# RR1tot: Ssk1(Total)
# RR2tot: Sln7(Total)

# The model parameters are listed in supplementary table S2 and are mostly derived from the following paper:
# Janiak-Spens F, Cook P F, West A H. “Kinetic Analysis of YPD1-Dependent Phosphotransfer Reactions in the Yeast Osmoregulatory”
# Phosphorelay System. Biochemistry 44, no. 1 (2005): 377-86. doi:10.1021/bi048433s.

# The result shown in Figure 5 of the main text are created by running a bifurcation analysis using parameter k1 as the bifurcation parameter.


####### REACTIONS ########
# H1-k1-> H2  # autophosphorylation
# H3-k3-> H4  # autophosphorylation
# H2-k2->H3   # internal phosphotransfer
# H3-k2r->H2  # internal reverse phosphotransfer

# H3+T<k7-k4>Tp+H1   	# forward and reverse phosphotransfer between HK and Hpt   
# H4+T<k3rb-k3b>Tp+H2	# forward and reverse phosphotransfer between HK and Hpt

# RR1+Tp<k12-k11>RR1p+T # forward and reverse phosphotransfer between Hpt and RR1
# RR2+Tp<k8-k5>RR2p+T 	# forward and reverse phosphotransfer between Hpt and RR2

# H3 -k9->H1 		# hydrolysis (i.e. dephosphorylation) from H3
# H4 -k10->H2		# hydrolysis (i.e. dephosphorylation) from H4 
# RR1p -k13-> RR1+p	# hydrolysis (i.e. dephosphorylation) from RR1
# RR2p -k14-> RR2+p	# hydrolysis (i.e. dephosphorylation) from RR2



###### PARAMETER VALUES ######
#input parameter
p k1=0.0001,m=1,k3=m*k1

# forward and reverse phosphorylation rates in the HK
p k2=160,k2r=0,k4=30,k7=230
p k3b=30,k3rb=230

# forward and reverse phosphorylation rates downstream of the HK
p k11=1.4,k12=0.4,k5=160,k8=0

# hydrolysis rate/autodephosphorylation rate
p k9=0.1,k10=1,k13=0.1,k14=0.1

#total protein concentrations in all layers
p Htot=0.25,Ttot=2,RR1tot=1.5,RR2tot=1.5



#### INITIAL CONDITIONS ####
init H1=0
init H3=0,H4=0.1
init Tp=0.5,RR1p=0.01,RR2p=0.01



###### SYSTEM OF EQUATIONS ######

dH1/dt=k4*H3*(Ttot-Tp)+k9*H3-(k1+k7*Tp)*H1

dH3/dt=k2*(Htot-H1-H3-H4)+k7*Tp*H1-(k3+k2r+k4*(Ttot-Tp)+k9)*H3

dH4/dt=k3*H3+k3rb*Tp*(Htot-H1-H3-H4)-(k3b*H4*(Ttot-Tp)+k10)*H4

dTp/dt=(k4*H3+k3b*H4+k12*RR1p+k8*RR2p)*(Ttot-Tp)-(k7*H1+k3rb*(Htot-H1-H3-H4)+k11*(RR1tot-RR1p)+k5*(RR2tot-RR2p))*Tp

dRR1p/dt=k11*Tp*(RR1tot-RR1p)-(k12*(Ttot-Tp)+k13)*RR1p

dRR2p/dt=k5*Tp*(RR2tot-RR2p)-(k8*(Ttot-Tp)+k14)*RR2p


@ dt=0.001,maxstor=100000,math=stiff,total=100,bound=1000
done