Rules that use TotalRate are scaled by the reaction center symmetry factor

[wshlavacek]

NFsim applies the reaction center symmetry factor to rules that use the TotalRate
keyword. Those rules run at a fraction of the rate the model asks for. For a
homodimer the fraction is one half.

TotalRate means the rate law gives the whole propensity of the rule. The symmetry
factor exists to correct a counting problem, namely that a symmetric pattern
matches the same reaction more than once. Under TotalRate there is no counting to
correct, because the rate is stated outright, so the factor should not be applied.

Repro

begin model
begin parameters
  k 0.02
end parameters
begin molecule types
  A(b)
end molecule types
begin species
  A(b) 400
end species
begin observables
  Molecules Afree A(b)
end observables
begin reaction rules
  A(b) + A(b) -> A(b!1).A(b!1)  k TotalRate
end reaction rules
end model
begin actions
writeXML()
end actions

BioNetGen 2.9.3 writes this rule with symmetry_factor="0.5" and totalrate="1".

The propensity should be k, which is 0.02, so about 20 reactions fire over 1000
time units. Each one consumes two A, so Afree should end near 360.

NFsim's own verbose output shows the problem on a single line:

ReactionClass: _R1  ( baseFunction=_rateLaw1()=0.02,  a=0.01, fired=0 times, TotalRate=on )

The rate law evaluates to 0.02 and the propensity comes out as 0.01.

Mean of Afree at t=1000 over 20 seeds, same input file both times:

master (5962ea9)   383.4
v1.14.3            360.3

Where it happens

The constructor folds the symmetry factor into baseRate at
src/NFcore/reactionClass.cpp:191. FunctionalRxnClass::update_a then multiplies by
that member at src/NFreactions/reactions/reaction.cpp:87, before the branch that
decides whether the reactant counts apply.

a *= this->baseRate;
...
if(!this->totalRateFlag) {
    for(unsigned int i=0; i<n_reactants; i++)
        a*=(double)getCorrectedReactantCount(i);
}

So the factor lands on the rate whether or not the rule states a total rate.

This came in with PR #89, which made the symmetry factor reach rate laws that do
not go through setBaseRate(). That change was correct for ordinary rate laws. The
TotalRate case looks like it was simply missed, because the same member is shared.

Scope

In practice only FunctionalRxnClass is affected. BioNetGen forces every TotalRate
rate law into a function, even when the rate is a bare constant (RateLaw.pm line
84 sets a force function flag when TotalRate is set), so a TotalRate rule never
reaches the elementary path where setBaseRate applies the factor. BioNetGen also
rejects TotalRate on Sat, MM and Hill rate laws, on Arrhenius, and on local
functions, so the other reaction classes cannot be reached with a total rate from
generated input.

Suggested fix

Skip the symmetry factor when totalRateFlag is set. The check belongs where the
propensity is computed. A check placed where the factor is folded into baseRate
will not fire, because NFinput calls setTotalRateFlag() well after both the
constructor and setBaseRate(), so the flag is still false at both of those points.

On the reading of TotalRate

Two other implementations agree that the symmetry factor should not apply under
TotalRate. BioNetGen's own C++ engine returns the bare rate constant for a
TotalRate reaction and applies its statistical factor only in the other branch.
RuleMonkey does the same. The field comment for the flag in BioNetGen's RxnRule.pm
describes it as a TotalRate as against a per site rate, and the symmetry factor is
named MultScale there, described as the factor by which to scale multiplicity.
Under TotalRate no multiplicity enters the calculation, so there is nothing for it
to scale.