BioNetGen Language Support¶
MCell4 uses BioNetGen rule-based definition formalism to define species and reaction rules. Documentation on the BioNetGen language (BNGL) can be found on the official BioNetGen website. This section contains a short overview of BNGL with described limitations and extensions.
Loading a BNGL File¶
A file in BioNetGen language can be loaded directly from MCell4 Python model without any previous conversions like this:
model = m.Model()
...
# load our BNGL file and set directory for BNGL observables
model.load_bngl('model.bngl', './react_data/seed_' + str(SEED).zfill(5) + '/')
Calling load_bngl invokes a BNGL parser and the parsed data are converted internally into MCell4 classes. For more details see Model.load_bngl (that loads the full BNGL file), or methods that load only partial information: Subsystem.load_bngl_molecule_types_and_reaction_rules, Instatiation.load_bngl_compartments_and_seed_species, Observables.load_bngl_observables, and bngl_utils.load_bngl_parameters.
Terminology¶
A component type is a template for components that defines name and all allowed states of a component. An elementary molecule type is a template for elementary molecules that defines name and all component types used by an elementary molecule.
Example in BNGL:
begin molecule types
# simple elementary molecule type 'A' that has no components
A
# elementary molecule type 'B' that uses a single component type 'c' with no states
B(c)
# elementary molecule type 'C' that uses a single component type 'b' with states 'A' and 'B'
C(b~A~B)
end molecule types
A BNGL complex is defined as a set of elementary molecules where each elementary molecule may contain components. These components may have state and bond declared. Through bonds, the elementary molecules are linked together.
A complex is either a pattern where not all components need to be specified or a fully-specified complex where all contained elementary molecules contain all components for their type, and their states are set.
# block seed species defines initial molecules in a simulation
# all complexes used here must be fully-qualified
begin seed species
# complex with a single elementary molecule A, the initial count is 10
A 10
# complex with a single elementary molecule C with component b state set to A
C(b~A)
# complex composed of elem. molecules B and C bound through components
B(c!1).C(b~A!1)
end molecule types
# block observables defines what patterns will be counted,
# the complexes used here do not have to be fully-qualified
begin observables
# pattern with a single elementary molecule A, name of the observable is obs_A
Molecules obs_A A
# pattern matching a molecule C with any component state, bound or unbound
Molecules obs_C C
# pattern matching a molecule C with component's b state A, component b must be unbound
Molecules obs_CbA C(b~A)
# pattern matching a molecule that has both B and C with any component states
# (in our example, the only complexes that match must be bound but the pattern does not require it)
Molecules obs_CB B.C
end observables
What we call elementary molecule is in BNGL usually called just molecule, but we needed to be precise because in the MCell context, the term molecule is a single simulated agent. Each MCell molecule has a species assigned that defines its type. A species is represented by a fully-specified complex.
Reaction Rules¶
Simple Reaction Rules¶
We will start with a simple example of a reversible bimolecular reaction.
begin reaction rules
A + B <-> C fwd_rate, rev_rate
end reaction rules
All reaction rules in MCell follow the law of mass action and this particular rule tells us that when molecule of species A hits a molecule of species B, then with a probability given by fwd_rate, a reaction can occur that transforms A and B into C. This is a reversible reaction rule, so also a unimolecular reaction is defined for C to split into A and B and the probability of this happening is given by rev_rate.
Graph Transformations¶
A reaction rule in BNGL defines a graph transformation that creates or removes bonds, changes states, or adds or removes elementary molecules. We will demonstrate this graph transformation on an example:
begin molecule types
# A has two identical component types
A(c0~R~S, c0~R~S)
# B has two different component types (although their states are named identically)
B(c2~U~V, c3~U~V)
end molecule types
begin reaction rules
A(c0~R) + B(c2) -> A(c0~S!1).B(c2!1) rate
end reaction rules
The reaction rule makes a bond between A’s component c0 with B’s c2 and also changes the state of A’s c0 from R to S as shown in the following image:
The way this transformation works is shown in the following image. Our inputs are reactants (2nd row). A mapping from each elementary molecules and each component from reactant patterns onto reactants is computed. All elementary molecules and components must be mapped otherwise the rule cannot be applied. There can be multiple mappings in general, but there is just one in our example.
The next step is to compute a mapping of elementary molecules and components from the reaction rule product pattern onto reactant patterns.
The difference between the reaction rule product pattern and the reactant patterns tells us what changes need to be made. In this instance, a bond between A’s c0 with state R and B’s c2 is created,d the same A’s component’s c0 state is changed to S.
Once we have our mappings, we can follow the arrows leading from the reaction rule product pattern to reactant patterns and then to reactants and do exactly these changes on the reactants, as shown in the row reactants changed into products.
Reaction Rates¶
MCell and BioNetGen use different units for bimolecular kinetic rates. In MCell, a volume-volume reaction (reaction between two molecules that are free to diffuse in 3D space) is M-1*s-1where M is the molar concentration (number of moles per liter). In BioNetGen, the user is not constricted to a specific unit but a usual unit is N-1*s-1where N is a number of molecules per compartment and the default compartment volume if a compartment is not specified is 1fl (= 1 um3).
To convert from BioNetGen to MCell units, one needs to multiply the BioNetGen rate by NA * V where NA is Avogardro’s constant and V is volume of the compartment is liters as derived here:
1/M = 1/(#moles/V) = 1/( (N/NA)/V) = NA * V * 1/N
The unimolecular reaction rates in MCell and BioNetGen both use unit s-1.
MCell Parameter Prefixes¶
To simplify interoperability between BioNetGen and MCell, special parameters with reserved prefixes and names are used to provide the information needed by MCell such as diffusion constants.
MCELL_REDEFINE_ - allows a parameter to have a different value in MCell and for BioNetGen. Example:
begin parameters
NA_um3 6.022e8
VOL_RXN 1
# we will redefine the parameter VOL_RXN for MCell
MCELL_REDEFINE_VOL_RXN NA_um3
# the BioNetGen value of the parameter kp1 is 1 but for MCell
# it is 6.022e8 * 1
kp1 VOL_RXN * 1
end parameters
MCELL_DIFFUSION_CONSTANT_3D_, MCELL_DIFFUSION_CONSTANT_2D_ - allows to specify a diffusion constant for a 3D (volume) or 2D (surface) elementary molecule type.
begin parameters
MCELL_DIFFUSION_CONSTANT_3D_V 1e-6
MCELL_DIFFUSION_CONSTANT_2D_S 1e-7
end parameters
begin molecule types
V
S
end molecule types
All these parameter prefixes are ignored by the BioNetGen tools and are interpreted as any other parameter.
Compartments and Orientations¶
BNGL compartments allow to define hierarchical volumes or surfaces where molecules are located. With compartments, one can define e.g. a transport reaction where a molecule that usually diffuses in volume is transported by a channel located in a membrane into another volume compartment. Let’s say we need to define a reaction where a volume molecule A in a volume/3D compartment reacts with a surface molecule T (transporter) in a surface compartment as shown in Fig. X7.
An example of compartments is shown in the following figure. EC is extracellular space, PM is plasma membrane, and CP is cytoplasm. A is a molecule that diffuses freely in 3D space, and T is located in a membrane.
In BNGL, a reaction that defines transport of A from compartment EC into CP is represented like this: A@EC + T@PM -> A@CP + T@PM. An issue arises when one needs to model multiple instances of cells or organelles. A compartment is a specific object in the simulation. If we wanted to simulate multiple cells, we would need to repeat the definition of this reaction rule for each cell (with PM1, PM2, …) e.g. like this: A@EC + T@PM1 -> A@CP1 + T@PM1, A@EC + T@PM2 -> A@CP2 + T@PM2, etc.
To avoid the repetition of reaction rules for each compartment and to keep the BNG language consistent, an extension to the BNG language that uses compartment classes @IN and @OUT was introduced. The original BNG reaction with specific compartments is more generally represented as A@OUT + T -> A@IN + T. Only bimolecular reactions with one volume a one surface reactant may use @IN or @OUT compartment classes because the compartment of the surface reactant defines the meaning of the @IN and @OUT compartment class of the volume reactant.
When this rule is applied to reactants A@EC and T@CP, we know that the compartment of T is PM, the compartment outside (parent) is EC, and inside is CP. So, we insert this information to the rule A@OUT + T -> A@IN + T and get A@EC + T@PM -> A@CP + T@PM (same as the example rule we started with).
One more related case is the position of the reaction product. Let’s have a reaction AT -> T + A. AT is a surface molecule that represents A bound to T. When it dissociates, A must be placed to the ‘inside’ or ‘outside’ compartment. When no compartment or compartment class is specified for the product volume molecule, the target compartment is selected randomly. A surface product molecule stays at the same surface compartment as the reactant.
MCell3 used molecule orientations to limit which surface and volume reactions can occur. A surface molecule S could be released or created as a product in two different states S’ (UP) and S, (DOWN). This behavior was deprecated in MCell4 and all surface molecules are by default created with orientation UP. If specific orientation is needed, one can use a component with state such as S(orient~UP~DOWN).
Limitations¶
This section lists the limitations of the MCell4 BNGL support compared to the BioNetGen tools.
Section functions is not supported, and loading a file with BNGL functions will report an error.
The only reaction kinetics supported is the law of mass action, BNG supports others through keywords but the MCell4 parser will report an error.
Actions such as generate_network or simulate are silently ignored.