Skip to content

Fix (and new test) for steady state calculations with conserved cycle's changing total mass. #1354

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

Merged
luciansmith merged 1 commit into from
Sep 3, 2026
Merged

Fix (and new test) for steady state calculations with conserved cycle's changing total mass. #1354

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
39 changes: 30 additions & 9 deletions source/llvm/LLVMModelGenerator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -312,25 +312,46 @@ namespace rrllvm {

std::vector<std::string> newSymbols = newModel->getRateRuleSymbols();

// Species whose value is derived from a conserved-moiety total (see
// LLVMExecutableModel::setFloatingSpeciesAmounts) must be transferred
// after every other floating species: setting one of them shifts the
// conserved total by the difference from its currently derived value,
// which is only correct once the other species in that conservation
// law hold their transferred values rather than the new model's
// freshly-generated init values.
auto transferFloatingSpecies = [&](int i, const std::string& id) {
double value = 0;
oldModel->getFloatingSpeciesAmounts(1, &i, &value);
try {
newModel->setValue(id, value);
}
catch (const exception& e) {
rrLog(Logger::LOG_WARNING) << "regenerateModel: failed to transfer current "
"value for floating species '" << id << "' (value " << value
<< "): " << e.what();
}
};

std::vector<int> deferredConservedMoietySpecies;

for (int i = 0; i < oldModel->getNumFloatingSpecies(); i++) {
std::string id = oldModel->getFloatingSpeciesId(i);
int index = newModel->getFloatingSpeciesIndex(id);

if (index >= 0) {
// new model has this species
double value = 0;
oldModel->getFloatingSpeciesAmounts(1, &i, &value);
try {
newModel->setValue(id, value);
}
catch (const exception& e) {
rrLog(Logger::LOG_WARNING) << "regenerateModel: failed to transfer current "
"value for floating species '" << id << "' (value " << value
<< "): " << e.what();
if (newModel->symbols->isConservedMoietySpecies(id)) {
deferredConservedMoietySpecies.push_back(i);
continue;
}
transferFloatingSpecies(i, id);
}
}

for (int i : deferredConservedMoietySpecies) {
transferFloatingSpecies(i, oldModel->getFloatingSpeciesId(i));
}


for (int i = 0; i < oldModel->getNumBoundarySpecies(); i++) {
std::string id = oldModel->getBoundarySpeciesId(i);
Expand Down
23 changes: 23 additions & 0 deletions test/model_analysis/model_analysis.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,29 @@ class ModelAnalysisTests : public RoadRunnerTest {
};


TEST_F(ModelAnalysisTests, issue1351_steadyStateAfterSimulate) {
// https://github.com/sys-bio/roadrunner/issues/1351
// This model (ATP <-> ADP <-> AMP) has one conserved moiety, so
// ATP + ADP + AMP must equal the same total (2.1 + 1.5 + 0.33 = 3.93)
// whether it's read right after construction, after simulate(), or
// after steadyState() runs following a simulate().
const double conservedTotal = 2.1 + 1.5 + 0.33;

RoadRunner ssOnly((modelAnalysisModelsDir / "steady_state_after_sim.xml").string());
ssOnly.steadyState();
EXPECT_NEAR(ssOnly.getValue("ATP") + ssOnly.getValue("ADP") + ssOnly.getValue("AMP"),
conservedTotal, 1e-6);

RoadRunner rr((modelAnalysisModelsDir / "steady_state_after_sim.xml").string());
rr.simulate(0, 100, 200);
EXPECT_NEAR(rr.getValue("ATP") + rr.getValue("ADP") + rr.getValue("AMP"),
conservedTotal, 1e-6);

rr.steadyState();
EXPECT_NEAR(rr.getValue("ATP") + rr.getValue("ADP") + rr.getValue("AMP"),
conservedTotal, 1e-6);
}

// https://github.com/sys-bio/roadrunner/issues/1339
// SimulateOptions is a long-lived object (RoadRunner::self.simulateOpt) that
// the C API mutates field-by-field across repeated simulate() calls. The
Expand Down
87 changes: 87 additions & 0 deletions test/models/ModelAnalysis/steady_state_after_sim.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,87 @@
<?xml version="1.0" encoding="UTF-8"?>
<!-- Created by libAntimony version v3.2.0 with libSBML version 5.21.1. -->
<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
<model metaid="__main" id="__main">
<listOfCompartments>
<compartment sboTerm="SBO:0000410" id="default_compartment" spatialDimensions="3" size="1" constant="true"/>
</listOfCompartments>
<listOfSpecies>
<species id="ATP" compartment="default_compartment" initialConcentration="2.1" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species id="AMP" compartment="default_compartment" initialConcentration="0.33" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
<species id="ADP" compartment="default_compartment" initialConcentration="1.5" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
</listOfSpecies>
<listOfParameters>
<parameter id="kf" value="432.9" constant="true"/>
<parameter id="kr" value="133.333" constant="true"/>
<parameter id="kc" value="3.2076" constant="true"/>
<parameter id="kg" value="5" constant="true"/>
</listOfParameters>
<listOfReactions>
<reaction id="vAK" reversible="true">
<listOfReactants>
<speciesReference species="ATP" stoichiometry="1" constant="true"/>
<speciesReference species="AMP" stoichiometry="1" constant="true"/>
</listOfReactants>
<listOfProducts>
<speciesReference species="ADP" stoichiometry="2" constant="true"/>
</listOfProducts>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<minus/>
<apply>
<times/>
<ci> kf </ci>
<ci> AMP </ci>
<ci> ATP </ci>
</apply>
<apply>
<times/>
<ci> kr </ci>
<apply>
<power/>
<ci> ADP </ci>
<cn type="integer"> 2 </cn>
</apply>
</apply>
</apply>
</math>
</kineticLaw>
</reaction>
<reaction id="vconsum" reversible="true">
<listOfReactants>
<speciesReference species="ATP" stoichiometry="1" constant="true"/>
</listOfReactants>
<listOfProducts>
<speciesReference species="ADP" stoichiometry="1" constant="true"/>
</listOfProducts>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<times/>
<ci> kc </ci>
<ci> ATP </ci>
</apply>
</math>
</kineticLaw>
</reaction>
<reaction id="vgen" reversible="true">
<listOfReactants>
<speciesReference species="ADP" stoichiometry="1" constant="true"/>
</listOfReactants>
<listOfProducts>
<speciesReference species="ATP" stoichiometry="1" constant="true"/>
</listOfProducts>
<kineticLaw>
<math xmlns="http://www.w3.org/1998/Math/MathML">
<apply>
<times/>
<ci> kg </ci>
<ci> ADP </ci>
</apply>
</math>
</kineticLaw>
</reaction>
</listOfReactions>
</model>
</sbml>
Loading