Hybrid particle-population model generator
From BioNetWiki
Contents |
Managing combinatorial complexity and Large populations
NFsim is a network-free kinetic Monte Carlo simulator for biochemical models. NFsim permits efficient simulation of models with large or infinite state spaces. NFsim is similar to agent-based simulators in that each molecule in the system is represented as an individual object. Each reactant site is tracked by NFsim and updated as the system evolves. Thus, NFsim does not need to ennumerate the state space, but only tracks the reactions possible in the system in the present state. One disadvantage of this approach is that memory usage increases linearly with the population of molecules in the system. For many models, the memory available on PC computing resources is sufficient. However, for large scale models memory usage may become a barrier to simulation.
Rule-based models of biochemical systems often feature combinatorially complex molecular aggregates that induce exponential blow-up of the state space. A receptor with multiple modification and binding sites is an example of a molecule with combinatorial complexity. Furthermore, it is common to find a subset of molecule types that appear in very large numbers but have relatively few states. A soluble ligand is an example of a molecule with only a few states (bound/unbound) and a large population. The presence of combinatorial features necessitates the use of a network-free simulator, however is seems a waste of resources to individually represent thousands to millions of simple ligand molecules. In this case a hybrid approach, in which a subset of molecular species are lumped as population variables and the remainder are treated indivdiually, may be the optimal simulation approach.
Hybrid model generator
The hybrid model generator transforms arbitrary BNGL models into a hybrid particle-population form with identical dynamics. The hybrid model is written to the file [model]_hybrid.bngl, by default. The user is responsible for selecting species to handle as populations. The latest versions of NFsim are able to simulate hybrid models, treating population species as lumped variables (small memory footprint) and other species as individual particle objects.
The hybrid model is created by expanding all reaction rules such that the action on population species as reactants is explicit. Additionally, a set of gathering rules are constructed that identify population species that are created a reaction products and then lump them into population variables. The combination of partial network expansion and gathering rules permits us to represent population species as simple, unstructured population variables while preserving the precise dynamics of the rule-based model.
How to create and simulate a hybrid model with BioNetGen and NFsim
0) Start with an arbitary BNGL model.
1) Choose the species that will be lumped into population variables. For example, in the TLBR model below:
- L(r,r,r)
- R(l,l)
2) Add a population maps block to the BNGL model file. This block defines the population species and associates each with a population variable name.
begin population maps L(r,r,r) -> L_pop k_gather R(l,l) -> R_pop k_gather end population maps
3) Define a gathering rate parameter "k_gather" in the parameter block. Assign a large value in comparison to the time scale of the model. Ideally this parameter would be infinite, but BNG/NFsim does not yet support an infinite rate. Small values may lead to hybrid simulation error.
begin parameters ... k_gather 1000000 end parameters
4) Add the "generate_hybrid_model" action to the BNGL model file:
generate_hybrid_model({ overwrite=>1, execute=>1, verbose=>1 \
actions=>['simulate_nf({suffix=>nf,complex=>1,\
t_end=>160,n_steps=>100})'] });
5) Execute the BNGL model file.
[BNGpath]/Perl2/BNG2.pl tlbr.bngl
How to get the hybrid generator
Download and install the latest testing distribution of BioNetGen and NFsim.
Tips: Selecting population Species
Choosing population species is part science and part art. Here are a few rules of thumb.
- The species population should be large throughout the simulation (>10,000). Profiling the model with small scale simulations may be useful for identifying these species. If the population is large at t=0 and shrinks as the the system evolves, the memory savings will be temporal.
- The ideal population species partipates in just a few reaction rules. The hybrid method requires a partial network expansion around population species. If a species participates in too many reaction rules the partial network may become large, resulting in slower simulations. A simple ligand is a good candidate, a receptor with many sites of modification may
Options for generate_hybrid_model
Arguments are passed to generate_hybrid_model in a comma separated list enclosed by curly braces. Each argument/value pair has the form argument=>value.
- verbose=>1 - enable verbose output during hybrid model generation.
- overwrite=>1 - automatically overwrite any existing hybrid model.
- execute=>1 - execute actions after generating the hybrid model.
- suffix=>string - set the model suffix (default='hybrid')
- action=>[] - define a list of actions for the hybrid model.
Hybrid model example: TLBR
The following BNGL model defines a trivalent ligand - bivalent receptor system, selects the free ligand and free receptor species as populations, then generates the hybrid model and simulates the hybrid system using NFsim.
Example: Trivalent ligand, Bivalent receptor model
# Trivalent ligand - bivalent receptor model used to verify and test
# the performance of NFsim. When running this model, you can turn on
# aggregate or complex bookkeeping to block ring formation as in the
# original problem specification by using the parameter "-cb" for
# complex bookkeeping.
begin model
begin parameters
## Peaking Trajectory
Lig_tot 50000
Rec_tot 3000
cTot 2.7
beta 16.8
koff 0.01
kp1 (cTot*koff)/(3.0*Lig_tot) # FREE BINDING RATE
kp2 (beta*koff)/Rec_tot # CROSSLINKING RATE
k_gather 100000
end parameters
begin molecule types
R(l,l)
L(r,r,r)
end molecule types
begin seed species
R(l,l) Rec_tot
L(r,r,r) Lig_tot
end seed species
begin observables
Molecules Rfree R(l,l)
Molecules Lfree L(r,r,r)
Molecules Rtot R()
Molecules Ltot L()
Molecules LRmotif L(r!0).R(l!0)
Species R2 R==2
Species R3 R==3
Species R4 R==4
Species R5 R>4
end observables
begin reaction rules
R(l!1).L(r!1) -> R(l) + L(r) koff
L(r,r,r) + R(l) -> L(r!1,r,r).R(l!1) kp1
L(r,r,r!+) + R(l) -> L(r!1,r,r!+).R(l!1) kp2
L(r,r!+,r!+) + R(l) -> L(r!1,r!+,r!+).R(l!1) kp2
end reaction rules
begin population maps
L(r,r,r) -> L_pop() k_gather
R(l,l) -> R_pop() k_gather
end population maps
end model
## actions ##
generate_hybrid_model({ overwrite=>1, execute=>1, verbose=>1, \
actions=>'simulate_nf({suffix=>nf,complex=>1,\
t_end=>160,n_steps=>100})'] });
Here is the model after hybrid transformation.
# Created by BioNetGen 2.1.8+
substanceUnits("Number");
begin model
begin parameters
Lig_tot 50000
Rec_tot 3000
cTot 2.7
beta 16.8
koff 0.01
kp1 (cTot*koff)/(3.0*Lig_tot)
kp2 (beta*koff)/Rec_tot
k_gather 100000
kp1_1 3*((cTot*koff)/(3.0*Lig_tot))
kp1_2 6*((cTot*koff)/(3.0*Lig_tot))
kp2_1 2*((beta*koff)/Rec_tot)
kp2_2 2*((beta*koff)/Rec_tot)
end parameters
begin molecule types
R_pop() population
L_pop() population
R(l,l)
L(r,r,r)
end molecule types
begin observables
Molecules Rfree R(l,l) R_pop() R_pop()
Molecules Lfree L(r,r,r) L_pop() L_pop() L_pop() L_pop() L_pop() L_pop()
Molecules Rtot R() R_pop()
Molecules Ltot L() L_pop()
Molecules LRmotif L(r!1).R(l!1)
Species R2 R()==2
Species R3 R()==3
Species R4 R()==4
Species R5 R()>4
end observables
begin species
R_pop() Rec_tot
L_pop() Lig_tot
end species
begin reaction rules
Rule1_v1: \
L(r!1).R(l!1) -> R(l) + L(r) koff
Rule2_v1: \
L_pop() + R(l) -> L(r!1,r,r).R(l!1) kp1_1
Rule2_v2: \
L_pop() + R_pop() -> L(r!1,r,r).R(l!1,l) kp1_2
Rule3_v1: \
L(r%T1!+,r%T2,r%T3) + R(l) -> L(r%T1!+,r%T2!1,r%T3).R(l!1) kp2
Rule3_v2: \
L(r%T1!+,r%T2,r%T3) + R_pop() -> L(r%T1!+,r%T2!1,r%T3).R(l!1,l) kp2_1
Rule4_v1: \
L(r%T1!+,r%T2!+,r%T3) + R(l) -> L(r%T1!+,r%T2!+,r%T3!1).R(l!1) kp2
Rule4_v2: \
L(r%T1!+,r%T2!+,r%T3) + R_pop() -> L(r%T1!+,r%T2!+,r%T3!1).R(l!1,l) kp2_2
MapRule0: \
L(r,r,r) -> L_pop() k_gather
MapRule1: \
R(l,l) -> R_pop() k_gather
end reaction rules
end model
### model actions ###
simulate_nf({suffix=>nf,complex=>1,t_end=>160,n_steps=>100});
