`-bscb` appears to be required for reactant-side molecularity enforcement; `-cb` / `complex=>1` is insufficient

[wshlavacek]

I ran into what appears to be an important gap between NFsim flag semantics and the behavior one would expect for reactant-side molecularity in bimolecular BNGL rules.

In short:

• complex=>1 in BNGL maps to NFsim -cb
• -cb turns on complex bookkeeping
• but the actual reactant-side distinct-complex molecularity check appears to require -bscb, not plain -cb

This matters for models where a top-level + on the LHS is intended to mean "these reactant patterns must come from distinct complexes", as in cyclic crosslinking models.

Why this is surprising

For a bimolecular rule such as

A(a) + B(b) -> A(a!1).B(b!1) k

or, more importantly, for cyclic crosslinking rules where an LHS has separate reactant patterns that should not be matched within the same pre-existing connected complex, I would expect reactant-side molecularity to be enforced when complex bookkeeping is enabled.

However, in current NFsim source, plain -cb appears to be insufficient for that.

Source-level evidence

1. complex=>1 maps to -cb

In BioNetGen 2.9.3:

• simulate_nf(..., complex=>1) maps to NFsim -cb
• complex=>1 is already the default for simulate_nf

2. NFsim distinguishes -cb from -bscb

In src/NFsim.cpp, -bscb sets blockSameComplexBinding = true, while -cb separately sets turnOnComplexBookkeeping = true.

3. The actual transformation-set molecularity gate is wired to blockSameComplexBinding

In src/NFinput/NFinput.cpp, when reaction rules are initialized, NFsim does:

ts->setComplexBookkeeping( blockSameComplexBinding );

not

ts->setComplexBookkeeping( system_is_using_complex_bookkeeping );

4. ReactionClass::fire() checks molecularity through TransformationSet::checkMolecularity()

and in TransformationSet::checkMolecularity(), the full distinct-complex check is only taken when that transformation-set complex_bookkeeping flag is true.

Otherwise NFsim falls back to a weaker collision / overlap-style check rather than a full "distinct reactant complexes" check.

Practical consequence

For BNGL/NFsim users, this means:

• complex=>1 / -cb is not enough if correctness depends on reactant-side distinct-complex enforcement
• param=>"-bscb" is required

This is not obvious from the BNGL side, because complex=>1 sounds like it should be the relevant switch.

Minimal reproducer pattern

The problem is especially visible in cyclic crosslinking models, e.g. bivalent ligand / bivalent receptor models where a rule like

L(r!+,r) + R(l) -> L(r!+,r!1).R(l!1) kpx

must not match the R(l) from the same pre-existing complex when the intended chemistry is intermolecular crosslinking rather than an intracomplex closure via a bimolecular reactant match.

Without -bscb, NFsim appears able to admit matches that violate that intended reactant-side molecularity restriction. With -bscb, the resulting dynamics match the corresponding analytical theory much better.

Request

Could you clarify what the intended behavior is?

If the intended behavior is:

• -cb = bookkeeping only
• -bscb = bookkeeping + reactant-side distinct-complex enforcement

then I think this needs much clearer documentation, because it is very easy for BNGL users to assume that complex=>1 is the molecularity-relevant setting.

If instead the intended behavior is that ordinary reactant-side molecularity should already be enforced under -cb, then this may be a bug in how TransformationSet::setComplexBookkeeping() is wired during rule initialization.

I would be happy to provide a small BNGL reproducer if that would be useful.

[jrfaeder]

Yes, it would be helpful to have a simple example for validation. I agree the current behavior is confusing for users and should be modified as well as better documented. Ideally, these flags would be set automatically for the user based on the BNGL model specification.

[wshlavacek]

I put together a reduced validation case for this issue.

It is a one-dose version of the bivalent ligand / bivalent receptor size-2-ring Posner et al. (1995) model:

• 2 molecule types (L, R)
• 5 rules
• 3 observables (Y1, Y2, Z2)
• no Species observables in the reduced reproducer

Local bundle:

• nfsim-issue51-reproducer/issue51_lowdose_reproducer.bngl
• nfsim-issue51-reproducer/validate_issue51.py

Reproduction:

/path/to/BioNetGen-2.9.3/BNG2.pl issue51_lowdose_reproducer.bngl

/path/to/nfsim/build/NFsim \
  -xml issue51_lowdose_reproducer_cb_default.xml \
  -sim 600 -oSteps 240 -seed 1 -ogf \
  -o cb_seed1.gdat

/path/to/nfsim/build/NFsim \
  -xml issue51_lowdose_reproducer_cb_default.xml \
  -sim 600 -oSteps 240 -seed 1 -ogf -bscb \
  -o bscb_seed1.gdat

With the same XML and same RNG seed, I get different final values:

• plain bookkeeping / default simulate_nf behavior: Y1 = 2, Y2 = 3, Z2 = 24
• with -bscb: Y1 = 4, Y2 = 3, Z2 = 18

So the flag changes the trajectory, not just observable bookkeeping.

I also ran a 50-seed ensemble on the same XML (seeds 3001..3050) and compared to the Posner low-dose theory:

• 50-seed mean with plain bookkeeping: Y1 = 3.68, Y2 = 3.46, Z2 = 20.88
• 50-seed mean with -bscb: Y1 = 3.44, Y2 = 2.62, Z2 = 21.2
• Posner ODE at t = 600: Y1 = 3.2087, Y2 = 2.1873, Z2 = 22.0768
• Posner equilibrium: Y1 = 3.1758, Y2 = 2.0689, Z2 = 22.4345

For this reduced case, -bscb materially improves agreement for Y2, which is the observable most directly affected by the second binding step

L(r~c!+,r~u) + R(l~u) -> L(r~c!+,r~c!1).R(l~c!1) kpx

where the top-level + is intended to require the R(l~u) reactant to come from a different complex.

I have not mixed Species observables into this reduced reproducer, so this example is intended to isolate the reactant-side molecularity question rather than any bookkeeping side effects from Species observables.

[wshlavacek]

nfsim-issue51-reproducer.zip