diff --git a/src/NFreactions/reactions/reaction.cpp b/src/NFreactions/reactions/reaction.cpp index 141a2dc0..48814d8a 100644 --- a/src/NFreactions/reactions/reaction.cpp +++ b/src/NFreactions/reactions/reaction.cpp @@ -84,7 +84,17 @@ double FunctionalRxnClass::update_a() { // through setBaseRate(), so this factor is the reaction center symmetry // correction and nothing else. Without it a symmetric rule with a // functional rate fires at 1/symmetryFactor times its intended rate. - a *= this->baseRate; + // + // Not under TotalRate, though. The symmetry factor corrects a counting + // problem -- a reactant pattern with a non-trivial automorphism matches the + // same reaction more than once -- and TotalRate states the whole propensity + // of the rule outright, so there is no count to correct. It has to be + // skipped here rather than where the factor is folded into baseRate: + // NFinput calls setTotalRateFlag() well after both this class' constructor + // and setBaseRate(), so totalRateFlag is still false at both of those + // points and a guard placed there would never fire. + if (!this->totalRateFlag) + { a *= this->baseRate; } if(a<0) { cout<<"Warning!! The function you provided for functional rxn: '"< symmetry_factor="0.5" +# Tasym asymmetric TotalRate, rate 2*kt -> symmetry_factor="1" (control) +# +# Under TotalRate the propensity is constant while both reactant lists are +# non-empty, so consumption is exactly linear in time and the expected survivor +# count is closed-form rather than exponential: +# +# Tsym fires kt*T times, consuming two A per firing -> X0 - 2*kt*T +# Tasym fires 2*kt*T times, consuming one B per firing -> X0 - 2*kt*T +# +# With X0=4000, kt=1.0 and T=1000 both land on 2000 survivors. A wrongly applied +# 0.5 halves Tsym's propensity to 500 firings and leaves 3000 -- and leaves the +# asymmetric control untouched, because its symmetry factor is 1. Both pools end +# at 2000 free molecules, so neither reactant list ever empties and the constant +# propensity holds for the whole run. +# +# This fixture deliberately covers only the TotalRate row. The non-TotalRate row +# -- the symmetry factor that must still be applied on every rate law -- is +# symmetry_factor_rate_laws.bngl, which is left untouched so that a fix which +# simply stopped applying the factor everywhere fails there rather than passing +# here. + +begin model +begin parameters + kt 1.0 + kt2 2.0 # == 2*kt: one B consumed per firing instead of two A + X0 4000 +end parameters + +begin molecule types + A(b) + B(b) + C(b) +end molecule types + +begin seed species + A(b) X0 + B(b) X0 + C(b) X0 +end seed species + +begin observables + Molecules Tsym_free A(b) + Molecules Tasym_free B(b) +end observables + +begin reaction rules + # Symmetric reaction center: BNG emits symmetry_factor="0.5" here. + Tsym: A(b) + A(b) -> A(b!1).A(b!1) kt TotalRate + # Asymmetric control: symmetry_factor="1", so no fix can change it. + Tasym: B(b) + C(b) -> B(b!1).C(b!1) kt2 TotalRate +end reaction rules +end model + +begin actions +writeXML() +end actions diff --git a/test/symmetry/symmetry_factor_total_rate.xml b/test/symmetry/symmetry_factor_total_rate.xml new file mode 100644 index 00000000..b1028cdd --- /dev/null +++ b/test/symmetry/symmetry_factor_total_rate.xml @@ -0,0 +1,211 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + kt + + + + + + kt2 + + + + diff --git a/validate/validate.py b/validate/validate.py index 1fc76b96..081431ce 100644 --- a/validate/validate.py +++ b/validate/validate.py @@ -577,6 +577,73 @@ def test_symmetry_factor_is_applied_on_every_rate_law(self): ), ) + def test_totalrate_rules_are_not_scaled_by_the_symmetry_factor(self): + # FunctionalRxnClass::update_a() scales the propensity by baseRate, which + # for that class is the reaction center symmetry factor and nothing else. + # It applied it whether or not the rule uses TotalRate, so a TotalRate + # rule with a symmetric reaction center ran at a fraction of the rate the + # model asks for, one half on a homodimer. + # + # TotalRate means the rate law gives the whole propensity of the rule. The + # symmetry factor corrects a match count, and under TotalRate there is no + # count to correct, so the factor must not be applied. + # + # Every BNG-generated TotalRate rule reaches this class and no other. + # BNG2.pl forces a TotalRate rate law into a Function even when it is a + # bare constant, which routes it past the Ele/setBaseRate() path, and it + # rejects TotalRate on Sat/MM/Hill, on Arrhenius, and on local functions. + # + # Under TotalRate the propensity is constant while both reactant lists are + # non-empty, so consumption is linear in time and the expected survivor + # count is closed form. Both pools are built to land on the same one: + # Tsym fires kt*T times consuming two A each, Tasym fires 2*kt*T times + # consuming one B each, so both leave X0 - 2*kt*T = 4000 - 2000 = 2000. + # A wrongly applied 0.5 halves Tsym's propensity and leaves 3000 instead. + xmlPath = os.path.join( + nfsimPrePath, "test", "symmetry", "symmetry_factor_total_rate.xml" + ) + headers, mean = self._mean_final_row(xmlPath, "-sim 1000 -oSteps 2", 5) + + expected = 2000.0 + halved = 3000.0 + # Firing counts are Poisson, so the symmetric pool's per-seed scatter is + # 2*sqrt(1000) ~= 63 counts and a 5-seed mean has sigma ~= 28. A +/-200 + # band is ~7 sigma wide and still leaves the halved value 800 counts + # outside it. + tolerance = 200.0 + + pools = [ + ("Tsym_free", "symmetric reaction center, TotalRate"), + # symmetry_factor is 1 here, so no placement of the factor can move + # this pool. It fails only if TotalRate handling broke some other way. + ("Tasym_free", "asymmetric reaction center, TotalRate (control)"), + ] + for name, description in pools: + got = mean[headers.index(name)] + self.assertAlmostEqual( + got, + expected, + delta=tolerance, + msg="{0} ({1}) ended at {2:.1f}, expected about {3:.1f}; " + "{4:.1f} means the symmetry factor was applied to a rule that " + "states its own total rate".format( + name, description, got, expected, halved + ), + ) + + # The sharpest form of the claim, needing no expected value at all: under + # TotalRate a rule's propensity is whatever it says it is, so it must not + # depend on whether its reactant pattern happens to be symmetric. + sym = mean[headers.index("Tsym_free")] + asym = mean[headers.index("Tasym_free")] + self.assertAlmostEqual( + sym, + asym, + delta=2 * tolerance, + msg="TotalRate: symmetric pool ended at {0:.1f}, asymmetric at " + "{1:.1f}".format(sym, asym), + ) + def test_michaelis_menten_symmetry_factor_scales_the_substrate_count(self): # Where the MM law is linear in the substrate match count, scaling that # count and scaling the finished propensity coincide, so the fixture above