/** @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:
Entity in the %Model |
Initial value |
Units |
cytosol | \f$1 \times 10^{-14}\f$ | litres |
S | \f$1 \times 10^{-20}\f$ | moles |
P | 0 | moles |
E | \f$5 \times 10^{-21}\f$ | moles |
ES | 0 | | moles
kon | 1 000 000 | litre/mole/second |
koff | 0.2 | 1/second |
kcat | | 1/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:
//
//
// 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