/** @page libsbml-example Tutorial: Creating a complete model, from start to finish This section presents a complete example program in C++ demonstrating how to create a complete SBML model programmatically and then writing the result to a file. @section characteristics General characteristics of the model to be created The model presented here is the one described in Section 7.1 of the SBML Level 2 Version 4 specification. It consists of four species identified simply as E, S, ES, and P; one compartment identified as cytosol; and two reactions (one reversible, one irreversible) described by the following equation: @image html enzymatic-reaction.jpg @image latex enzymatic-reaction.jpg From these equations, it becomes clear that the model also contains three parameters: kon, koff, and kcat. We will need to assign initial values to these parameters, as well as to the size of the compartment and the initial quantities of the species. Since this model is only meant to illustrate the principles of programming with libSBML, and does not represent an actual biological entity, we hope readers will forgive us for not worrying about justifying the parameter values on the basis of biological plausibility. The table of values for all the quantities we will use is shown below:
moles
Entity in the %Model Initial value Units
cytosol\f$1 \times 10^{-14}\f$litres
S\f$1 \times 10^{-20}\f$moles
P0moles
E\f$5 \times 10^{-21}\f$moles
ES0
kon1 000 000litre/mole/second
koff0.21/second
kcat1/second
For variety, we will use substance units instead of concentration units for the initial values of the species. This is reflected in the use of moles instead of moles/litre as the species quantity units in the table above. @section code The program for creating the model To begin, our program has to include certain header files needed by libSBML under C++. For convenience, we will also set the default C++ namespace to std so that we do not have to prefix common C++ functions used throughout our code. @verbatim /** * File: example.cpp * Description: Create an example model. * Author: Akiya Jouraku */ #include #include using namespace std; @endverbatim The next few lines of our program will not be used until the end of the program, but the relevant definitions are best put near the beginning as a principle of good programming. Specifically, this is the definition of a program name and version number, to be inserted by libSBML as XML comments into the final SBML file. We will do something very simple here and define the relevant values as global parameters; other approaches are possible. @verbatim // // These variables are used in writeExampleSBML when writing an SBML // document. They are handed to libSBML functions in order to include // the program information into comments within the SBML file. // const static string ProgramName = "createExample"; const static string ProgramVersion = "1.0.0"; @endverbatim We need a few forward declarations for helper functions. Here they are: @verbatim // // Forward declarations for functions defined below. // SBMLDocument* createExampleEnzymaticReaction(); bool validateExampleSBML(SBMLDocument *sbmlDoc); bool writeExampleSBML(const SBMLDocument *sbmlDoc, const string& filename); @endverbatim Now we can define a main routine for the program. This essentially acts as scaffolding for calling worker functions and dealing with overall errors. @verbatim // // Main routine // int main (int argc, char *argv[]) { SBMLDocument* sbmlDoc = 0; bool SBMLok = false; try { sbmlDoc = createExampleEnzymaticReaction(); SBMLok = validateExampleSBML(sbmlDoc); if (SBMLok) writeExampleSBML(sbmlDoc, "enzymaticreaction.xml"); // Avoid memory leaks. delete sbmlDoc; // Return non-zero to indicate if a problem occurred. if (!SBMLok) return 1; } catch (std::bad_alloc& e) { cerr << e.what() << ": Unable to allocate memory." << endl; return 1; } catch (...) { cerr << "Unexpected exceptional condition encountered." << endl; return 1; } // A 0 return status is the standard Unix/Linux way to say "all ok". return 0; } @endverbatim The code above has three key functions, createExampleEnzymaticReaction(), validateExampleSBML, and writeExampleSBML(), which do all the real work in the program. Their verbose names hopefully imply their functions. They are described and defined in the next three subsections. @subsection creator Function for creating the model We begin this function by creating an SBMLDocument object, creating a Model object within the document object, and setting the Model object's identifier. @verbatim /** * Create an SBML model represented in "7.1 A Simple example application * of SBML" in the SBML Level 2 Version 4 Specification. */ SBMLDocument* createExampleEnzymaticReaction() { // Create an SBMLDocument object in Level 2 Version 4 format: SBMLDocument* sbmlDoc = new SBMLDocument(2, 4); // Create a Model object inside the SBMLDocument object and set its // identifier: Model* model = sbmlDoc->createModel(); model->setId("EnzymaticReaction"); @endverbatim Next, we need to define some units in order to properly characterize the quantities used in the model. The libSBML API for unit definitions mirrors the SBML object model, which means there is an object class for UnitDefinition and an object class for the Unit objects that can be placed inside a UnitDefinition object to compose whatever unit is desired. The resulting code for creating units is a little bit tedious because there is a lot to set up, but on the other hand, it is quite straightforward and hopefully easy to follow: @verbatim // Temporary pointers (reused more than once below). UnitDefinition* unitdef; Unit* unit; // (UnitDefinition 1) Create a UnitDefinition object for "per_second". unitdef = model->createUnitDefinition(); unitdef->setId("per_second"); // Create an Unit inside the UnitDefinition object above: unit = unitdef->createUnit(); unit->setKind(UNIT_KIND_SECOND); unit->setExponent(-1); // (UnitDefinition 2) Create a UnitDefinition object for // "litre_per_mole_per_second". Note that we can reuse the pointers // 'unitdef' and 'unit' because the actual UnitDefinition object // (along with the Unit objects within it) is already attached to // the Model object. (Refer to the first two lines above.) unitdef = model->createUnitDefinition(); unitdef->setId("litre_per_mole_per_second"); // Create individual unit objects that will be put inside // the UnitDefinition to compose "litre_per_mole_per_second". unit = unitdef->createUnit(); unit->setKind(UNIT_KIND_MOLE); unit->setExponent(-1); unit = unitdef->createUnit(); unit->setKind(UNIT_KIND_LITRE); unit->setExponent(1); unit = unitdef->createUnit(); unit->setKind(UNIT_KIND_SECOND); unit->setExponent(-1); @endverbatim Before going on, we digress to explain one libSBML API idiom illustrated in the code above and in the rest of the program below. LibSBML provides two main mechanisms for creating objects: class constructors (e.g., @link Unit::Unit() Unit() @endlink), and createObject() methods provided by certain object classes such as UnitDefinition. The multiple mechanisms are provided by libSBML for flexibility and to support different use-cases, but they also have different implications for the overall model structure. In general, the recommended approach is to use the createObject() methods, as in the calls to UnitDefinition::createUnit() above. These methods create an instance and immediately link it into the parent, as well as return a pointer to the object created. By contrast, the other main way of creating an object involves first using the class constructor (e.g., @link Unit::Unit() Unit() @endlink) to create an object instance, and then calling an addObject(...) method to add the instance to another object. However, the addObject(...) methods in libSBML @em copy the object instance given to them. The result is that the caller has to be careful about making changes to the object: changes to the original will not affect the copy added to the model by the addObject(...) method, and moreover, the original must be explicitly deleted at some point or else the program will leak memory. In summary, unless special circumstances apply, most applications are best served using the various createObject() methods. Back to the example now, we create a compartment object: @verbatim // Create a string for the identifier of the compartment. This // will be used again later, so it's worth defining a constant here. const string compName = "cytosol"; // Create a Compartment object ("cytosol") Compartment* comp = model->createCompartment(); comp->setId(compName); // Set the "size" attribute of the Compartment object. We are not // setting the units on the compartment size explicitly, so the units // of this Compartment object will be the default SBML units of volume, // which are liters. comp->setSize(1e-14); @endverbatim Next, we turn to the species definitions in the model. For each species, the program needs to set an identifier, the compartment in which it's located, and optionally the initial quantity of the species in the compartment. @verbatim // Temporary pointer (reused more than once below). Species *sp; // Create the Species objects inside the Model object. // (Species1 ) Create a Species object for species "ES", and set its // compartment attribute to be the "cytosol" compartment in the model. sp = model->createSpecies(); sp->setCompartment(compName); sp->setId("ES"); // Set the "initialAmount" attribute of the Species object. // // In SBML, the units of a Species object's initial quantity are // determined by two attributes, "substanceUnits" and // "hasOnlySubstanceUnits", and the "spatialDimensions" attribute // of the Compartment object ("cytosol") in which the species // object is located. Here, we are using the default values for // "substanceUnits" (which is "mole") and "hasOnlySubstanceUnits" // (which is "false"). The compartment in which the species is // located uses volume units of liters, so the units of these // species (when the species appear in numerical formulas in the // model) will be moles/liter. // sp->setInitialAmount(0); // (Species 2) Create a Species object for "P" sp = model->createSpecies(); sp->setCompartment(compName); sp->setId("P"); sp->setInitialAmount(0); // (Species3) Create a Species object for "S" sp = model->createSpecies(); sp->setCompartment(compName); sp->setId("S"); sp->setInitialAmount(1e-20); // (Species4) Create a Species object for "E" sp = model->createSpecies(); sp->setCompartment(compName); sp->setId("E"); sp->setInitialAmount(5e-21 ); @endverbatim It remains only to create objects representing the reactions in the model. In this example, we will create the system as one reversible and one irreversible reaction; alternatively, the system could have been created as three irreversible reactions. This program has hard-coded versions of the reaction rate formulas only for the purposes of this example. Of course, in a more practical real-life situation, the formulas would not be fixed in this way—one would have methods to read/write/translate expressions originating from different sources (users, files, other programs etc.) to/from the Abstract Syntax Trees (ASTs) used by libSBML. @verbatim // Temporary pointers reused below. Reaction* reaction; SpeciesReference* spr; Parameter* para; KineticLaw* kl; // Create Reactant objects inside the Reaction object ("veq"). reaction = model->createReaction(); reaction->setId("veq"); // (Reactant 1) Create a Reactant object that references Species "E" // in the model. The object will be created within the reaction in the // SBML . spr = reaction->createReactant(); spr->setSpecies("E"); // (Reactant 2) Create a Reactant object that references Species "S" // in the model. spr = reaction->createReactant(); spr->setSpecies("S"); // Create a Product object that references Species "ES" in the model. spr = reaction->createProduct(); spr->setSpecies("ES"); // Create a KineticLaw object inside the Reaction object ("veq"). kl = reaction->createKineticLaw(); // Next, we will create an Abstract Syntax Tree representation of the // following MathML expression: // // // // // cytosol // // // // // kon // E // S // // // // koff // ES // // // // // First, create nodes representing the variables: ASTNode* astCytosol = new ASTNode(AST_NAME); astCytosol->setName("cytosol"); ASTNode* astKon = new ASTNode(AST_NAME); astKon->setName("kon"); ASTNode* astKoff = new ASTNode(AST_NAME); astKoff->setName("koff"); ASTNode* astE = new ASTNode(AST_NAME); astE->setName("E"); ASTNode* astS = new ASTNode(AST_NAME); astS->setName("S"); ASTNode* astES = new ASTNode(AST_NAME); astES->setName("ES"); // Next, we build up the tree bottom-up, beginning with the // following part: // // // // koff // ES // ASTNode *astTimes1 = new ASTNode(AST_TIMES); astTimes1->addChild(astKoff); astTimes1->addChild(astES); // Continuing, we now create nodes representing the following part: // // // // kon // E // S // // // Note: it is worth pointing out that the expression above has 3 // elements inside the , but libSBML's ASTs are binary trees, // and only accept two operands. The following code illustrates how // to construct such an expression using binary ASTNode objects. ASTNode *astTimes2 = new ASTNode(AST_TIMES); astTimes2->addChild(astE); astTimes2->addChild(astS); ASTNode *astTimes = new ASTNode(AST_TIMES); astTimes->addChild(astKon); astTimes->addChild(astTimes2); // Now for a node representing the following part: // // // // // // kon // E // S // // // // koff // ES // // ASTNode *astMinus = new ASTNode(AST_MINUS); astMinus->addChild(astTimes); astMinus->addChild(astTimes1); // And finally, the top node, representing the complete expression: // // // // cytosol // // // // // kon // E // S // // // // koff // ES // // // ASTNode* astMath = new ASTNode(AST_TIMES); astMath->addChild(astCytosol); astMath->addChild(astMinus); // Now set the element in the KineticLaw object: kl->setMath(astMath); // KineticLaw::setMath(const ASTNode*) sets the math of the // KineticLaw object to a COPY of the given ASTNode. The caller // (i.e., this program) should delete the original ASTNode object to // avoid a memory leak. delete astMath; // Create local Parameter objects for "kon" and "koff" inside the // KineticLaw object. Note how we set the units on the parameter // values. para = kl->createParameter(); para->setId("kon"); para->setValue(1000000); para->setUnits("litre_per_mole_per_second"); para = kl->createParameter(); para->setId("koff"); para->setValue(0.2); para->setUnits("per_second"); @endverbatim That was one reaction. One more reaction to go. This time around, for variety and to illustrate another one of libSBML's facilities, we create the AST for the reaction rate formula starting from a string in MathML form. @verbatim // (Reaction 2) Create a Reaction object ("vcat") . reaction = model->createReaction(); reaction->setId("vcat"); reaction->setReversible(false); // Create Reactant objects inside the Reaction object ("vcat"). // (Reactant 1) Create a Reactant object that references Species // "ES" in the model. spr = reaction->createReactant(); spr->setSpecies("ES"); // (Product 1) Create a Product object that references Species "E" // in the model. spr = reaction->createProduct(); spr->setSpecies("E"); // (Product 2) Create a Product object that references Species "P" // in the model. spr = reaction->createProduct(); spr->setSpecies("P"); // Create a KineticLaw object inside the Reaction object ("vcat"). kl = reaction->createKineticLaw(); // To create mathematical expressions, one would typically construct // an ASTNode tree as the above example code which creates a math of // another KineticLaw object. Here, to save some space and // illustrate another approach of doing it, we will write out the // formula in MathML form and then use a libSBML convenience // function to create the ASTNode tree for us. (However, this is a // bit dangerous; it's very easy to make mistakes when writing // MathML by hand, so in a real program, we would not really want to // do it this way. We do this here mostly to illustrate this approach.) string mathXMLString = "" " " " " " cytosol " " kcat " " ES " " " ""; // The following reuses 'astMath' declared earlier in this program. // Note also how this is careful to delete the structure after handing // it to setMath(). As noted earlier, setMath() copies the structure // given to it, so if the caller is finished using it, the structure // should be deleted in order to avoid leaking memory. astMath = readMathMLFromString(mathXMLString.c_str()); kl->setMath(astMath); delete astMath; // Create a local Parameter object for "kcat" inside the KineticLaw // object. para = kl->createParameter(); para->setId("kcat"); para->setValue(0.1); para->setUnits("per_second"); @endverbatim That concludes all the model creation code. Our createExampleEnzymaticReaction() function now only needs to return the SBMLDocument object that it created. @verbatim // Return the created SBMLDocument object. // The returned object must be explicitly deleted by the caller, // otherwise a memory leak will happen. return sbmlDoc; } @endverbatim @subsection validating Function for verifying the model One of the most important features of libSBML is its ability to perform SBML validation to ensure that a model adheres to the SBML specification for whatever Level+Version combination the model uses. Although the example program presented here is short and we can verify its correctness by inspection, in more general programs, it is difficult to ensure that everything in a model is constructed properly. The ability to perform automatic checks then becomes very useful. We therefore take this opportunity to illustrate the use of libSBML's validation and checking facilities. The rationale behind libSBML's consistency-checking and validation functions is the following. In libSBML versions up to and including the version 3.3.x series, the individual methods for creating and setting attributes and other components are quite lenient, and allow a caller to compose SBML entities that may not, in the end, represent valid SBML. This allows applications the freedom to do things such as save incomplete models (which is useful when models are being developed over long periods of time), but at the same time, it means that a separate validation step is necessary when a calling program finally wants to finish a complete SBML document. LibSBML implements verification of SBML in two steps, represented by the two methods SBMLDocument::checkInternalConsistency() and SBMLDocument::checkConsistency(). The former verifies the basic internal consistency and syntax of an SBML document, and the latter implements more elaborate validation rules (both those defined by the SBML specifications, as well as additional rules offered by libSBML). When an application builds up a model programmatically and is finally ready to say "yes, this model is finished", it should call both of these methods to help ensure the correctness and consistency of the finished result. Conversely, if it is reading a model from a file or data stream, libSBML will automatically verify the syntax and schema consistency of the document. The application only needs to check the results by interrogating the error log on the SBMLDocument object, then call SBMLDocument::checkConsistency() for checking the model against the SBML validation rules. The following code illustrates a complete sequence of checking a newly-created model using both SBMLDocument::checkInternalConsistency() and SBMLDocument::checkConsistency(), and handling the reporting of problems to the user. @verbatim bool validateExampleSBML (SBMLDocument* sbmlDoc) { if (!sbmlDoc) { cerr << "validateExampleSBML: given a null SBML Document" << endl; return false; } string consistencyMessages; string validationMessages; bool noProblems = true; unsigned int numCheckFailures = 0; unsigned int numConsistencyErrors = 0; unsigned int numConsistencyWarnings = 0; unsigned int numValidationErrors = 0; unsigned int numValidationWarnings = 0; // LibSBML 3.3 is lenient when generating models from scratch using the // API for creating objects. Once the whole model is done and before it // gets written out, it's important to check that the whole model is in // fact complete, consistent and valid. numCheckFailures = sbmlDoc->checkInternalConsistency(); if ( numCheckFailures > 0 ) { noProblems = false; for (unsigned int i = 0; i < numCheckFailures; i++) { const SBMLError* sbmlErr = sbmlDoc->getError(i); if ( sbmlErr->isFatal() || sbmlErr->isError() ) { ++numConsistencyErrors; } else { ++numConsistencyWarnings; } } ostringstream oss; sbmlDoc->printErrors(oss); consistencyMessages = oss.str(); } // If the internal checks fail, it makes little sense to attempt // further validation, because the model may be too compromised to // be properly interpreted. if (numConsistencyErrors > 0) { consistencyMessages += "Further validation aborted."; } else { numCheckFailures = sbmlDoc->checkConsistency(); if ( numCheckFailures > 0 ) { noProblems = false; for (unsigned int i = 0; i < numCheckFailures; i++) { const SBMLError* sbmlErr = sbmlDoc->getError(i); if ( sbmlErr->isFatal() || sbmlErr->isError() ) { ++numValidationErrors; } else { ++numValidationWarnings; } } ostringstream oss; sbmlDoc->printErrors(oss); validationMessages = oss.str(); } } if (noProblems) return true; else { if (numConsistencyErrors > 0) { cout << "ERROR: encountered " << numConsistencyErrors << " consistency error" << (numConsistencyErrors == 1 ? "" : "s") << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl; } if (numConsistencyWarnings > 0) { cout << "Notice: encountered " << numConsistencyWarnings << " consistency warning" << (numConsistencyWarnings == 1 ? "" : "s") << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl; } cout << endl << consistencyMessages; if (numValidationErrors > 0) { cout << "ERROR: encountered " << numValidationErrors << " validation error" << (numValidationErrors == 1 ? "" : "s") << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl; } if (numValidationWarnings > 0) { cout << "Notice: encountered " << numValidationWarnings << " validation warning" << (numValidationWarnings == 1 ? "" : "s") << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl; } cout << endl << validationMessages; return (numConsistencyErrors == 0 && numValidationErrors == 0); } } @endverbatim @subsection writing Function for writing the model The final step of writing out a model is relatively trivial. LibSBML provides a class of object (SBMLWriter) that provides the necessary functionality. All that is left to do is to call it with a desired file name. @verbatim /** * * Writes the given SBMLDocument to the given file. * */ bool writeExampleSBML(const SBMLDocument* sbmlDoc, const string& filename) { SBMLWriter sbmlWriter; bool result = sbmlWriter.writeSBML(sbmlDoc, filename); if (result) { cout << "Wrote file \"" << filename << "\"" << endl; return true; } else { cerr << "Failed to write \"" << filename << "\"" << endl; return false; } } @endverbatim @section source Source code An expanded version of this example is supplied with the libSBML distribution as the file named createExampleSBML.cpp in the directory examples/c++. */