RuleWorld / bionetgen

Rule-based modeling framework
https://bionetgen.org/
MIT License
59 stars 25 forks source link

Reactions not found #205

Closed iustin94 closed 6 years ago

iustin94 commented 6 years ago

I am having an issue with the network generation, I suspect.

The message I get is :

BioNetGen` version 2.3.1 Reading from file /home/iustinian/RuleBender-workspace/Dynamin_I/Dynamin/main.bngl (level 0) Read 10 parameters. Read 7 molecule types. Read 4 species. Read 6 observable(s). Read 4 functions. Read 6 reaction rule(s). ACTION: generate_network( main ) Iteration 0: 4 species 0 rxns 0.00e+00 CPU s Iteration 1: 4 species 0 rxns 0.00e+00 CPU s Cumulative CPU time for each rule Rule 1: 0 reactions 0.00e+00 CPU s 0.00e+00 CPU s/rxn Rule 2: 0 reactions 0.00e+00 CPU s 0.00e+00 CPU s/rxn Rule 3: 0 reactions 0.00e+00 CPU s 0.00e+00 CPU s/rxn Rule 4: 0 reactions 0.00e+00 CPU s 0.00e+00 CPU s/rxn Rule 5: 0 reactions 0.00e+00 CPU s 0.00e+00 CPU s/rxn Rule 6: 0 reactions 0.00e+00 CPU s 0.00e+00 CPU s/rxn Total : 0 reactions 0.00e+00 CPU s 0.00e+00 CPU s/rxn ABORT: writeFile() was asked to write the network, but no reactions were found. Did you remember to call generate_network() before writing network output? at line 92

I am not sure why this is happening. The error is given at the generate_network() line which is a bit confusing given that I am calling the method but I am told that I didn't? My simulation script is the following:

begin model

begin parameters
NaV 6.02e8 # Conversion constant: M -> #/um^3
Vcell 1000 # Typical eukaryotic cell volume ~ 1000 um^3
Vec   1000*Vcell # Volume of extracellular space around each cell (1/cell density)
lig_conc 1e-8

# Get forward binding rates for Dynamin binding sites
kp1 1/(NaV*Vec) # Forward binding rate constant for Phopshorilation
km1 1/(NaV*Vec) # Reverse binding rate constant for L-R

#   Initial concentrations
DYN1_0 1e-8*NaV*Vec# uM -> molecules/cell
#ACTA1_0 0.2*NaV*Vcell # uM -> molecules/cell
#FNBP17_0 0.05*NaV*Vcell # 
#GRB2_0 0.4*NaV*Vcell 
P_0 1.2e-8*NaV*Vec
CaN_0 1e8*NaV*Vec
Cdk5_0 1e8*NaV*Vec

end parameters

# Proteins that interact together only in the cell 
begin molecule types

 # Dynamin I homo-sapiens
 # GTPase~O unphosphorilated
 # GTPase~P phosphorilated
DYN1(P,GTPase~O~P,PRD)

CaN(d)
Cdk5(d)

ACTA1(DYN1) #Actin alpha skeletal muscleW
FNBP17(SH3,FBAR,REM1) #Formin-binding protein 1
GRB2(SH2,SH3,SH2) #Growth factor receptor bound protein
P(l) #Phosphate
end molecule types

begin seed species
DYN1(P,GTPase~P,PRD) DYN1_0
P(l) DYN1_0
CaN(d) DYN1_0 
Cdk5(d) DYN1_0

#DYN1(P,GTPase~O,PRD!1).Cdk5(d!1) DYN
# ACTA1(DYN1) ACTA1_0
#FNBP17(SH3,FBAR,REM1) FNBP17_0
#GRB2(SH2,SH3,SH2) GRB2_0

end seed species

begin observables 
Molecules DU DYN1(P,GTPase~P,PRD)
Molecules DP DYN1(P!1,GTPase~O,PRD).P(l!1)
Molecules CD DYN1(P,GTPase~O,PRD!2).CaN(d!2)
Molecules CdD DYN1(P,GTPase~O,PRD!1).Cdk5(d!1)
Molecules Cdk Cdk5(d)
Molecules Can CaN(d)
end observables

#Simple interdependency to make some waves
begin functions
pFunction()= (DP/DU)*1e8*NaV*Vec
uFunction()= (DU/DP)*1e8*NaV*Vec
cdk5BindR()= (DP+1/DU)*Cdk*1e8*NaV*Vec
canBindR()= uFunction()*Can
end functions

begin reaction rules
#Phosphorilation
Cdk5Bingind: DYN1(P,GTPase~O,PRD)+Cdk5(d)->DYN1(P,GTPase~O,PRD!1).Cdk5(d!1) cdk5BindR()
RePhopshorilation:DYN1(P,GTPase~O,PRD!1).Cdk5(d!1)+P(l)->DYN1(P!1,GTPase~P,PRD!2).P(l!1).Cdk5(d!2) pFunction()
CdkUnbund: DYN1(P!1,GTPase~P,PRD!2).P(l!1).Cdk5(d!2)-> DYN1(P!1,GTPase~P,PRD).P(l!1)+ Cdk5(d) cdk5BindR()*1e-1 

#Dephosphorilation
CaNBinging:DYN1(P!1,GTPase~P,PRD).P(l!1)+CaN(d)-> DYN1(P!1,GTPase~P,PRD!2).P(l!1).CaN(d!2) canBindR()
DePhosphorilation: DYN1(P!1,GTPase~P,PRD!2).P(l!1).CaN(d!2)-> DYN1(P,GTPase~O,PRD!2).CaN(d!2)+ P(l) uFunction()
CaNUnbinding: DYN1(P,GTPase~O,PRD!2).CaN(d!2) -> DYN1(P,GTPase~O,PRD) + CaN(d) canBindR()*1e-1

#Actin_bind: DYN1(P,GTPase~P,Middle,PH,GED,PRD)+ ACTA1(DYN1) ->DYN1(P,GTPase~P,Middle!1,PH,GED,PRD).ACTA1(DYN1!1) DA_kd
#FNBP17_bind: DYN1(P,GTPase~P,Middle,PH,GED,PRD)+ FNBP17(SH3) -> DYN1(P,GTPase~P,Middle,PH,GED,PRD!1).FNBP17(SH3!1) DF_kd
#GRB2_bind: DYN1(P,GTPase~P,Middle,PH,GED,PRD)+GRB2(SH3) -> DYN1(P,GTPase~P,Middle,PH,GED,PRD!1).GRB2(SH3!1) DG_kd
end reaction rules
end model

#Actions
generate_network({overwrite=>1})
writeModel()
simulate({method=>"ssa",t_end=>1e-5,n_steps=>200})
jrfaeder commented 6 years ago

What is happening is that there are no feasible reactions that can occur from the seed set of species as you have defined them. I believe you probably meant to define the initial state of DYN1(GTPase) as 'O', but you have defined it as 'P' instead: DYN1(P,GTPase~P,PRD) DYN1_0 in the seed species block. When I change this to DYN1(P,GTPase~O,PRD) DYN1_0 a network with 9 species and 6 reactions is generated.