1 /**
2 * @file createExampleSBML.cpp
3 * @brief Creates example SBML models presented in the SBML specification.
4 * @author Akiya Jouraku
5 * @author Michael Hucka
6 * @author Sarah Keating
7 *
8 * <!--------------------------------------------------------------------------
9 * This sample program is distributed under a different license than the rest
10 * of libSBML. This program uses the open-source MIT license, as follows:
11 *
12 * Copyright (c) 2013-2018 by the California Institute of Technology
13 * (California, USA), the European Bioinformatics Institute (EMBL-EBI, UK)
14 * and the University of Heidelberg (Germany), with support from the National
15 * Institutes of Health (USA) under grant R01GM070923. All rights reserved.
16 *
17 * Permission is hereby granted, free of charge, to any person obtaining a
18 * copy of this software and associated documentation files (the "Software"),
19 * to deal in the Software without restriction, including without limitation
20 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
21 * and/or sell copies of the Software, and to permit persons to whom the
22 * Software is furnished to do so, subject to the following conditions:
23 *
24 * The above copyright notice and this permission notice shall be included in
25 * all copies or substantial portions of the Software.
26 *
27 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
28 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
29 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
30 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
31 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
32 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
33 * DEALINGS IN THE SOFTWARE.
34 *
35 * Neither the name of the California Institute of Technology (Caltech), nor
36 * of the European Bioinformatics Institute (EMBL-EBI), nor of the University
37 * of Heidelberg, nor the names of any contributors, may be used to endorse
38 * or promote products derived from this software without specific prior
39 * written permission.
40 * ------------------------------------------------------------------------ -->
41 */
42
43 #include <iostream>
44 #include <sbml/SBMLTypes.h>
45
46 using namespace std;
47 LIBSBML_CPP_NAMESPACE_USE
48
49 //
50 // Functions for creating the Example SBML documents.
51 //
52 SBMLDocument* createExampleEnzymaticReaction(); /* 7.1 */
53 SBMLDocument* createExampleInvolvingUnits(); /* 7.2 */
54 SBMLDocument* createExampleInvolvingFunctionDefinitions(); /* 7.8 */
55
56 //
57 // Helper functions for validating and writing the SBML documents created.
58 //
59 bool validateExampleSBML(SBMLDocument *sbmlDoc);
60 bool writeExampleSBML(const SBMLDocument *sbmlDoc, const string& filename);
61
62 //
63 // These variables are used in writeExampleSBML when writing an SBML
64 // document. They are handed to libSBML functions in order to include
65 // the program information into comments within the SBML file.
66 //
67 const static string ProgramName = "createExampleModels";
68 const static string ProgramVersion = "1.0.0";
69
70 //
71 // The SBML Level and Version of the example SBML models.
72 //
73 const static unsigned int Level = 2;
74 const static unsigned int Version = 4;
75
76
77 //===============================================================================
78 //
79 // Main routine
80 //
81 // Creates SBML models represented in "Example models expressed in XML using
82 // SBML" in Section 7 of the SBML Level 2 Version 4 specification(*).
83 //
84 // (*) The specification document is available at the following URL:
85 // http://sbml.org/Documents/Specifications
86 //
87 //===============================================================================
88 //
89 int
main(int argc,char * argv[])90 main (int argc, char *argv[])
91 {
92 SBMLDocument* sbmlDoc = 0;
93 bool SBMLok = false;
94
95 try
96 {
97 //-------------------------------------------------
98 // 7.1 A Simple example application of SBML
99 //-------------------------------------------------
100
101 sbmlDoc = createExampleEnzymaticReaction();
102 SBMLok = validateExampleSBML(sbmlDoc);
103 if (SBMLok) writeExampleSBML(sbmlDoc, "enzymaticreaction.xml");
104 delete sbmlDoc;
105 if (!SBMLok) return 1;
106
107 //-------------------------------------------------
108 // 7.2 Example involving units
109 //-------------------------------------------------
110
111 sbmlDoc = createExampleInvolvingUnits();
112 SBMLok = validateExampleSBML(sbmlDoc);
113 if (SBMLok) writeExampleSBML(sbmlDoc, "units.xml");
114 delete sbmlDoc;
115 if (!SBMLok) return 1;
116
117 //-------------------------------------------------
118 // 7.8 Example involving function definitions
119 //-------------------------------------------------
120
121 sbmlDoc = createExampleInvolvingFunctionDefinitions();
122 SBMLok = validateExampleSBML(sbmlDoc);
123 if (SBMLok) writeExampleSBML(sbmlDoc, "functiondef.xml");
124 delete sbmlDoc;
125 if (!SBMLok) return 1;
126
127 }
128 catch (std::bad_alloc& e)
129 {
130 cerr << e.what() << ": Unable to allocate memory." << endl;
131 return 1;
132 }
133 catch (...)
134 {
135 cerr << "Unexpected exceptional condition encountered." << endl;
136 return 1;
137 }
138
139 // A 0 return status is the standard Unix/Linux way to say "all ok".
140 return 0;
141
142 }
143
144
145 //===============================================================================
146 //
147 //
148 // Functions for creating the Example SBML documents.
149 //
150 //
151 //===============================================================================
152
153
154 /**
155 *
156 * Creates an SBML model represented in "7.1 A Simple example application of SBML"
157 * in the SBML Level 2 Version 4 Specification.
158 *
159 */
createExampleEnzymaticReaction()160 SBMLDocument* createExampleEnzymaticReaction()
161 {
162 const unsigned int level = Level;
163 const unsigned int version = Version;
164
165 //---------------------------------------------------------------------------
166 //
167 // Creates an SBMLDocument object
168 //
169 //---------------------------------------------------------------------------
170
171 SBMLDocument* sbmlDoc = new SBMLDocument(level,version);
172
173 //---------------------------------------------------------------------------
174 //
175 // Creates a Model object inside the SBMLDocument object.
176 //
177 //---------------------------------------------------------------------------
178
179 Model* model = sbmlDoc->createModel();
180 model->setId("EnzymaticReaction");
181
182 //---------------------------------------------------------------------------
183 //
184 // Creates UnitDefinition objects inside the Model object.
185 //
186 //---------------------------------------------------------------------------
187
188 // Temporary pointers (reused more than once below).
189
190 UnitDefinition* unitdef;
191 Unit* unit;
192
193 //---------------------------------------------------------------------------
194 // (UnitDefinition1) Creates an UnitDefinition object ("per_second")
195 //---------------------------------------------------------------------------
196
197 unitdef = model->createUnitDefinition();
198 unitdef->setId("per_second");
199
200 // Creates an Unit inside the UnitDefinition object
201
202 unit = unitdef->createUnit();
203 unit->setKind(UNIT_KIND_SECOND);
204 unit->setExponent(-1);
205
206 //--------------------------------------------------------------------------------
207 // (UnitDefinition2) Creates an UnitDefinition object ("litre_per_mole_per_second")
208 //--------------------------------------------------------------------------------
209
210 // Note that we can reuse the pointers 'unitdef' and 'unit' because the
211 // actual UnitDefinition object (along with the Unit objects within it)
212 // is already attached to the Model object.
213
214 unitdef = model->createUnitDefinition();
215 unitdef->setId("litre_per_mole_per_second");
216
217 // Creates an Unit inside the UnitDefinition object ("litre_per_mole_per_second")
218
219 unit = unitdef->createUnit();
220 unit->setKind(UNIT_KIND_MOLE);
221 unit->setExponent(-1);
222
223 // Creates an Unit inside the UnitDefinition object ("litre_per_mole_per_second")
224
225 unit = unitdef->createUnit();
226 unit->setKind(UNIT_KIND_LITRE);
227 unit->setExponent(1);
228
229 // Creates an Unit inside the UnitDefinition object ("litre_per_mole_per_second")
230
231 unit = unitdef->createUnit();
232 unit->setKind(UNIT_KIND_SECOND);
233 unit->setExponent(-1);
234
235
236 //---------------------------------------------------------------------------
237 //
238 // Creates a Compartment object inside the Model object.
239 //
240 //---------------------------------------------------------------------------
241
242 Compartment* comp;
243 const string compName = "cytosol";
244
245 // Creates a Compartment object ("cytosol")
246
247 comp = model->createCompartment();
248 comp->setId(compName);
249
250 // Sets the "size" attribute of the Compartment object.
251 //
252 // We are not setting the units on the compartment size explicitly, so
253 // the units of this Compartment object will be the default SBML units of
254 // volume, which are liters.
255 //
256 comp->setSize(1e-14);
257
258
259 //---------------------------------------------------------------------------
260 //
261 // Creates Species objects inside the Model object.
262 //
263 //---------------------------------------------------------------------------
264
265 // Temporary pointer (reused more than once below).
266
267 Species *sp;
268
269 //---------------------------------------------------------------------------
270 // (Species1) Creates a Species object ("ES")
271 //---------------------------------------------------------------------------
272
273 // Create the Species objects inside the Model object.
274
275 sp = model->createSpecies();
276 sp->setId("ES");
277 sp->setName("ES");
278
279 // Sets the "compartment" attribute of the Species object to identify the
280 // compartment in which the Species object is located.
281
282 sp->setCompartment(compName);
283
284 // Sets the "initialAmount" attribute of the Species object.
285 //
286 // In SBML, the units of a Species object's initial quantity are
287 // determined by two attributes, "substanceUnits" and
288 // "hasOnlySubstanceUnits", and the "spatialDimensions" attribute
289 // of the Compartment object ("cytosol") in which the species
290 // object is located. Here, we are using the default values for
291 // "substanceUnits" (which is "mole") and "hasOnlySubstanceUnits"
292 // (which is "false"). The compartment in which the species is
293 // located uses volume units of liters, so the units of these
294 // species (when the species appear in numerical formulas in the
295 // model) will be moles/liters.
296 //
297 sp->setInitialAmount(0);
298
299 //---------------------------------------------------------------------------
300 // (Species2) Creates a Species object ("P")
301 //---------------------------------------------------------------------------
302
303 sp = model->createSpecies();
304 sp->setCompartment(compName);
305 sp->setId("P");
306 sp->setName("P");
307 sp->setInitialAmount(0);
308
309 //---------------------------------------------------------------------------
310 // (Species3) Creates a Species object ("S")
311 //---------------------------------------------------------------------------
312
313 sp = model->createSpecies();
314 sp->setCompartment(compName);
315 sp->setId("S");
316 sp->setName("S");
317 sp->setInitialAmount(1e-20);
318
319 //---------------------------------------------------------------------------
320 // (Species4) Creates a Species object ("E")
321 //---------------------------------------------------------------------------
322
323 sp = model->createSpecies();
324 sp->setCompartment(compName);
325 sp->setId("E");
326 sp->setName("E");
327 sp->setInitialAmount(5e-21);
328
329
330 //---------------------------------------------------------------------------
331 //
332 // Creates Reaction objects inside the Model object.
333 //
334 //---------------------------------------------------------------------------
335
336 // Temporary pointers.
337
338 Reaction* reaction;
339 SpeciesReference* spr;
340 KineticLaw* kl;
341
342 //---------------------------------------------------------------------------
343 // (Reaction1) Creates a Reaction object ("veq").
344 //---------------------------------------------------------------------------
345
346 reaction = model->createReaction();
347 reaction->setId("veq");
348
349 // (Reactant1) Creates a Reactant object that references Species "E"
350 // in the model. The object will be created within the reaction in the
351 // SBML <listOfReactants>.
352
353 spr = reaction->createReactant();
354 spr->setSpecies("E");
355
356 // (Reactant2) Creates a Reactant object that references Species "S"
357 // in the model.
358
359 spr = reaction->createReactant();
360 spr->setSpecies("S");
361
362 //---------------------------------------------------------------------------
363 // (Product1) Creates a Product object that references Species "ES" in
364 // the model.
365 //---------------------------------------------------------------------------
366
367 spr = reaction->createProduct();
368 spr->setSpecies("ES");
369
370 //---------------------------------------------------------------------------
371 // Creates a KineticLaw object inside the Reaction object ("veq").
372 //---------------------------------------------------------------------------
373
374 kl = reaction->createKineticLaw();
375
376 //---------------------------------------------------------------------------
377 // Creates an ASTNode object which represents the following math of the
378 // KineticLaw.
379 //
380 // <math xmlns="http://www.w3.org/1998/Math/MathML">
381 // <apply>
382 // <times/>
383 // <ci> cytosol </ci>
384 // <apply>
385 // <minus/>
386 // <apply>
387 // <times/>
388 // <ci> kon </ci>
389 // <ci> E </ci>
390 // <ci> S </ci>
391 // </apply>
392 // <apply>
393 // <times/>
394 // <ci> koff </ci>
395 // <ci> ES </ci>
396 // </apply>
397 // </apply>
398 // </apply>
399 // </math>
400 //
401 //---------------------------------------------------------------------------
402
403 //------------------------------------------
404 //
405 // create nodes representing the variables
406 //
407 //------------------------------------------
408
409 ASTNode* astCytosol = new ASTNode(AST_NAME);
410 astCytosol->setName("cytosol");
411
412 ASTNode* astKon = new ASTNode(AST_NAME);
413 astKon->setName("kon");
414
415 ASTNode* astKoff = new ASTNode(AST_NAME);
416 astKoff->setName("koff");
417
418 ASTNode* astE = new ASTNode(AST_NAME);
419 astE->setName("E");
420
421 ASTNode* astS = new ASTNode(AST_NAME);
422 astS->setName("S");
423
424 ASTNode* astES = new ASTNode(AST_NAME);
425 astES->setName("ES");
426
427
428 //--------------------------------------------
429 //
430 // create node representing
431 // <apply>
432 // <times/>
433 // <ci> koff </ci>
434 // <ci> ES </ci>
435 // </apply>
436 //
437 //--------------------------------------------
438
439 ASTNode *astTimes1 = new ASTNode(AST_TIMES);
440 astTimes1->addChild(astKoff);
441 astTimes1->addChild(astES);
442
443 //--------------------------------------------
444 //
445 // create node representing
446 // <apply>
447 // <times/>
448 // <ci> kon </ci>
449 // <ci> E </ci>
450 // <ci> S </ci>
451 // </apply>
452 //
453 //
454 // (NOTES)
455 //
456 // Since there is a restriction with an ASTNode of "<times/>" operation
457 // such that the ASTNode is a binary class and thus only two operands can
458 // be directly added, the following code in this comment block is invalid
459 // because the code directly adds three <ci> ASTNodes to <times/> ASTNode.
460 //
461 // ASTNode *astTimes = new ASTNode(AST_TIMES);
462 // astTimes->addChild(astKon);
463 // astTimes->addChild(astE);
464 // astTimes->addChild(astS);
465 //
466 // The following valid code after this comment block creates the ASTNode
467 // as a binary tree.
468 //
469 // Please see "Converting between ASTs and text strings" described
470 // at http://sbml.org/Software/libSBML/docs/cpp-api/class_a_s_t_node.html
471 // for the detailed information.
472 //
473 //--------------------------------------------
474
475 ASTNode *astTimes2 = new ASTNode(AST_TIMES);
476 astTimes2->addChild(astE);
477 astTimes2->addChild(astS);
478
479 ASTNode *astTimes = new ASTNode(AST_TIMES);
480 astTimes->addChild(astKon);
481 astTimes->addChild(astTimes2);
482
483 //--------------------------------------------
484 //
485 // create node representing
486 // <apply>
487 // <minus/>
488 // <apply>
489 // <times/>
490 // <ci> kon </ci>
491 // <ci> E </ci>
492 // <ci> S </ci>
493 // </apply>
494 // <apply>
495 // <times/>
496 // <ci> koff </ci>
497 // <ci> ES </ci>
498 // </apply>
499 // </apply>
500 //
501 //--------------------------------------------
502
503 ASTNode *astMinus = new ASTNode(AST_MINUS);
504 astMinus->addChild(astTimes);
505 astMinus->addChild(astTimes1);
506
507
508 //--------------------------------------------
509 //
510 // create node representing
511 // <apply>
512 // <times/>
513 // <ci> cytosol </ci>
514 // <apply>
515 // <minus/>
516 // <apply>
517 // <times/>
518 // <ci> kon </ci>
519 // <ci> E </ci>
520 // <ci> S </ci>
521 // </apply>
522 // <apply>
523 // <times/>
524 // <ci> koff </ci>
525 // <ci> ES </ci>
526 // </apply>
527 // </apply>
528 // </apply>
529 //
530 //--------------------------------------------
531
532 ASTNode* astMath = new ASTNode(AST_TIMES);
533 astMath->addChild(astCytosol);
534 astMath->addChild(astMinus);
535
536 //---------------------------------------------
537 //
538 // set the Math element
539 //
540 //------------------------------------------------
541
542 kl->setMath(astMath);
543
544 // KineticLaw::setMath(const ASTNode*) sets the math of the KineticLaw object
545 // to a copy of the given ASTNode, and thus basically the caller should delete
546 // the original ASTNode object if the caller has the ownership of the object to
547 // avoid memory leak.
548
549 delete astMath;
550
551
552 //---------------------------------------------------------------------------
553 // Creates local Parameter objects inside the KineticLaw object.
554 //---------------------------------------------------------------------------
555
556 // Creates a Parameter ("kon")
557
558 Parameter* para = kl->createParameter();
559 para->setId("kon");
560 para->setValue(1000000);
561 para->setUnits("litre_per_mole_per_second");
562
563 // Creates a Parameter ("koff")
564
565 para = kl->createParameter();
566 para->setId("koff");
567 para->setValue(0.2);
568 para->setUnits("per_second");
569
570
571 //---------------------------------------------------------------------------
572 // (Reaction2) Creates a Reaction object ("vcat") .
573 //---------------------------------------------------------------------------
574
575 reaction = model->createReaction();
576 reaction->setId("vcat");
577 reaction->setReversible(false);
578
579 //---------------------------------------------------------------------------
580 // Creates Reactant objects inside the Reaction object ("vcat").
581 //---------------------------------------------------------------------------
582
583 // (Reactant1) Creates a Reactant object that references Species "ES" in the
584 // model.
585
586 spr = reaction->createReactant();
587 spr->setSpecies("ES");
588
589 //---------------------------------------------------------------------------
590 // Creates a Product object inside the Reaction object ("vcat").
591 //---------------------------------------------------------------------------
592
593 // (Product1) Creates a Product object that references Species "E" in the model.
594
595 spr = reaction->createProduct();
596 spr->setSpecies("E");
597
598 // (Product2) Creates a Product object that references Species "P" in the model.
599
600 spr = reaction->createProduct();
601 spr->setSpecies("P");
602
603 //---------------------------------------------------------------------------
604 // Creates a KineticLaw object inside the Reaction object ("vcat").
605 //---------------------------------------------------------------------------
606
607 kl = reaction->createKineticLaw();
608
609 //---------------------------------------------------------------------------
610 // Sets a math (ASTNode object) to the KineticLaw object.
611 //---------------------------------------------------------------------------
612
613 // To create mathematical expressions, one would typically construct
614 // an ASTNode tree as the above example code which creates a math of another
615 // KineticLaw object. Here, to save some space and illustrate another approach
616 // of doing it, we will write out the formula in MathML form and then use a
617 // libSBML convenience function to create the ASTNode tree for us.
618 // (This is a bit dangerous; it's very easy to make mistakes when writing MathML
619 // by hand, so in a real program, we would not really want to do it this way.)
620
621 string mathXMLString = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">"
622 " <apply>"
623 " <times/>"
624 " <ci> cytosol </ci>"
625 " <ci> kcat </ci>"
626 " <ci> ES </ci>"
627 " </apply>"
628 "</math>";
629
630 astMath = readMathMLFromString(mathXMLString.c_str());
631 kl->setMath(astMath);
632 delete astMath;
633
634 //---------------------------------------------------------------------------
635 // Creates local Parameter objects inside the KineticLaw object.
636 //---------------------------------------------------------------------------
637
638 // Creates a Parameter ("kcat")
639
640 para = kl->createParameter();
641 para->setId("kcat");
642 para->setValue(0.1);
643 para->setUnits("per_second");
644
645
646 // Returns the created SBMLDocument object.
647 // The returned object must be explicitly deleted by the caller,
648 // otherwise a memory leak will happen.
649
650 return sbmlDoc;
651
652 }
653
654
655 /**
656 *
657 * Creates an SBML model represented in "7.2 Example involving units"
658 * in the SBML Level 2 Version 4 Specification.
659 *
660 */
createExampleInvolvingUnits()661 SBMLDocument* createExampleInvolvingUnits()
662 {
663 const unsigned int level = Level;
664 const unsigned int version = Version;
665
666 //---------------------------------------------------------------------------
667 //
668 // Creates an SBMLDocument object
669 //
670 //---------------------------------------------------------------------------
671
672 SBMLDocument* sbmlDoc = new SBMLDocument(level,version);
673
674 // Adds the namespace for XHTML to the SBMLDocument object. We need this
675 // because we will add notes to the model. (By default, the SBML document
676 // created by SBMLDocument only declares the SBML XML namespace.)
677
678 sbmlDoc->getNamespaces()->add("http://www.w3.org/1999/xhtml", "xhtml");
679
680 //---------------------------------------------------------------------------
681 //
682 // Creates a Model object inside the SBMLDocument object.
683 //
684 //---------------------------------------------------------------------------
685
686 Model* model = sbmlDoc->createModel();
687 model->setId("unitsExample");
688
689 //---------------------------------------------------------------------------
690 //
691 // Creates UnitDefinition objects inside the Model object.
692 //
693 //---------------------------------------------------------------------------
694
695 // Temporary pointers (reused more than once below).
696
697 UnitDefinition* unitdef;
698 Unit *unit;
699
700 //---------------------------------------------------------------------------
701 // (UnitDefinition1) Creates an UnitDefinition object ("substance").
702 //
703 // This has the effect of redefining the default unit of subtance for the
704 // whole model.
705 //---------------------------------------------------------------------------
706
707 unitdef = model->createUnitDefinition();
708 unitdef->setId("substance");
709
710 // Creates an Unit inside the UnitDefinition object
711
712 unit = unitdef->createUnit();
713 unit->setKind(UNIT_KIND_MOLE);
714 unit->setScale(-3);
715
716 //--------------------------------------------------------------------------------
717 // (UnitDefinition2) Creates an UnitDefinition object ("mmls")
718 //--------------------------------------------------------------------------------
719
720 // Note that we can reuse the pointers 'unitdef' and 'unit' because the
721 // actual UnitDefinition object (along with the Unit objects within it)
722 // is already attached to the Model object.
723
724 unitdef = model->createUnitDefinition();
725 unitdef->setId("mmls");
726
727 // Creates an Unit inside the UnitDefinition object ("mmls")
728
729 unit = unitdef->createUnit();
730 unit->setKind(UNIT_KIND_MOLE);
731 unit->setScale(-3);
732
733 // Creates an Unit inside the UnitDefinition object ("mmls")
734
735 unit = unitdef->createUnit();
736 unit->setKind(UNIT_KIND_LITRE);
737 unit->setExponent(-1);
738
739 // Creates an Unit inside the UnitDefinition object ("mmls")
740
741 unit = unitdef->createUnit();
742 unit->setKind(UNIT_KIND_SECOND);
743 unit->setExponent(-1);
744
745 //--------------------------------------------------------------------------------
746 // (UnitDefinition3) Creates an UnitDefinition object ("mml")
747 //--------------------------------------------------------------------------------
748
749 unitdef = model->createUnitDefinition();
750 unitdef->setId("mml");
751
752 // Creates an Unit inside the UnitDefinition object ("mml")
753
754 unit = unitdef->createUnit();
755 unit->setKind(UNIT_KIND_MOLE);
756 unit->setScale(-3);
757
758 // Creates an Unit inside the UnitDefinition object ("mml")
759
760 unit = unitdef->createUnit();
761 unit->setKind(UNIT_KIND_LITRE);
762 unit->setExponent(-1);
763
764
765 //---------------------------------------------------------------------------
766 //
767 // Creates a Compartment object inside the Model object.
768 //
769 //---------------------------------------------------------------------------
770
771 Compartment* comp;
772 const string compName = "cell";
773
774 // Creates a Compartment object ("cell")
775
776 comp = model->createCompartment();
777 comp->setId(compName);
778
779 // Sets the "size" attribute of the Compartment object.
780 //
781 // The units of this Compartment object is the default SBML
782 // units of volume (litre), and thus we don't have to explicitly invoke
783 // setUnits("litre") function to set the default units.
784 //
785 comp->setSize(1);
786
787
788 //---------------------------------------------------------------------------
789 //
790 // Creates Species objects inside the Model object.
791 //
792 //---------------------------------------------------------------------------
793
794 // Temporary pointer (reused more than once below).
795
796 Species *sp;
797
798 //---------------------------------------------------------------------------
799 // (Species1) Creates a Species object ("x0")
800 //---------------------------------------------------------------------------
801
802 sp = model->createSpecies();
803 sp->setId("x0");
804
805 // Sets the "compartment" attribute of the Species object to identify the
806 // compartnet in which the Species object located.
807
808 sp->setCompartment(compName);
809
810 // Sets the "initialConcentration" attribute of the Species object.
811 //
812 // The units of this Species object is determined by two attributes of this
813 // Species object ("substanceUnits" and "hasOnlySubstanceUnits") and the
814 // "spatialDimensions" attribute of the Compartment object ("cytosol") in which
815 // this species object is located.
816 // Since the default values are used for "substanceUnits" (substance (mole))
817 // and "hasOnlySubstanceUnits" (false) and the value of "spatialDimension" (3)
818 // is greater than 0, the units of this Species object is moles/liters .
819 //
820 sp->setInitialConcentration(1);
821
822 //---------------------------------------------------------------------------
823 // (Species2) Creates a Species object ("x1")
824 //---------------------------------------------------------------------------
825
826 sp = model->createSpecies();
827 sp->setId("x1");
828 sp->setCompartment(compName);
829 sp->setInitialConcentration(1);
830
831 //---------------------------------------------------------------------------
832 // (Species3) Creates a Species object ("s1")
833 //---------------------------------------------------------------------------
834
835 sp = model->createSpecies();
836 sp->setCompartment(compName);
837 sp->setId("s1");
838 sp->setInitialConcentration(1);
839
840 //---------------------------------------------------------------------------
841 // (Species4) Creates a Species object ("s2")
842 //---------------------------------------------------------------------------
843
844 sp = model->createSpecies();
845 sp->setCompartment(compName);
846 sp->setId("s2");
847 sp->setInitialConcentration(1);
848
849 //---------------------------------------------------------------------------
850 //
851 // Creates global Parameter objects inside the Model object.
852 //
853 //---------------------------------------------------------------------------
854
855 Parameter* para;
856
857 // Creates a Parameter ("vm")
858
859 para = model->createParameter();
860 para->setId("vm");
861 para->setValue(2);
862 para->setUnits("mmls");
863
864 // Creates a Parameter ("km")
865
866 para = model->createParameter();
867 para->setId("km");
868 para->setValue(2);
869 para->setUnits("mml");
870
871
872 //---------------------------------------------------------------------------
873 //
874 // Creates Reaction objects inside the Model object.
875 //
876 //---------------------------------------------------------------------------
877
878 // Temporary pointers.
879
880 Reaction* reaction;
881 SpeciesReference* spr;
882 KineticLaw* kl;
883
884 //---------------------------------------------------------------------------
885 // (Reaction1) Creates a Reaction object ("v1").
886 //---------------------------------------------------------------------------
887
888 reaction = model->createReaction();
889 reaction->setId("v1");
890
891 //---------------------------------------------------------------------------
892 // Creates Reactant objects inside the Reaction object ("v1").
893 //---------------------------------------------------------------------------
894
895 // (Reactant1) Creates a Reactant object that references Species "x0"
896 // in the model.
897
898 spr = reaction->createReactant();
899 spr->setSpecies("x0");
900
901 //---------------------------------------------------------------------------
902 // Creates a Product object inside the Reaction object ("v1").
903 //---------------------------------------------------------------------------
904
905 // Creates a Product object that references Species "s1" in the model.
906
907 spr = reaction->createProduct();
908 spr->setSpecies("s1");
909
910 //---------------------------------------------------------------------------
911 // Creates a KineticLaw object inside the Reaction object ("v1").
912 //---------------------------------------------------------------------------
913
914 kl = reaction->createKineticLaw();
915
916 // Creates a <notes> element in the KineticLaw object.
917 // Here we illustrate how to do it using a literal string. This requires
918 // known the required syntax of XHTML and the requirements for SBML <notes>
919 // elements. Later below, we show how to create notes using objects instead
920 // of strings.
921
922 string notesString = "<xhtml:p> ((vm * s1)/(km + s1)) * cell </xhtml:p>";
923 kl->setNotes(notesString);
924
925 //---------------------------------------------------------------------------
926 // Creates an ASTNode object which represents the following KineticLaw object.
927 //
928 // <math xmlns=\"http://www.w3.org/1998/Math/MathML\">
929 // <apply>
930 // <times/>
931 // <apply>
932 // <divide/>
933 // <apply>
934 // <times/>
935 // <ci> vm </ci>
936 // <ci> s1 </ci>
937 // </apply>
938 // <apply>
939 // <plus/>
940 // <ci> km </ci>
941 // <ci> s1 </ci>
942 // </apply>
943 // </apply>
944 // <ci> cell </ci>
945 // </apply>
946 // </math>
947 //---------------------------------------------------------------------------
948
949 //
950 // In the following code, ASTNode objects, which construct an ASTNode tree
951 // of the above math, are created and added in the order of preorder traversal
952 // of the tree (i.e. the order corresponds to the nested structure of the above
953 // MathML elements), and thus the following code maybe a bit more efficient but
954 // maybe a bit difficult to read.
955 //
956
957 ASTNode* astMath = new ASTNode(AST_TIMES);
958
959 astMath->addChild(new ASTNode(AST_DIVIDE));
960 ASTNode* astDivide = astMath->getLeftChild();
961
962 astDivide->addChild(new ASTNode(AST_TIMES));
963 ASTNode* astTimes = astDivide->getLeftChild();
964
965 astTimes->addChild(new ASTNode(AST_NAME));
966 astTimes->getLeftChild()->setName("vm");
967
968 astTimes->addChild(new ASTNode(AST_NAME));
969 astTimes->getRightChild()->setName("s1");
970
971 astDivide->addChild(new ASTNode(AST_PLUS));
972 ASTNode* astPlus = astDivide->getRightChild();
973
974 astPlus->addChild(new ASTNode(AST_NAME));
975 astPlus->getLeftChild()->setName("km");
976
977 astPlus->addChild(new ASTNode(AST_NAME));
978 astPlus->getRightChild()->setName("s1");
979
980
981 astMath->addChild(new ASTNode(AST_NAME));
982 astMath->getRightChild()->setName("cell");
983
984 //---------------------------------------------
985 //
986 // set the Math element
987 //
988 //------------------------------------------------
989
990 kl->setMath(astMath);
991 delete astMath;
992
993
994 //---------------------------------------------------------------------------
995 // (Reaction2) Creates a Reaction object ("v2").
996 //---------------------------------------------------------------------------
997
998 reaction = model->createReaction();
999 reaction->setId("v2");
1000
1001 //---------------------------------------------------------------------------
1002 // Creates Reactant objects inside the Reaction object ("v2").
1003 //---------------------------------------------------------------------------
1004
1005 // (Reactant2) Creates a Reactant object that references Species "s1"
1006 // in the model.
1007
1008 spr = reaction->createReactant();
1009 spr->setSpecies("s1");
1010
1011 //---------------------------------------------------------------------------
1012 // Creates a Product object inside the Reaction object ("v2").
1013 //---------------------------------------------------------------------------
1014
1015 // Creates a Product object that references Species "s2" in the model.
1016
1017 spr = reaction->createProduct();
1018 spr->setSpecies("s2");
1019
1020 //---------------------------------------------------------------------------
1021 // Creates a KineticLaw object inside the Reaction object ("v2").
1022 //---------------------------------------------------------------------------
1023
1024 kl = reaction->createKineticLaw();
1025
1026 // Sets a notes (by XMLNode) to the KineticLaw object.
1027 //
1028 // The following code is an alternative to using setNotes(const string&).
1029 // The equivalent code would be like this:
1030 //
1031 // notesString = "<xhtml:p>((vm * s2)/(km + s2))*cell</xhtml:p>";
1032 // kl->setNotes(notesString);
1033
1034 // Creates an XMLNode of start element (<xhtml:p>) without attributes.
1035
1036 XMLNode notesXMLNode(XMLTriple("p", "", "xhtml"), XMLAttributes());
1037
1038 // Adds a text element to the start element.
1039
1040 notesXMLNode.addChild(XMLNode(" ((vm * s2)/(km + s2)) * cell "));
1041
1042 // Adds it to the kineticLaw object.
1043
1044 kl->setNotes(¬esXMLNode);
1045
1046 //---------------------------------------------------------------------------
1047 // Sets a math (ASTNode object) to the KineticLaw object.
1048 //---------------------------------------------------------------------------
1049
1050 // To create mathematical expressions, one would typically construct
1051 // an ASTNode tree as the above example code which creates a math of another
1052 // KineticLaw object. Here, to save some space and illustrate another approach
1053 // of doing it, we will write out the formula in MathML form and then use a
1054 // libSBML convenience function to create the ASTNode tree for us.
1055 // (This is a bit dangerous; it's very easy to make mistakes when writing MathML
1056 // by hand, so in a real program, we would not really want to do it this way.)
1057
1058 string mathXMLString = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">"
1059 " <apply>"
1060 " <times/>"
1061 " <apply>"
1062 " <divide/>"
1063 " <apply>"
1064 " <times/>"
1065 " <ci> vm </ci>"
1066 " <ci> s2 </ci>"
1067 " </apply>"
1068 " <apply>"
1069 " <plus/>"
1070 " <ci> km </ci>"
1071 " <ci> s2 </ci>"
1072 " </apply>"
1073 " </apply>"
1074 " <ci> cell </ci>"
1075 " </apply>"
1076 "</math>";
1077
1078 astMath = readMathMLFromString(mathXMLString.c_str());
1079 kl->setMath(astMath);
1080 delete astMath;
1081
1082
1083 //---------------------------------------------------------------------------
1084 // (Reaction3) Creates a Reaction object ("v3").
1085 //---------------------------------------------------------------------------
1086
1087 reaction = model->createReaction();
1088 reaction->setId("v3");
1089
1090 //---------------------------------------------------------------------------
1091 // Creates Reactant objects inside the Reaction object ("v3").
1092 //---------------------------------------------------------------------------
1093
1094 // (Reactant2) Creates a Reactant object that references Species "s2"
1095 // in the model.
1096
1097 spr = reaction->createReactant();
1098 spr->setSpecies("s2");
1099
1100 //---------------------------------------------------------------------------
1101 // Creates a Product object inside the Reaction object ("v3").
1102 //---------------------------------------------------------------------------
1103
1104 // Creates a Product object that references Species "x1" in the model.
1105
1106 spr = reaction->createProduct();
1107 spr->setSpecies("x1");
1108
1109
1110 //---------------------------------------------------------------------------
1111 // Creates a KineticLaw object inside the Reaction object ("v3").
1112 //---------------------------------------------------------------------------
1113
1114 kl = reaction->createKineticLaw();
1115
1116 // Sets a notes (by string) to the KineticLaw object.
1117
1118 notesString = "<xhtml:p> ((vm * x1)/(km + x1)) * cell </xhtml:p>";
1119 kl->setNotes(notesString);
1120
1121 //---------------------------------------------------------------------------
1122 // Sets a math (ASTNode object) to the KineticLaw object.
1123 //---------------------------------------------------------------------------
1124
1125 mathXMLString = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">"
1126 " <apply>"
1127 " <times/>"
1128 " <apply>"
1129 " <divide/>"
1130 " <apply>"
1131 " <times/>"
1132 " <ci> vm </ci>"
1133 " <ci> x1 </ci>"
1134 " </apply>"
1135 " <apply>"
1136 " <plus/>"
1137 " <ci> km </ci>"
1138 " <ci> x1 </ci>"
1139 " </apply>"
1140 " </apply>"
1141 " <ci> cell </ci>"
1142 " </apply>"
1143 "</math>";
1144
1145 astMath = readMathMLFromString(mathXMLString.c_str());
1146 kl->setMath(astMath);
1147 delete astMath;
1148
1149
1150 // Returns the created SBMLDocument object.
1151 // The returned object must be explicitly deleted by the caller,
1152 // otherwise memory leak will happen.
1153
1154 return sbmlDoc;
1155
1156 }
1157
1158
1159
1160 /**
1161 *
1162 * Creates an SBML model represented in "7.8 Example involving function definitions"
1163 * in the SBML Level 2 Version 4 Specification.
1164 *
1165 */
createExampleInvolvingFunctionDefinitions()1166 SBMLDocument* createExampleInvolvingFunctionDefinitions()
1167 {
1168 const unsigned int level = Level;
1169 const unsigned int version = Version;
1170
1171 //---------------------------------------------------------------------------
1172 //
1173 // Creates an SBMLDocument object
1174 //
1175 //---------------------------------------------------------------------------
1176
1177 SBMLDocument* sbmlDoc = new SBMLDocument(level,version);
1178
1179 //---------------------------------------------------------------------------
1180 //
1181 // Creates a Model object inside the SBMLDocument object.
1182 //
1183 //---------------------------------------------------------------------------
1184
1185 Model* model = sbmlDoc->createModel();
1186 model->setId("functionExample");
1187
1188 //---------------------------------------------------------------------------
1189 //
1190 // Creates a FunctionDefinition object inside the Model object.
1191 //
1192 //---------------------------------------------------------------------------
1193
1194 FunctionDefinition* fdef = model->createFunctionDefinition();
1195 fdef->setId("f");
1196
1197 // Sets a math (ASTNode object) to the FunctionDefinition object.
1198
1199 string mathXMLString = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">"
1200 " <lambda>"
1201 " <bvar>"
1202 " <ci> x </ci>"
1203 " </bvar>"
1204 " <apply>"
1205 " <times/>"
1206 " <ci> x </ci>"
1207 " <cn> 2 </cn>"
1208 " </apply>"
1209 " </lambda>"
1210 "</math>";
1211
1212 ASTNode* astMath = readMathMLFromString(mathXMLString.c_str());
1213 fdef->setMath(astMath);
1214 delete astMath;
1215
1216
1217 //---------------------------------------------------------------------------
1218 //
1219 // Creates a Compartment object inside the Model object.
1220 //
1221 //---------------------------------------------------------------------------
1222
1223 Compartment* comp;
1224 const string compName = "compartmentOne";
1225
1226 // Creates a Compartment object ("compartmentOne")
1227
1228 comp = model->createCompartment();
1229 comp->setId(compName);
1230
1231 // Sets the "size" attribute of the Compartment object.
1232 //
1233 // The units of this Compartment object is the default SBML
1234 // units of volume (litre), and thus we don't have to explicitly invoke
1235 // setUnits("litre") function to set the default units.
1236 //
1237 comp->setSize(1);
1238
1239
1240 //---------------------------------------------------------------------------
1241 //
1242 // Creates Species objects inside the Model object.
1243 //
1244 //---------------------------------------------------------------------------
1245
1246 Species* sp;
1247
1248 //---------------------------------------------------------------------------
1249 // (Species1) Creates a Species object ("S1")
1250 //---------------------------------------------------------------------------
1251
1252 sp = model->createSpecies();
1253 sp->setId("S1");
1254
1255 // Sets the "compartment" attribute of the Species object to identify the
1256 // compartnet in which the Species object located.
1257
1258 sp->setCompartment(compName);
1259
1260 // Sets the "initialConcentration" attribute of the Species object.
1261 //
1262 // The units of this Species object is determined by two attributes of this
1263 // Species object ("substanceUnits" and "hasOnlySubstanceUnits") and the
1264 // "spatialDimension" attribute of the Compartment object ("cytosol") in which
1265 // this species object located.
1266 // Since the default values are used for "substanceUnits" (substance (mole))
1267 // and "hasOnlySubstanceUnits" (false) and the value of "spatialDimension" (3)
1268 // is greater than 0, the units of this Species object is mole/litre .
1269 //
1270
1271 sp->setInitialConcentration(1);
1272
1273 //---------------------------------------------------------------------------
1274 // (Species2) Creates a Species object ("S2")
1275 //---------------------------------------------------------------------------
1276
1277 sp = model->createSpecies();
1278 sp->setId("S2");
1279 sp->setCompartment(compName);
1280 sp->setInitialConcentration(0);
1281
1282
1283 //---------------------------------------------------------------------------
1284 //
1285 // Creates a global Parameter object inside the Model object.
1286 //
1287 //---------------------------------------------------------------------------
1288
1289 Parameter* para;
1290
1291 // Creates a Parameter ("t")
1292
1293 para = model->createParameter();
1294 para->setId("t");
1295 para->setValue(1);
1296 para->setUnits("second");
1297
1298
1299 //---------------------------------------------------------------------------
1300 //
1301 // Creates Reaction objects inside the Model object.
1302 //
1303 //---------------------------------------------------------------------------
1304
1305 // Temporary pointers.
1306
1307 Reaction* reaction;
1308 SpeciesReference* spr;
1309 KineticLaw* kl;
1310
1311 //---------------------------------------------------------------------------
1312 // (Reaction1) Creates a Reaction object ("reaction_1").
1313 //---------------------------------------------------------------------------
1314
1315 reaction = model->createReaction();
1316 reaction->setId("reaction_1");
1317 reaction->setReversible(false);
1318
1319 //---------------------------------------------------------------------------
1320 // Creates Reactant objects inside the Reaction object ("reaction_1").
1321 //---------------------------------------------------------------------------
1322
1323 // (Reactant1) Creates a Reactant object that references Species "S1"
1324 // in the model.
1325
1326 spr = reaction->createReactant();
1327 spr->setSpecies("S1");
1328
1329 //---------------------------------------------------------------------------
1330 // Creates a Product object inside the Reaction object ("reaction_1").
1331 //---------------------------------------------------------------------------
1332
1333 // Creates a Product object that references Species "S2" in the model.
1334
1335 spr = reaction->createProduct();
1336 spr->setSpecies("S2");
1337
1338
1339 //---------------------------------------------------------------------------
1340 // Creates a KineticLaw object inside the Reaction object ("reaction_1").
1341 //---------------------------------------------------------------------------
1342
1343 kl = reaction->createKineticLaw();
1344
1345 //---------------------------------------------------------------------------
1346 // Sets a math (ASTNode object) to the KineticLaw object.
1347 //---------------------------------------------------------------------------
1348
1349 mathXMLString = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">"
1350 " <apply>"
1351 " <divide/>"
1352 " <apply>"
1353 " <times/>"
1354 " <apply>"
1355 " <ci> f </ci>"
1356 " <ci> S1 </ci>"
1357 " </apply>"
1358 " <ci> compartmentOne </ci>"
1359 " </apply>"
1360 " <ci> t </ci>"
1361 " </apply>"
1362 "</math>";
1363
1364 astMath = readMathMLFromString(mathXMLString.c_str());
1365 kl->setMath(astMath);
1366 delete astMath;
1367
1368
1369 // Returns the created SBMLDocument object.
1370 // The returned object must be explicitly deleted by the caller,
1371 // otherwise memory leak will happen.
1372
1373 return sbmlDoc;
1374 }
1375
1376
1377 //===============================================================================
1378 //
1379 //
1380 // Helper functions for writing/validating the given SBML documents.
1381 //
1382 //
1383 //===============================================================================
1384
1385 /**
1386 *
1387 * Validates the given SBMLDocument.
1388 *
1389 * This function is based on validateSBML.cpp implemented by
1390 * Sarah Keating, Ben Bornstein, and Michael Hucka.
1391 *
1392 */
validateExampleSBML(SBMLDocument * sbmlDoc)1393 bool validateExampleSBML (SBMLDocument* sbmlDoc)
1394 {
1395 if (!sbmlDoc)
1396 {
1397 cerr << "validateExampleSBML: given a null SBML Document" << endl;
1398 return false;
1399 }
1400
1401 string consistencyMessages;
1402 string validationMessages;
1403 bool noProblems = true;
1404 unsigned int numCheckFailures = 0;
1405 unsigned int numConsistencyErrors = 0;
1406 unsigned int numConsistencyWarnings = 0;
1407 unsigned int numValidationErrors = 0;
1408 unsigned int numValidationWarnings = 0;
1409
1410 // LibSBML 3.3 is lenient when generating models from scratch using the
1411 // API for creating objects. Once the whole model is done and before it
1412 // gets written out, it's important to check that the whole model is in
1413 // fact complete, consistent and valid.
1414
1415 numCheckFailures = sbmlDoc->checkInternalConsistency();
1416 if ( numCheckFailures > 0 )
1417 {
1418 noProblems = false;
1419 for (unsigned int i = 0; i < numCheckFailures; i++)
1420 {
1421 const SBMLError* sbmlErr = sbmlDoc->getError(i);
1422 if ( sbmlErr->isFatal() || sbmlErr->isError() )
1423 {
1424 ++numConsistencyErrors;
1425 }
1426 else
1427 {
1428 ++numConsistencyWarnings;
1429 }
1430 }
1431 ostringstream oss;
1432 sbmlDoc->printErrors(oss);
1433 consistencyMessages = oss.str();
1434 }
1435
1436 // If the internal checks fail, it makes little sense to attempt
1437 // further validation, because the model may be too compromised to
1438 // be properly interpreted.
1439
1440 if (numConsistencyErrors > 0)
1441 {
1442 consistencyMessages += "Further validation aborted.";
1443 }
1444 else
1445 {
1446 numCheckFailures = sbmlDoc->checkConsistency();
1447 if ( numCheckFailures > 0 )
1448 {
1449 noProblems = false;
1450 for (unsigned int i = 0; i < numCheckFailures; i++)
1451 {
1452 const SBMLError* sbmlErr = sbmlDoc->getError(i);
1453 if ( sbmlErr->isFatal() || sbmlErr->isError() )
1454 {
1455 ++numValidationErrors;
1456 }
1457 else
1458 {
1459 ++numValidationWarnings;
1460 }
1461 }
1462 ostringstream oss;
1463 sbmlDoc->printErrors(oss);
1464 validationMessages = oss.str();
1465 }
1466 }
1467
1468 if (noProblems)
1469 return true;
1470 else
1471 {
1472 if (numConsistencyErrors > 0)
1473 {
1474 cout << "ERROR: encountered " << numConsistencyErrors
1475 << " consistency error" << (numConsistencyErrors == 1 ? "" : "s")
1476 << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl;
1477 }
1478 if (numConsistencyWarnings > 0)
1479 {
1480 cout << "Notice: encountered " << numConsistencyWarnings
1481 << " consistency warning" << (numConsistencyWarnings == 1 ? "" : "s")
1482 << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl;
1483 }
1484 cout << endl << consistencyMessages;
1485
1486 if (numValidationErrors > 0)
1487 {
1488 cout << "ERROR: encountered " << numValidationErrors
1489 << " validation error" << (numValidationErrors == 1 ? "" : "s")
1490 << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl;
1491 }
1492 if (numValidationWarnings > 0)
1493 {
1494 cout << "Notice: encountered " << numValidationWarnings
1495 << " validation warning" << (numValidationWarnings == 1 ? "" : "s")
1496 << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl;
1497 }
1498 cout << endl << validationMessages;
1499
1500 return (numConsistencyErrors == 0 && numValidationErrors == 0);
1501 }
1502 }
1503
1504
1505 /**
1506 *
1507 * Writes the given SBMLDocument to the given file.
1508 *
1509 */
writeExampleSBML(const SBMLDocument * sbmlDoc,const string & filename)1510 bool writeExampleSBML(const SBMLDocument* sbmlDoc, const string& filename)
1511 {
1512 SBMLWriter sbmlWriter;
1513
1514 bool result = sbmlWriter.writeSBML(sbmlDoc, filename);
1515
1516 if (result)
1517 {
1518 cout << "Wrote file \"" << filename << "\"" << endl;
1519 return true;
1520 }
1521 else
1522 {
1523 cerr << "Failed to write \"" << filename << "\"" << endl;
1524 return false;
1525 }
1526 }
1527