Skip to content

Do not apply the reaction center symmetry factor under TotalRate #92

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wshlavacek wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Do not apply the reaction center symmetry factor under TotalRate #92

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 11 additions & 1 deletion src/NFreactions/reactions/reaction.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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: '"<<name<<"' evaluates\n";
Expand Down
72 changes: 72 additions & 0 deletions test/symmetry/symmetry_factor_total_rate.bngl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
# The reaction center symmetry factor must not be applied under TotalRate.
#
# TotalRate states the whole propensity of a rule outright. The symmetry factor
# corrects a counting problem -- a reactant pattern with a non-trivial
# automorphism matches the same reaction more than once -- so where nothing is
# being counted there is nothing to correct, and the factor must not be applied.
#
# BNG2.pl forces every TotalRate rate law into a Function, even a bare constant
# (RateLaw.pm: `my $force_fcn = $totalRate ? 1 : 0;`), and rejects TotalRate on
# Sat/MM/Hill, Arrhenius, and local functions. FunctionalRxnClass is therefore
# the only reaction class a BNG-generated TotalRate rule can reach, and it is
# also the class that "Apply the reaction center symmetry factor on every rate
# law, not just Ele" taught to scale by baseRate -- which for that class carries
# the symmetry factor and nothing else. That is how the two features collided.
#
# Tsym symmetric TotalRate, rate kt -> 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
211 changes: 211 additions & 0 deletions test/symmetry/symmetry_factor_total_rate.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,211 @@
<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by BioNetGen 2.9.3 -->
<sbml xmlns="http://www.sbml.org/sbml/level3" level="3" version="1">
<model id="symmetry_factor_total_rate">
<ListOfParameters>
<Parameter id="kt" type="Constant" value="1" expr="1.0"/>
<Parameter id="kt2" type="Constant" value="2" expr="2.0"/>
<Parameter id="X0" type="Constant" value="4000" expr="4000"/>
</ListOfParameters>
<ListOfMoleculeTypes>
<MoleculeType id="A">
<ListOfComponentTypes>
<ComponentType id="b"/>
</ListOfComponentTypes>
</MoleculeType>
<MoleculeType id="B">
<ListOfComponentTypes>
<ComponentType id="b"/>
</ListOfComponentTypes>
</MoleculeType>
<MoleculeType id="C">
<ListOfComponentTypes>
<ComponentType id="b"/>
</ListOfComponentTypes>
</MoleculeType>
</ListOfMoleculeTypes>
<ListOfCompartments>
</ListOfCompartments>
<ListOfSpecies>
<Species id="S1" concentration="X0" name="A(b)">
<ListOfMolecules>
<Molecule id="S1_M1" name="A">
<ListOfComponents>
<Component id="S1_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</Species>
<Species id="S2" concentration="X0" name="B(b)">
<ListOfMolecules>
<Molecule id="S2_M1" name="B">
<ListOfComponents>
<Component id="S2_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</Species>
<Species id="S3" concentration="X0" name="C(b)">
<ListOfMolecules>
<Molecule id="S3_M1" name="C">
<ListOfComponents>
<Component id="S3_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</Species>
</ListOfSpecies>
<ListOfReactionRules>
<ReactionRule id="RR1" name="Tsym" symmetry_factor="0.5">
<ListOfReactantPatterns>
<ReactantPattern id="RR1_RP1">
<ListOfMolecules>
<Molecule id="RR1_RP1_M1" name="A">
<ListOfComponents>
<Component id="RR1_RP1_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</ReactantPattern>
<ReactantPattern id="RR1_RP2">
<ListOfMolecules>
<Molecule id="RR1_RP2_M1" name="A">
<ListOfComponents>
<Component id="RR1_RP2_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</ReactantPattern>
</ListOfReactantPatterns>
<ListOfProductPatterns>
<ProductPattern id="RR1_PP1">
<ListOfMolecules>
<Molecule id="RR1_PP1_M1" name="A">
<ListOfComponents>
<Component id="RR1_PP1_M1_C1" name="b" numberOfBonds="1"/>
</ListOfComponents>
</Molecule>
<Molecule id="RR1_PP1_M2" name="A">
<ListOfComponents>
<Component id="RR1_PP1_M2_C1" name="b" numberOfBonds="1"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
<ListOfBonds>
<Bond id="RR1_PP1_B1" site1="RR1_PP1_M1_C1" site2="RR1_PP1_M2_C1"/>
</ListOfBonds>
</ProductPattern>
</ListOfProductPatterns>
<RateLaw id="RR1_RateLaw" type="Function" name="_rateLaw1" totalrate="1">
<ListOfArguments>
</ListOfArguments>
</RateLaw>
<Map>
<MapItem sourceID="RR1_RP1_M1" targetID="RR1_PP1_M1"/>
<MapItem sourceID="RR1_RP1_M1_C1" targetID="RR1_PP1_M1_C1"/>
<MapItem sourceID="RR1_RP2_M1" targetID="RR1_PP1_M2"/>
<MapItem sourceID="RR1_RP2_M1_C1" targetID="RR1_PP1_M2_C1"/>
</Map>
<ListOfOperations>
<AddBond site1="RR1_RP1_M1_C1" site2="RR1_RP2_M1_C1"/>
</ListOfOperations>
</ReactionRule>
<ReactionRule id="RR2" name="Tasym" symmetry_factor="1">
<ListOfReactantPatterns>
<ReactantPattern id="RR2_RP1">
<ListOfMolecules>
<Molecule id="RR2_RP1_M1" name="B">
<ListOfComponents>
<Component id="RR2_RP1_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</ReactantPattern>
<ReactantPattern id="RR2_RP2">
<ListOfMolecules>
<Molecule id="RR2_RP2_M1" name="C">
<ListOfComponents>
<Component id="RR2_RP2_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</ReactantPattern>
</ListOfReactantPatterns>
<ListOfProductPatterns>
<ProductPattern id="RR2_PP1">
<ListOfMolecules>
<Molecule id="RR2_PP1_M1" name="B">
<ListOfComponents>
<Component id="RR2_PP1_M1_C1" name="b" numberOfBonds="1"/>
</ListOfComponents>
</Molecule>
<Molecule id="RR2_PP1_M2" name="C">
<ListOfComponents>
<Component id="RR2_PP1_M2_C1" name="b" numberOfBonds="1"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
<ListOfBonds>
<Bond id="RR2_PP1_B1" site1="RR2_PP1_M1_C1" site2="RR2_PP1_M2_C1"/>
</ListOfBonds>
</ProductPattern>
</ListOfProductPatterns>
<RateLaw id="RR2_RateLaw" type="Function" name="_rateLaw2" totalrate="1">
<ListOfArguments>
</ListOfArguments>
</RateLaw>
<Map>
<MapItem sourceID="RR2_RP1_M1" targetID="RR2_PP1_M1"/>
<MapItem sourceID="RR2_RP1_M1_C1" targetID="RR2_PP1_M1_C1"/>
<MapItem sourceID="RR2_RP2_M1" targetID="RR2_PP1_M2"/>
<MapItem sourceID="RR2_RP2_M1_C1" targetID="RR2_PP1_M2_C1"/>
</Map>
<ListOfOperations>
<AddBond site1="RR2_RP1_M1_C1" site2="RR2_RP2_M1_C1"/>
</ListOfOperations>
</ReactionRule>
</ListOfReactionRules>
<ListOfObservables>
<Observable id="O1" name="Tsym_free" type="Molecules">
<ListOfPatterns>
<Pattern id="O1_P1">
<ListOfMolecules>
<Molecule id="O1_P1_M1" name="A">
<ListOfComponents>
<Component id="O1_P1_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</Pattern>
</ListOfPatterns>
</Observable>
<Observable id="O2" name="Tasym_free" type="Molecules">
<ListOfPatterns>
<Pattern id="O2_P1">
<ListOfMolecules>
<Molecule id="O2_P1_M1" name="B">
<ListOfComponents>
<Component id="O2_P1_M1_C1" name="b" numberOfBonds="0"/>
</ListOfComponents>
</Molecule>
</ListOfMolecules>
</Pattern>
</ListOfPatterns>
</Observable>
</ListOfObservables>
<ListOfFunctions>
<Function id="_rateLaw1">
<ListOfReferences>
<Reference name="kt" type="Constant"/>
</ListOfReferences>
<Expression> kt </Expression>
</Function>
<Function id="_rateLaw2">
<ListOfReferences>
<Reference name="kt2" type="Constant"/>
</ListOfReferences>
<Expression> kt2 </Expression>
</Function>
</ListOfFunctions>
</model>
</sbml>
67 changes: 67 additions & 0 deletions validate/validate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading