-
Notifications
You must be signed in to change notification settings - Fork 3
/
example5_fit.bngl
202 lines (140 loc) · 5.17 KB
/
example5_fit.bngl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
# A simple model
#
begin model
begin parameters
K2RT__FREE__ 1.12188953e+01
KD1__FREE__ 5.92599705e-01
kdephos__FREE__ 1.59234811e+00
km1__FREE__ 7.83140139e-03
km2__FREE__ 0.0179107743424136
kphos__FREE__ 3.02102168e+00
# simulation parameters
# fraction of a single cell to be considered in a stochastic simulation
f 0.1 # [=] dimensionless, 0<f<=1
# physical constants
# Avogadro constant
NA 6.02214e23 # [=] number of molecules per mol
# Scenario 1: Cells are adherant and grown to confluence in a cell culture dish
# volume of extracellular fluid per cell
# Assumptions:
# - 15 cm culture dishes (cells 95% confluent)
# - 1.2e5 cells per cm^2
# - growth surface area of 140 cm^2
# - 1.7e7 cells per dish
# - 0.25 mL of medium per square cm of growth surface
# - 35 mL of medium per dish
Vecf 2.1e-9*f # [=] L per cell
# ligand concentration
EGF_conc_nM 1.0 # [=] nM
# number of ligands per cell (derived)
# Convert from nM to Molar
EGF_conc = EGF_conc_nM * 1e-9 # M
EGF_copy_number = EGF_conc*(NA*Vecf) # [=] number of molecules per cell
# receptor copy number
EGFR_copy_number f*1.0e5 # [=] number of molecules per cell
# binding parameters
# a free ligand reversibly binds a free receptor site
# equilibrium dissociation constant
# KD1 is specified as being free. It has a matching option in the .conf file.
KD1 KD1__FREE__ # [=] nM
# reverse rate constant (derived)
#km1 is specified as being free. It has a matching option in the .conf file.
km1 km1__FREE__ # [=] /s
# forward rate constant (derived)
kp1=(km1/KD1) # [=] /nM/s
# receptor-receptor interaction
# dimensionless equilibrium association constant
# K2RT is specified as being free. It has a matching option in the .conf file.
K2RT K2RT__FREE__ # [=] dimensionless
# a value of 0.1 => weak crosslinking
# a value of 1.0 => moderate crosslinking
# a value of 10.0 => strong crosslinking
# reverse rate constant (derived)
km2 km2__FREE__ # [=] /s
# forward rate constant (derived)
kp2=K2RT*km2/EGFR_copy_number # [=] /nM/s
kp3=kp2*0.2
# phosphorylation rate constant
# kphos is specified as being free. It has a matching option in the .conf file.
kphos kphos__FREE__ # [=] /s
# dephosphorylation rate constant
# kdephos is specified as being free. It has a matching option in the .conf file.
kdephos kdephos__FREE__ # [=] /s
# A boolean value indicating whether or not EGF is present. This
# is used to set the forward rate constant for ligand binding in
# the "func()" function. In pre-equilibration this value is set to
# 0 to set the rate of ligand binding to 0. After equilibration
# the value will be set to 1 to let ligand bind.
Ligand_isPresent 0
end parameters
begin molecule types
# ligand
L(r)
# receptor
R(l,r,Y~0~P)
end molecule types
#bound ligand and p
begin seed species
L(r) EGF_copy_number
R(l,r,Y~0) EGFR_copy_number
end seed species
begin observables
# total number of ligands
Molecules Ltot L()
# number of free ligands
Species freeL L(r)
# total number of receptors
Molecules Rtot R()
# number of bound ligands = Ltot - freeL
# number of free receptors
#Species freeR R(l,r)
# number of monomeric (unclustered) receptors
#Species Rmon R==1
# number of receptor dimers
Species Rdim R==2
# number of ligand-induced receptor aggregates
# = number of receptor clusters
# = number of complexes containing more than 1 receptor
#Species n_agg_gt1 R>1
# number of ligand-receptor bonds
# = number of ligand-occupied receptor sites
# = number of receptor-occupied ligand sites
Molecules RLbonds L(r!1).R(l!1) # = R(l!+) = L(r!+)
# number of receptors in clusters = Rtot - R1
# average size of a receptor cluster (of size >1)
# = (# of receptors in clusters)/n_agg_gt1
# = (Rtot - R1)/n_agg_gt1
# number of phosphorylated receptors
Molecules pR R(Y~P)
end observables
begin functions
# This function evaluates to 0 if Ligand_isPresent = 0. The function evaluates to kp1 if Ligand_isPresent = 1
# The function is used in the reaction rule governing ligand binding. After equilibration, Ligand_isPresent is
# set to 1 to enable ligand binding.
func() kp1 * Ligand_isPresent
end functions
begin reaction rules
# ligand capture
# a free ligand binds a receptor with a free site
L(r)+R(l)<->L(r!1).R(l!1) func(),km1
# receptor dimerization
R(l!+,r)+R(l!+,r)->R(l!+,r!1).R(l!+,r!1) kp2
# receptor-receptor bond dissociation
R(r!1).R(r!1)->R(r)+R(r) km2
# receptor dimerization (no ligand)
R(l,r) + R(l,r) -> R(l,r!1).R(l,r!1) kp3
# dimer-mediated receptor phosphorylation
R(r!+,Y~0)->R(r!+,Y~P) kphos
# dephosphorylation
R(Y~P)->R(Y~0) kdephos
end reaction rules
end model
# actions
# Generate the network
generate_network();
# Simulate for 600 seconds to reach equilibrium
simulate({method=>"ode",t_start=>0,t_end=>600,n_steps=>1,suffix=>"equil"})
# Set Ligand_isPresent to 1 to enable ligand binding by means of the "func()" function
setParameter("Ligand_isPresent",1)
# Simulate for 60 seconds. This simulation output is fit to the data in example3.exp
simulate({method=>"ode",t_start=>0,t_end=>60,sample_times=>[0,5,10,15,20,25,30,35,40,45,50,55,60],suffix=>"example5"})