From 67b2eefc3390aa6fb1deaccd7f6a3fd50cdd4de6 Mon Sep 17 00:00:00 2001
From: Sarah Keating <s.keating@ucl.ac.uk>
Date: Sun, 8 Dec 2024 14:47:16 +0000
Subject: [PATCH] adding these for now but not quite finished

---
 src/sbml/conversion/ExpressionAnalyser.cpp    | 911 ++++++++++++++----
 src/sbml/conversion/ExpressionAnalyser.h      |  78 +-
 src/sbml/conversion/SBMLRateRuleConverter.cpp |  57 +-
 src/sbml/conversion/SBMLRateRuleConverter.h   |   6 +-
 .../test/TestExpressionAnalyser.cpp           | 555 +++++++++++
 .../test/TestSBMLRateRuleConverter.cpp        |  68 +-
 .../conversion/test/test-data/mraterules5.xml |   7 +-
 .../conversion/test/test-data/mraterules7.xml |  36 +
 src/sbml/conversion/test/test-data/test       | 122 +++
 .../conversion/test/test-data/valid_54_bio.c  | 103 ++
 .../conversion/test/test-data/valid_54_rr.c   |  79 ++
 11 files changed, 1785 insertions(+), 237 deletions(-)
 create mode 100644 src/sbml/conversion/test/TestExpressionAnalyser.cpp
 create mode 100644 src/sbml/conversion/test/test-data/mraterules7.xml
 create mode 100644 src/sbml/conversion/test/test-data/test
 create mode 100644 src/sbml/conversion/test/test-data/valid_54_bio.c
 create mode 100644 src/sbml/conversion/test/test-data/valid_54_rr.c

diff --git a/src/sbml/conversion/ExpressionAnalyser.cpp b/src/sbml/conversion/ExpressionAnalyser.cpp
index b548e67835..a8a4541476 100644
--- a/src/sbml/conversion/ExpressionAnalyser.cpp
+++ b/src/sbml/conversion/ExpressionAnalyser.cpp
@@ -6,7 +6,7 @@
  *
  * <!--------------------------------------------------------------------------
  * This file is part of libSBML.  Please visit http://sbml.org for more
- * information about SBML, and the latest version of libSBML.
+
  *
  * Copyright (C) 2013-2018 jointly by the following organizations:
  *     1. California Institute of Technology, Pasadena, CA, USA
@@ -40,6 +40,7 @@
 #include <string>
 #include <vector>
 #include <map>
+#include <iostream>
 #include <sbml/SBMLTypes.h>
 #include <sbml/math/ASTNodeType.h>
 #include <sbml/conversion/SBMLRateRuleConverter.h>
@@ -55,7 +56,8 @@ ExpressionAnalyser::ExpressionAnalyser()
     : mModel (NULL), 
       mODEs (),
     mNewVarName("newVar"),
-    mNewVarCount(1)
+    mNewVarCount(1),
+    mHiddenSpecies (NULL)
 {
 }
 
@@ -65,7 +67,8 @@ ExpressionAnalyser::ExpressionAnalyser(Model * m, pairODEs odes)
     : mModel(m),
     mODEs(odes),
     mNewVarName("newVar"),
-    mNewVarCount(1)
+    mNewVarCount(1),
+    mHiddenSpecies(NULL) 
 {
   SBMLTransforms::mapComponentValues(mModel);
   mModel->populateAllElementIdList();
@@ -75,7 +78,8 @@ ExpressionAnalyser::ExpressionAnalyser(const ExpressionAnalyser& orig) :
   mModel( orig.mModel),
     mODEs(orig.mODEs),
     mNewVarName(orig.mNewVarName),
-    mNewVarCount(orig.mNewVarCount)
+    mNewVarCount(orig.mNewVarCount),
+    mHiddenSpecies(orig.mHiddenSpecies)
 {
     SBMLTransforms::mapComponentValues(mModel);
     mModel->populateAllElementIdList();
@@ -93,6 +97,7 @@ ExpressionAnalyser::operator=(const ExpressionAnalyser& rhs)
     mODEs = rhs.mODEs;
     mNewVarName = rhs.mNewVarName;
     mNewVarCount = rhs.mNewVarCount;
+    mHiddenSpecies = rhs.mHiddenSpecies;
   }
   SBMLTransforms::mapComponentValues(mModel);
   mModel->populateAllElementIdList();
@@ -122,6 +127,7 @@ ExpressionAnalyser::~ExpressionAnalyser ()
   }
   mODEs.clear();
   SBMLTransforms::clearComponentValues(mModel);
+  mHiddenSpecies = NULL;
 }
 
 /*
@@ -147,6 +153,93 @@ ExpressionAnalyser::setModel(Model* model)
   return LIBSBML_OPERATION_SUCCESS;
 }
 
+void ExpressionAnalyser::substituteParametersForExpressions(List* hiddenSpecies)
+{
+    if (hiddenSpecies == NULL || mExpressions.empty())
+        return;
+
+    // need to actually address the expressions in the correct order
+   // replace k-x-y first with newParam type=TYPE_K_MINUS_X_MINUS_Y
+   // and then k+v-x-y with newParam+v  type=TYPE_K_PLUS_V_MINUS_X_MINUS_Y
+   // and then k-x+w-y with newParam+w  type=TYPE_K_MINUS_X_PLUS_W_MINUS_Y
+    // and then k-x with newParam2      type=TYPE_K_MINUS_X
+    // and then k+v-x with newParam2+v  type=TYPE_K_PLUS_V_MINUS_X  
+
+
+    // need to go through the expressions and make sure that they are substituting the correct thing
+    // this needs to be done before any parameters are substituted I think
+
+    for (unsigned int j = 0; j < mExpressions.size(); j++)
+    {
+        SubstitutionValues_t* exp = mExpressions[j];
+        if (exp->type == TYPE_K_MINUS_X_MINUS_Y && exp->z_value.empty())
+        {
+            std::string zName = getUniqueNewParameterName();
+            exp->z_value = zName;
+            mNewVarCount++;
+        }
+        else if (j > 0 &&
+            exp->type == TYPE_K_PLUS_V_MINUS_X_MINUS_Y &&
+            exp->z_expression == NULL)
+        {
+            // we have a value for  k-x-y
+            ASTNode* replacement = new ASTNode(AST_PLUS);
+            ASTNode* z = new ASTNode(AST_NAME);
+            z->setName(mExpressions[j - 1]->z_value.c_str());
+            ASTNode* v = exp->v_expression->deepCopy();
+            replacement->addChild(z);
+            replacement->addChild(v);
+            exp->z_expression = replacement->deepCopy();
+
+        }
+    }
+}
+
+void ExpressionAnalyser::substituteParameters(List* hiddenSpecies, SubstitutionValues_t* exp)
+{
+    for (unsigned int j = 0; j < mODEs.size(); j++)
+    {
+        std::pair<std::string, ASTNode*> ode = mODEs.at(j);
+        ASTNode* odeRHS = ode.second;
+        bool index = isParameterAlreadyCreated(exp->z_value);
+        if (index == true)
+        {
+            replaceExpressionWithNewParameter(odeRHS, exp);
+        }
+        else
+        {
+            std::string zName = getUniqueNewParameterName();
+            exp->z_value = zName;
+            mNewVarCount++;
+//            replaceExpressionWithNewParameter(odeRHS, exp);
+//            addParametersAndRateRules(hiddenSpecies, exp);
+        }
+        //cout << "ode in main: " << SBML_formulaToL3String(odeRHS) << endl;
+    }
+}
+
+SubstitutionValues_t* ExpressionAnalyser::getSubstitutionValuesByType(ExpressionType_t type,
+    size_t index)
+{
+    SubstitutionValues_t* exp = NULL;
+
+    while (exp == NULL && index < mExpressions.size())
+    {
+        exp = mExpressions[index];
+        if (exp->type == type)
+        {
+            index++;
+            return exp;
+        }
+        else
+        {
+            exp = NULL;
+            index++;
+        }
+    }
+    return NULL;
+}
+
 /*
 * Check whether two SubstitutionValues_t are identical.
  */
@@ -206,6 +299,177 @@ ExpressionAnalyser::areIdenticalSubstitutionValues(SubstitutionValues_t* values1
     }
 }
 
+/*
+* Check whether two SubstitutionValues_t match the values that we expect if we need to add them
+* based on the type of the SubstitutionValues_t
+ */
+bool ExpressionAnalyser::expressionExists(SubstitutionValues_t* current, 
+    SubstitutionValues_t* mightAdd)
+{
+    bool alreadyExists = false;
+    // all expressions will have K and X
+    alreadyExists = matchesK(current, mightAdd) &&
+        matchesXValue(current, mightAdd) &&
+        matchesDxdtExpression(current, mightAdd);
+
+    //switch (current->type)
+    //{
+    //case TYPE_K_MINUS_X_MINUS_Y:
+    //    switch (mightAdd->type)
+    //    {
+    //    case TYPE_K_MINUS_X_MINUS_Y:
+    //        alreadyExists = matchesKValue(current, mightAdd) && 
+    //            matchesXValue(current, mightAdd) && 
+    //            matchesYValue(current, mightAdd) && 
+    //            matchesDxdtExpression(current, mightAdd) && 
+    //            matchesDydtExpression(current, mightAdd) ;
+    //        break;
+    //    case TYPE_K_PLUS_V_MINUS_X_MINUS_Y:
+    //        alreadyExists = matchesKValue(current, mightAdd) && 
+    //            matchesVExpression(current, mightAdd) && 
+    //            matchesXValue(current, mightAdd) && 
+    //            matchesYValue(current, mightAdd);
+    //        break;
+    //    case TYPE_K_MINUS_X_PLUS_W_MINUS_Y:
+    //        alreadyExists = matchesKValue(current, mightAdd) && 
+    //            matchesWExpression(current, mightAdd) && 
+    //            matchesXValue(current, mightAdd) && 
+    //            matchesYValue(current, mightAdd);
+    //        break;
+    //    case TYPE_K_MINUS_X:
+    //        alreadyExists = matchesKValue(current, mightAdd) && 
+    //            matchesXValue(current, mightAdd) &&
+    //            matchesDxdtExpression(current, mightAdd);
+    //         break;
+    //    case TYPE_K_PLUS_V_MINUS_X:
+    //        alreadyExists = matchesKValue(current, mightAdd) && 
+    //            matchesVExpression(current, mightAdd) && 
+    //            matchesXValue(current, mightAdd);
+    //        break;
+    //    default:
+    //        break;
+    //    }
+    //    break;
+    //case TYPE_K_PLUS_V_MINUS_X_MINUS_Y:
+    //    switch (mightAdd->type)
+    //    { // if we already have k+v-x-y then we need to record that we have k-x-y as well
+    //        // so we only record a match value if one of the values doesn't match
+    //    case TYPE_K_MINUS_X_MINUS_Y:
+    //        alreadyExists = !matchesKValue(current, mightAdd) ||
+    //            !matchesXValue(current, mightAdd) ||
+    //            !matchesYValue(current, mightAdd) ||
+    //            !matchesDxdtExpression(current, mightAdd) ||
+    //            !matchesDydtExpression(current, mightAdd);
+    //        break;
+    //    case TYPE_K_PLUS_V_MINUS_X_MINUS_Y:
+    //        alreadyExists = matchesKValue(current, mightAdd) && 
+    //            matchesVExpression(current, mightAdd) && 
+    //            matchesXValue(current, mightAdd) && 
+    //            matchesYValue(current, mightAdd);
+    //        break;
+    //    case TYPE_K_MINUS_X_PLUS_W_MINUS_Y:
+    //        alreadyExists = matchesKValue(current, mightAdd) && 
+    //            matchesWExpression(current, mightAdd) && 
+    //            matchesXValue(current, mightAdd) &&
+    //            matchesYValue(current, mightAdd);
+    //        break;
+    //    case TYPE_K_MINUS_X:
+    //        alreadyExists = matchesKValue(current, mightAdd) &&
+    //            matchesXValue(current, mightAdd) &&
+    //            matchesDxdtExpression(current, mightAdd);
+    //        break;
+    //    case TYPE_K_PLUS_V_MINUS_X:
+    //        alreadyExists = matchesKValue(current, mightAdd) &&
+    //            matchesXValue(current, mightAdd) &&
+    //            matchesDxdtExpression(current, mightAdd);
+    //        break;
+    //    default:
+    //        break;
+    //    }
+    //    break;
+    //case TYPE_K_MINUS_X_PLUS_W_MINUS_Y:
+    //    printSubstitutionValues(current);
+    //    printSubstitutionValues(mightAdd);
+    //    alreadyExists = matchesKValue(current, mightAdd) && 
+    //        matchesWExpression(current, mightAdd) && 
+    //        matchesXValue(current, mightAdd) && 
+    //        matchesDxdtExpression(current, mightAdd) &&
+    //        matchesYValue(current, mightAdd) &&
+    //        matchesDydtExpression(current, mightAdd);
+    //    break;
+    //case TYPE_K_MINUS_X:
+    //    alreadyExists = matchesKValue(current, mightAdd) && 
+    //        matchesXValue(current, mightAdd);
+    //    break;
+    //case TYPE_K_PLUS_V_MINUS_X:
+    //    alreadyExists = matchesKValue(current, mightAdd) && 
+    //        matchesVExpression(current, mightAdd) && 
+    //        matchesXValue(current, mightAdd);
+    //    break;
+    //default:
+    //    break;
+    //}
+
+    return alreadyExists;
+} 
+
+bool ExpressionAnalyser::matchesK(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return matchesKValue(values1, values2) || matchesKRealValue(values1, values2);
+}
+
+bool ExpressionAnalyser::matchesKValue(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return values1->k_value == values2->k_value;
+}
+
+bool ExpressionAnalyser::matchesKRealValue(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return ((util_isNaN(values1->k_real_value) && util_isNaN(values2->k_real_value)) ||
+        util_isEqual(values1->k_real_value, values2->k_real_value));
+}
+
+bool ExpressionAnalyser::matchesXValue(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return values1->x_value == values2->x_value;
+}
+
+bool ExpressionAnalyser::matchesYValue(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return values1->y_value == values2->y_value;
+}
+
+bool ExpressionAnalyser::matchesVExpression(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return (values1->v_expression != NULL && values2->v_expression != NULL &&
+        values1->v_expression->exactlyEqual(*(values2->v_expression)) == true);
+}
+
+bool ExpressionAnalyser::matchesWExpression(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return (values1->w_expression != NULL && values2->w_expression != NULL &&
+        values1->w_expression->exactlyEqual(*(values2->w_expression)) == true);
+}
+
+bool ExpressionAnalyser::matchesDxdtExpression(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return (values1->dxdt_expression != NULL && values2->dxdt_expression != NULL &&
+        values1->dxdt_expression->exactlyEqual(*(values2->dxdt_expression)) == true);
+}
+
+bool ExpressionAnalyser::matchesDydtExpression(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return (values1->dydt_expression != NULL && values2->dydt_expression != NULL &&
+        values1->dydt_expression->exactlyEqual(*(values2->dydt_expression)) == true);
+}
+
+bool ExpressionAnalyser::matchesCurrentNode(SubstitutionValues_t* values1, SubstitutionValues_t* values2)
+{
+    return (values1->current != NULL && values2->current != NULL &&
+        values1->current->exactlyEqual(*(values2->current)) == true);
+}
+
+
 void
 ExpressionAnalyser::printSubstitutionValues(const SubstitutionValues_t* values)
 {
@@ -230,6 +494,10 @@ ExpressionAnalyser::printSubstitutionValues(const SubstitutionValues_t* values)
             cout << "w_expression: " << SBML_formulaToL3String(values->w_expression) << " " << values->w_expression << endl;
         else
             cout << "w_expression: NULL" << endl;
+        if (values->z_expression != NULL)
+            cout << "z_expression: " << SBML_formulaToL3String(values->z_expression) << " " << values->z_expression << endl;
+        else
+            cout << "z_expression: NULL" << endl;
         cout << "type: " << values->type << endl;
         if (values->current != NULL) 
             cout << "current: " << SBML_formulaToL3String(values->current) << " " << values->current << endl;
@@ -238,7 +506,39 @@ ExpressionAnalyser::printSubstitutionValues(const SubstitutionValues_t* values)
         cout << "odeIndex: " << values->odeIndex << endl;
 }
 
+SubstitutionValues_t* ExpressionAnalyser::createBlankSubstitutionValues()
+{
+    SubstitutionValues_t* values = new SubstitutionValues_t;
+    values->type = TYPE_UNKNOWN;
+    values->k_real_value = util_NaN();
+    values->dxdt_expression = NULL;
+    values->dydt_expression = NULL;
+    values->v_expression = NULL;
+    values->w_expression = NULL;
+    values->z_expression = NULL;
+    values->k_value = "";
+    values->x_value = "";
+    values->y_value = "";
+    values->z_value = "";
+    values->odeIndex = 0;
+    values->current = NULL;
+    return values;
+}
+
+unsigned int
+ExpressionAnalyser::getNumExpressions()
+{
+    return mExpressions.size();
+}
 
+SubstitutionValues_t* ExpressionAnalyser::getExpression(unsigned int index)
+{
+    if (index < mExpressions.size())
+    {
+        return mExpressions.at(index);
+    }
+    return nullptr;
+}
 
 /*
 * Check whether the expression has a parent expression which may already have been analysed 
@@ -246,23 +546,15 @@ ExpressionAnalyser::printSubstitutionValues(const SubstitutionValues_t* values)
 * e.g. if we have k-x-y do not need to analyse k-x
 */
 bool
-ExpressionAnalyser::hasExpressionAlreadyRecorded(SubstitutionValues_t* value)
+ExpressionAnalyser::hasExpressionAlreadyBeenRecorded(SubstitutionValues_t* value)
 {
-    bool found = false;
-  for (unsigned int i = mExpressions.size(); i > 0; i--)
-  {
-    SubstitutionValues_t* exp = mExpressions.at(i - 1);
+  bool found = false;
+  size_t size = mExpressions.size();
 
-    // if we already have parent eg K+v-x-y dont record children eg K+v-x
-    std::pair<ASTNode*, int> parent = getParentNode(value->current, exp->current);
-    if (parent.first != NULL)
-    {
-    found = true;
-    }
-    else
-    {
-        found = found || areIdenticalSubstitutionValues(value, exp);
-    }
+  while (size > 0 && !found)
+  {
+      found = expressionExists(mExpressions.at(size - 1), value);
+      size--;
   }
   return found  ;
 }
@@ -271,130 +563,170 @@ ExpressionAnalyser::hasExpressionAlreadyRecorded(SubstitutionValues_t* value)
 bool
 ExpressionAnalyser::analyseNode(ASTNode* node, SubstitutionValues_t *value)
 {
-  unsigned int numChildren = node->getNumChildren();
-  ASTNodeType_t type = node->getType();
-  ASTNode* rightChild = node->getRightChild();
-  ASTNode* leftChild = node->getLeftChild();
-  bool isNumber = false;
-  switch (type)
-  {
-  case AST_PLUS:
-    //  -x+y node binary; plus; left child type minus; rightchild var/const
-    //           +
-    //        -     y
-    //        x
-    if (numChildren != 2 || rightChild->getType() != AST_NAME
-      || leftChild->getType() != AST_MINUS
-      || leftChild->getNumChildren() != 1)
-      return false;
-
-    // if we get to this point, the only thing left to check is 
-    // whether the ->left->right grandchild (the x in -x+y) is a variable species.
-    if (isVariableSpeciesOrParameter(leftChild->getChild(0)))
-    {
-      value->x_value = leftChild->getChild(0)->getName();
-      value->y_value = rightChild->getName();
-      value->dydt_expression = getODEFor(rightChild->getName());
-      value->dxdt_expression = getODEFor(leftChild->getChild(0)->getName());
-      value->type = TYPE_MINUS_X_PLUS_Y;
-      value->current = node;
-      return true;
-    }
-    break;
-
-  case AST_MINUS:
-    //          -                -               
-    //        k   x           -     y    
-    //                      k   x      
-    //  k-x or k-x-y node binary; right child (x,y) is variable
-    //          -                    -                          -
-    //      +      x            -         y              +            y
-    //   k     v            +      x                 -       w
-    //                  k      v                 k       x
-    //  k+v-x; right child x var; left child plus node with left child k constant
-    //  k+v-x-y; rightchild var y; left child minus == k+v-x
-    //  k-x+w-y; rightchild var y; left child plus with left child == k-x 
-    if (numChildren != 2 || !isVariableSpeciesOrParameter(rightChild))
-      return false;
-    // if left child is  numerical constant or a parameter and right child variable, it IS k-x
-    if (isNumericalConstantOrConstantParameter(leftChild, isNumber)
-      && isVariableSpeciesOrParameter(rightChild))
-    {
-     
-      if (isNumber)
-      {
-          value->k_value = "number";
-          value->k_real_value = leftChild->getValue();
-      }
-      else
-      {
-          value->k_value = leftChild->getName();
-      }
-      value->x_value = rightChild->getName();
-      value->dxdt_expression = getODEFor(rightChild->getName());
-      value->type = TYPE_K_MINUS_X;
-      value->current = node;
-      return true;
-    }
-    // left child + with it's left child const we have k+v-x
-    // left child + with it's left child k-x we have k-x+w-y
-    else if (leftChild->getType() == AST_PLUS)
+    unsigned int numChildren = node->getNumChildren();
+    ASTNodeType_t type = node->getType();
+    ASTNode* rightChild = node->getRightChild();
+    ASTNode* leftChild = node->getLeftChild();
+      cout << "RIGHT CHILD: " << SBML_formulaToL3String(rightChild) << endl;
+      cout << "left CHILD: " << SBML_formulaToL3String(leftChild) << endl;
+
+    if (isTypeKminusX(numChildren, rightChild, leftChild, type, value) ||
+        isTypeKminusXminusY(numChildren, rightChild, leftChild, type, value) ||
+        isTypeKplusVminusX(numChildren, rightChild, leftChild, type, value) ||
+        isTypeKplusVminusXminusY(numChildren, rightChild, leftChild, type, value) ||
+        isTypeKminusXplusWminusY(numChildren, rightChild, leftChild, type, value))
     {
-
-        // TO DO fix this for k Or w being a number
-      if (isNumericalConstantOrConstantParameter(leftChild->getChild(0), isNumber))
-      {
-        value->k_value = leftChild->getChild(0)->getName();
-        value->x_value = rightChild->getName();
-        value->dxdt_expression = getODEFor(rightChild->getName());
-        value->v_expression = leftChild->getChild(1);
-        value->type = TYPE_K_PLUS_V_MINUS_X;
-        value->current = node;
-        return true;
-      }
-      else if (analyseNode(leftChild->getChild(0), value) && value->type == TYPE_K_MINUS_X)
-      {
-        value->y_value = rightChild->getName();
-        value->dydt_expression = getODEFor(rightChild->getName());
-        value->w_expression = leftChild->getChild(1);
-        value->type = TYPE_K_MINUS_X_PLUS_W_MINUS_Y;
         value->current = node;
         return true;
-      }
-    }
-    else if (leftChild->getType() == AST_MINUS
-      && isVariableSpeciesOrParameter(leftChild->getRightChild()))
-    {
-      // if left child is k+v-x we have k+v-x-y
-      // or if left child is k-x we have k-x-y
-      if (analyseNode(leftChild, value))
-      {
-        if (value->type == TYPE_K_PLUS_V_MINUS_X)
-        {
-          value->type = TYPE_K_PLUS_V_MINUS_X_MINUS_Y;
-          value->y_value = rightChild->getName();
-          value->dydt_expression = getODEFor(value->y_value);
-          value->current = node;
-          return true;
-        }
-        else if (value->type == TYPE_K_MINUS_X)
-        {
-          value->y_value = rightChild->getName();
-          value->dydt_expression = getODEFor(rightChild->getName());
-          value->type = TYPE_K_MINUS_X_MINUS_Y;
-          value->current = node;
-          return true;
-        }
-      }
-      return false;
     }
-    break;
-  default:
     return false;
-  }
-  return false;
 }
 
+
+
+  //switch (type)
+  //{
+  ////case AST_PLUS:
+  //  //  -x+y node binary; plus; left child type minus; rightchild var/const
+  //  //           +
+  //  //        -     y
+  //  //        x
+  //  // we might have encountered w+(k-x)as part of k-x+w-y but it might not be part of it
+  //    if (rightChild->getNumChildren() == 2 && analyseNode(rightChild, value) && value->type == TYPE_K_MINUS_X)
+  //    {
+  //        value->w_expression = leftChild;
+  //        value->type = TYPE_K_MINUS_X_PLUS_W_MINUS_Y;
+  //        value->current = node;
+  //        return true;
+  //    }
+  //  if (numChildren != 2 || rightChild->getType() != AST_NAME
+  //    || leftChild->getType() != AST_MINUS
+  //    || leftChild->getNumChildren() != 1)
+  //    return false;
+
+  //  // if we get to this point, the only thing left to check is 
+  //  // whether the ->left->right grandchild (the x in -x+y) is a variable species.
+  //  if (isVariableSpeciesOrParameter(leftChild->getChild(0)))
+  //  {
+  //    value->x_value = leftChild->getChild(0)->getName();
+  //    value->y_value = rightChild->getName();
+  //    value->dydt_expression = getODEFor(rightChild->getName());
+  //    value->dxdt_expression = getODEFor(leftChild->getChild(0)->getName());
+  //    value->type = TYPE_MINUS_X_PLUS_Y;
+  //    value->current = node;
+  //    return true;
+  //  }
+  //  break;
+
+  //case AST_MINUS:
+  //  //          -                -               
+  //  //        k   x           -     y    
+  //  //                      k   x      
+  //  //  k-x or k-x-y node binary; right child (x,y) is variable
+  //  //          -                    -                          -
+  //  //      +      x            -         y              +            y
+  //  //   k     v            +      x                 -       w
+  //  //                  k      v                 k       x
+  //  //    
+  //  // //  k+v-x; right child x var; left child plus node with left child k constant
+  //  //  k+v-x-y; rightchild var y; left child minus == k+v-x
+  //  //  k-x+w-y; rightchild var y; left child plus with left child == k-x 
+  //  if (numChildren != 2 || !isVariableSpeciesOrParameter(rightChild))
+  //    return false;
+  //  // if left child is  numerical constant or a parameter and right child variable, it IS k-x
+  //  if (isNumericalConstantOrConstantParameter(leftChild, isNumber)
+  //    && isVariableSpeciesOrParameter(rightChild))
+  //  {
+  //   
+  //    if (isNumber)
+  //    {
+  //        value->k_value = "number";
+  //        value->k_real_value = leftChild->getValue();
+  //    }
+  //    else
+  //    {
+  //        value->k_value = leftChild->getName();
+  //    }
+  //    value->x_value = rightChild->getName();
+  //    value->dxdt_expression = getODEFor(rightChild->getName());
+  //    value->type = TYPE_K_MINUS_X;
+  //    value->current = node;
+  //    return true;
+  //  }
+  //  // left child + with it's left child const we have k+v-x
+  //  // left child + with it's left child k-x+w-y we have already finished
+  //  else if (leftChild->getType() == AST_PLUS)
+  //  {
+
+  //      // TO DO fix this for k Or w being a number
+  //      if (value->type == TYPE_K_MINUS_X_PLUS_W_MINUS_Y)
+  //      {
+  //          return true;
+  //      }
+  //      else if (analyseNode(leftChild, value))
+  //      {
+  //          if (value->type == TYPE_K_MINUS_X_PLUS_W_MINUS_Y)
+  //          {
+  //              value->y_value = rightChild->getName();
+  //              value->dydt_expression = getODEFor(rightChild->getName());
+  //              value->current = node;
+  //              return true;
+  //          }
+  //          else if (value->type == TYPE_K_MINUS_X)
+  //          {
+  //              value->y_value = rightChild->getName();
+  //              value->dydt_expression = getODEFor(rightChild->getName());
+  //              value->w_expression = leftChild->getChild(1);
+  //              value->type = TYPE_K_MINUS_X_PLUS_W_MINUS_Y;
+  //              value->current = node;
+  //              return true;
+
+  //          }
+  //      }
+  //      else if (isNumericalConstantOrConstantParameter(leftChild->getChild(0), isNumber))
+  //      {
+  //          value->k_value = leftChild->getChild(0)->getName();
+  //          value->x_value = rightChild->getName();
+  //          value->dxdt_expression = getODEFor(rightChild->getName());
+  //          value->v_expression = leftChild->getChild(1);
+  //          value->type = TYPE_K_PLUS_V_MINUS_X;
+  //          value->current = node;
+  //          return true;
+  //      }
+  //  }
+  //  else if (leftChild->getType() == AST_MINUS
+  //    && isVariableSpeciesOrParameter(leftChild->getRightChild()))
+  //  {
+  //    // if left child is k+v-x we have k+v-x-ymean you'll have to come againthat you can pay3
+  //    // or if left child is k-x we have k-x-y
+  //    if (analyseNode(leftChild, value))
+  //    {
+  //      if (value->type == TYPE_K_PLUS_V_MINUS_X)
+  //      {
+  //        value->type = TYPE_K_PLUS_V_MINUS_X_MINUS_Y;
+  //        value->y_value = rightChild->getName();
+  //        value->dydt_expression = getODEFor(value->y_value);
+  //        value->current = node;
+  //        return true;
+  //      }
+  //      else if (value->type == TYPE_K_MINUS_X)
+  //      {
+  //        value->y_value = rightChild->getName();
+  //        value->dydt_expression = getODEFor(rightChild->getName());
+  //        value->type = TYPE_K_MINUS_X_MINUS_Y;
+  //        value->current = node;
+  //        return true;
+  //      }
+  //    }
+  //    return false;
+  //  }
+  //  break;
+  //default:
+  //  return false;
+  //}
+  //return false;
+
+
 /*
 * Return the ODE for the given variable 
 * or an ASTNode representing zero if there is no time derivative
@@ -414,6 +746,69 @@ ExpressionAnalyser::getODEFor(std::string name)
   zero->setValue(0.0);
   return zero->deepCopy();
 }
+
+void
+ExpressionAnalyser::detect_minusXPlusYOnly()
+{
+    for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++)
+    {
+        std::pair<std::string, ASTNode*> ode = mODEs.at(odeIndex);
+        ASTNode* odeRHS = ode.second;
+        //odeRHS->decompose();
+        odeRHS->reduceToBinary();
+        odeRHS->decompose();
+        List* operators = odeRHS->getListOfNodes((ASTNodePredicate)ASTNode_isOperator);
+        ListIterator it = operators->begin();
+
+        while (it != operators->end())
+        {
+            ASTNode* currentNode = (ASTNode*)*it;
+                cout << "current node in not -x+y analyze: " << SBML_formulaToL3String(currentNode) << endl;
+            if (currentNode->getType() != AST_PLUS)
+            {
+                it++;
+                continue;
+            }
+            SubstitutionValues_t* value = createBlankSubstitutionValues();
+            unsigned int numChildren = currentNode->getNumChildren();
+            ASTNode* rightChild = currentNode->getRightChild();
+            ASTNode* leftChild = currentNode->getLeftChild();
+            //  -x+y node binary; plus; left child type minus; rightchild var/const
+            //           +
+            //        -     y
+            //        x
+            if (numChildren != 2 || rightChild->getType() != AST_NAME
+                || leftChild->getType() != AST_MINUS
+                || leftChild->getNumChildren() != 1)
+            {
+                it++;
+                continue;
+            }
+
+            // if we get to this point, the only thing left to check is 
+            // whether the ->left->right grandchild (the x in -x+y) is a variable species.
+            if (isVariableSpeciesOrParameter(leftChild->getChild(0)))
+            {
+                value->x_value = leftChild->getChild(0)->getName();
+                value->y_value = rightChild->getName();
+                value->dydt_expression = getODEFor(rightChild->getName());
+                value->dxdt_expression = getODEFor(leftChild->getChild(0)->getName());
+                value->type = TYPE_MINUS_X_PLUS_Y;
+                value->current = currentNode;
+                value->odeIndex = odeIndex;
+                if (!hasExpressionAlreadyBeenRecorded(value))
+                {
+                    printSubstitutionValues(value);
+                    mExpressions.push_back(value);
+                }
+
+            }
+
+            it++;
+        }
+    }
+
+}
 void
 ExpressionAnalyser::analyse(bool minusXPlusYOnly)
 {
@@ -421,6 +816,7 @@ ExpressionAnalyser::analyse(bool minusXPlusYOnly)
   {
     std::pair<std::string, ASTNode*> ode = mODEs.at(odeIndex);
     ASTNode* odeRHS = ode.second;
+    odeRHS->decompose();
     odeRHS->reduceToBinary();
     List* operators = odeRHS->getListOfNodes((ASTNodePredicate)ASTNode_isOperator);
     ListIterator it = operators->begin();
@@ -428,26 +824,24 @@ ExpressionAnalyser::analyse(bool minusXPlusYOnly)
     while (it != operators->end())
     {
       ASTNode* currentNode = (ASTNode*)*it;
-      //cout << "current node " << SBML_formulaToL3String(currentNode)  << endl;
+      if (minusXPlusYOnly)
+        cout << "current node in -x+y analyze: " << SBML_formulaToL3String(currentNode)  << endl;
+      else
+        cout << "current node in not -x+y analyze: " << SBML_formulaToL3String(currentNode) << endl;
       if (minusXPlusYOnly && currentNode->getType() != AST_PLUS)
       {
         it++;
         continue;
       }
-      SubstitutionValues_t* value = new SubstitutionValues_t;
-      value->type = TYPE_UNKNOWN;
-      value->k_real_value = util_NaN();
-      value->dxdt_expression = NULL;
-      value->dydt_expression = NULL;
-      value->v_expression = NULL;
-      value->w_expression = NULL;
+      SubstitutionValues_t* value = createBlankSubstitutionValues();
+
       if (analyseNode(currentNode, value))
       {
         value->odeIndex = odeIndex;
-        if (!hasExpressionAlreadyRecorded(value))
+        if (!hasExpressionAlreadyBeenRecorded(value))
         {
             //printSubstitutionValues(value);
-          mExpressions.push_back(value);
+            mExpressions.push_back(value);
         }
       }
       it++;
@@ -455,17 +849,35 @@ ExpressionAnalyser::analyse(bool minusXPlusYOnly)
   }
 }
 
+void ExpressionAnalyser::orderExpressions()
+{
+  for (unsigned int i = 0; i < mExpressions.size(); i++)
+  {
+    for (unsigned int j = i + 1; j < mExpressions.size(); j++)
+    {
+      if (mExpressions[i]->type > mExpressions[j]->type)
+      {
+        SubstitutionValues_t* temp = mExpressions[i];
+        mExpressions[i] = mExpressions[j];
+        mExpressions[j] = temp;
+      }
+    }
+  }
+}
+
 void
 ExpressionAnalyser::detectHiddenSpecies(List * hiddenSpecies)
 {
   // find -x+y and replace with y-x 
-  analyse(true);
+  // actually don't have to do this if we decompose the AST
+  // 
+  //detect_minusXPlusYOnly();
   //for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++)
   //{
-  //    cout << mODEs[odeIndex].first << ": " << SBML_formulaToL3String(mODEs[odeIndex].second) << endl;
+  //    cout << mODEs[odeIndex].first << ": " << m(mODEs[odeIndex].second) << endl;
   //}
 
-  reorderMinusXPlusYIteratively();
+  //reorderMinusXPlusYIteratively();
   mExpressions.clear();
   
   // find cases of k-x/k-x-y/k+v-x/k+v-x-y/k-x+w-y
@@ -474,31 +886,7 @@ ExpressionAnalyser::detectHiddenSpecies(List * hiddenSpecies)
   //{
   //    cout << mODEs[odeIndex].first << ": " << SBML_formulaToL3String(mODEs[odeIndex].second) << endl;
   //}
-
-
-  for (unsigned int i = 0; i < mExpressions.size(); i++)
-  {
-    SubstitutionValues_t *exp = mExpressions.at(i);
-    for (unsigned int j = 0; j < mODEs.size(); j++)
-    {
-        std::pair<std::string, ASTNode*> ode = mODEs.at(j);
-        ASTNode* odeRHS = ode.second;
-        bool index = isParameterAlreadyCreated(exp->z_value);
-        if (index == true)
-        {
-            replaceExpressionWithNewParameter(odeRHS, exp);
-        }
-        else
-        {
-            std::string zName = getUniqueNewParameterName();
-            exp->z_value = zName;
-            mNewVarCount++;
-            replaceExpressionWithNewParameter(odeRHS, exp);
-            addParametersAndRateRules(hiddenSpecies, exp);
-        }
-        //cout << "ode in main: " << SBML_formulaToL3String(odeRHS) << endl;
-    }
-  }
+  substituteParametersForExpressions(hiddenSpecies);
 }
 
 /*
@@ -656,12 +1044,13 @@ ExpressionAnalyser::replaceExpressionWithNewParameter(ASTNode* ode, Substitution
 {
   if (exp->type == TYPE_K_MINUS_X || exp->type == TYPE_K_MINUS_X_MINUS_Y)
   {
+      printSubstitutionValues(exp);
     replaceExpressionInNodeWithVar(ode, exp->current, exp->z_value);
     //cout << "ode in new param var: " << SBML_formulaToL3String(ode) << endl;
     for (unsigned int i = 0; i < mExpressions.size(); i++)
     {
       SubstitutionValues_t *thisexp = mExpressions.at(i);
-
+      printSubstitutionValues(thisexp);
       if (thisexp->dxdt_expression != NULL)
       {
         replaceExpressionInNodeWithVar(thisexp->dxdt_expression, exp->current, exp->z_value);
@@ -755,6 +1144,178 @@ ExpressionAnalyser::isParameterAlreadyCreated(std::string& name)
     }
 }
 
+bool ExpressionAnalyser::isTypeKminusXminusY(unsigned int numChildren, ASTNode* rightChild, 
+    ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value)
+{
+    // for this type of expression, k-x-y, the left child is the expression k-x, the right child is the y variable
+    // the type is MINUS and the number of children is 2
+    if (numChildren != 2 || type != AST_MINUS)
+    {
+        return false;
+    }
+
+
+    if (isTypeKminusX(leftChild->getNumChildren(), leftChild->getRightChild(), 
+        leftChild->getLeftChild(), leftChild->getType(), value) && 
+        isVariableSpeciesOrParameter(rightChild))
+    {
+        value->type = TYPE_K_MINUS_X_MINUS_Y;
+        value->y_value = rightChild->getName();
+        value->dydt_expression = getODEFor(rightChild->getName());
+        return true;
+    }
+    return false;
+}
+
+bool ExpressionAnalyser::isTypeKminusX(unsigned int numChildren, ASTNode* rightChild, 
+    ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value)
+{   
+    // for this type of expression, k-x, the left child is the parameter k, the right child is the x variable
+    // the type is MINUS and the number of children is 2
+    bool isNumber = false;
+    if (numChildren != 2 || type != AST_MINUS)
+    {
+        return false;
+    }
+
+
+    if (isNumericalConstantOrConstantParameter(leftChild, isNumber)
+      && isVariableSpeciesOrParameter(rightChild))
+    {   
+      if (isNumber)
+      {
+          value->k_value = "number";
+          value->k_real_value = leftChild->getValue();
+      }
+      else
+      {
+          value->k_value = leftChild->getName();
+      }
+      value->x_value = rightChild->getName();
+      value->dxdt_expression = getODEFor(rightChild->getName());
+      value->type = TYPE_K_MINUS_X;
+      return true;
+    }
+    return false;
+}
+
+bool ExpressionAnalyser::isTypeKplusVminusX(unsigned int numChildren, ASTNode* rightChild, 
+    ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value)
+{
+    // for this type of expression, k+v-x, the left child is the expression k+v, the right child is the x variable
+    // the type is MINUS and the number of children is 2
+    if (numChildren != 2 || type != AST_MINUS)
+    {
+        return false;
+    }
+
+
+    if (isTypeKplusV(leftChild->getNumChildren(), leftChild->getRightChild(),
+        leftChild->getLeftChild(), leftChild->getType(), value) &&
+        isVariableSpeciesOrParameter(rightChild))
+    {
+        value->type = TYPE_K_PLUS_V_MINUS_X;
+        value->x_value = rightChild->getName();
+        value->dxdt_expression = getODEFor(rightChild->getName());
+        return true;
+    }
+    return false;
+}
+
+bool ExpressionAnalyser::isTypeKplusV(unsigned int numChildren, ASTNode* rightChild, 
+    ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value)
+{
+    // for this type of expression, k+v, the left child is the constant k, the right child is the v variable
+    // the type is PLUS and the number of children is 2
+    if (numChildren != 2 || type != AST_PLUS)
+    {
+        return false;
+    }
+
+    bool isNumber = false;
+    if (isNumericalConstantOrConstantParameter(leftChild, isNumber))
+    {   
+      if (isNumber)
+      {
+          value->k_value = "number";
+          value->k_real_value = leftChild->getValue();
+      }
+      else
+      {
+          value->k_value = leftChild->getName();
+      }
+      value->v_expression = rightChild;
+      return true;
+    }
+    return false;
+}
+
+bool ExpressionAnalyser::isTypeKplusVminusXminusY(unsigned int numChildren, ASTNode* rightChild, ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value)
+{
+    // for this type of expression, k+v-x-y, the left child is the expression k+v-x, the right child is the y variable
+    // the type is MINUS and the number of children is 2
+    if (numChildren != 2 || type != AST_MINUS)
+    {
+        return false;
+    }
+
+
+    if (isTypeKplusVminusX(leftChild->getNumChildren(), leftChild->getRightChild(),
+        leftChild->getLeftChild(), leftChild->getType(), value) &&
+        isVariableSpeciesOrParameter(rightChild))
+    {
+        value->type = TYPE_K_PLUS_V_MINUS_X_MINUS_Y;
+        value->y_value = rightChild->getName();
+        value->dydt_expression = getODEFor(rightChild->getName());
+        return true;
+    }
+    return false;
+}
+
+bool ExpressionAnalyser::isTypeKminusXplusWminusY(unsigned int numChildren, ASTNode* rightChild, 
+    ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value)
+{
+    // for this type of expression, k-x+w-y, the left child is the expression w+k-x, the right child is the y variable
+    // the type is MINUS and the number of children is 2
+    if (numChildren != 2 || type != AST_MINUS)
+    {
+        return false;
+    }
+
+
+    if (isTypeWplusKminusX(leftChild->getNumChildren(), leftChild->getRightChild(),
+        leftChild->getLeftChild(), leftChild->getType(), value) &&
+        isVariableSpeciesOrParameter(rightChild))
+    {
+        value->y_value = rightChild->getName();
+        value->dydt_expression = getODEFor(rightChild->getName());
+        return true;
+    }
+    return false;
+}
+
+bool ExpressionAnalyser::isTypeWplusKminusX(unsigned int numChildren, ASTNode* rightChild, 
+    ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value)
+{
+    // for this type of expression, w+k-x, the left child is the expression w, the right child is the k-x expression
+    // the type is PLUS and the number of children is 2
+    if (numChildren != 2 || type != AST_PLUS)
+    {
+        return false;
+    }
+
+
+    if (isTypeKminusX(rightChild->getNumChildren(), rightChild->getRightChild(),
+        rightChild->getLeftChild(), rightChild->getType(), value))
+    {
+        value->type = TYPE_K_MINUS_X_PLUS_W_MINUS_Y;
+        value->w_expression = leftChild;
+        return true;
+    }
+    return false;
+}
+
+
 /*
  * Check whether for node is a name node representing species or a non constant parameter
 */
diff --git a/src/sbml/conversion/ExpressionAnalyser.h b/src/sbml/conversion/ExpressionAnalyser.h
index 6d356cd38d..b569ccb3a2 100644
--- a/src/sbml/conversion/ExpressionAnalyser.h
+++ b/src/sbml/conversion/ExpressionAnalyser.h
@@ -75,8 +75,8 @@ typedef enum
 } ExpressionType_t;
 
 /*
-*
-|*/
+* the structure to contain the substitution values
+*/
 struct SubstitutionValues_t {
   std::string k_value;
   double k_real_value;
@@ -89,6 +89,7 @@ struct SubstitutionValues_t {
   ExpressionType_t type;
   ASTNode* current;
   std::string z_value;
+  ASTNode* z_expression;
   unsigned int odeIndex;
 };
 
@@ -155,23 +156,60 @@ class LIBSBML_EXTERN ExpressionAnalyser
   int setModel(Model* m);
 
 
+  SubstitutionValues_t* createBlankSubstitutionValues();
+  unsigned int getNumExpressions();
+  SubstitutionValues_t* getExpression(unsigned int index);
+
+
   void detectHiddenSpecies(List * hiddenSpecies);
 
+  bool analyseNode(ASTNode* node, SubstitutionValues_t* value);
+
+  void analyse(bool minusXPlusYOnly = false);
+
+  void orderExpressions();
+
+  void printSubstitutionValues(const SubstitutionValues_t* values1);
+  bool areIdenticalSubstitutionValues(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+
+  bool expressionExists(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+
+  void substituteParametersForExpressions(List* hiddenSpecies);
+
+  /*
+  * identify instances of - x + y within formula and create expressions
+  * 
+  */
+  void detect_minusXPlusYOnly();
+
+  void reorderMinusXPlusYIteratively();
 
 private:
   /** @cond doxygenLibsbmlInternal */
+    bool matchesK(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+
+    bool matchesKValue(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    bool matchesKRealValue(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    bool matchesXValue(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    bool matchesYValue(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    bool matchesVExpression(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    bool matchesWExpression(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    bool matchesDxdtExpression(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    bool matchesDydtExpression(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    bool matchesCurrentNode(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+
+
+    void substituteParameters(List* hiddenSpecies, SubstitutionValues_t* values);
 
-    bool areIdenticalSubstitutionValues(SubstitutionValues_t* values1, SubstitutionValues_t* values2);
+    SubstitutionValues_t* getSubstitutionValuesByType(ExpressionType_t type, size_t index = 0);
 
-    void printSubstitutionValues(const SubstitutionValues_t* values1);
+
+    //void printSubstitutionValues(const SubstitutionValues_t* values1);
 
   // functions that represents steps of algo 3.1
 
-  void reorderMinusXPlusYIteratively();
 
 
-  void analyse(bool minusXPlusYOnly = false);
-  bool analyseNode(ASTNode* node, SubstitutionValues_t* value);
   
   /*
   * Return the ODE for the given variable
@@ -199,7 +237,7 @@ class LIBSBML_EXTERN ExpressionAnalyser
   /*
    * Loops through expressions already recorded and checks for exact matches
    */
-  bool hasExpressionAlreadyRecorded(SubstitutionValues_t* value);
+  bool hasExpressionAlreadyBeenRecorded(SubstitutionValues_t* value);
 
   /**
    * Searches for a node's parent and its index as the parent's child in a one-directional tree (nodes know their children, but not their parent).
@@ -235,12 +273,36 @@ class LIBSBML_EXTERN ExpressionAnalyser
   */
   bool isParameterAlreadyCreated(std::string& name);
 
+  bool isTypeKminusXminusY(unsigned int numChildren, ASTNode* rightChild, 
+      ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value);
+
+  bool isTypeKminusX(unsigned int numChildren, ASTNode* rightChild,
+      ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value);
+
+  bool isTypeKplusVminusX(unsigned int numChildren, ASTNode* rightChild,
+      ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value);
+ 
+  bool isTypeKplusV(unsigned int numChildren, ASTNode* rightChild,
+      ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value);
+
+  bool isTypeKplusVminusXminusY(unsigned int numChildren, ASTNode* rightChild,
+      ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value);
+
+  bool isTypeKminusXplusWminusY(unsigned int numChildren, ASTNode* rightChild,
+      ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value);
+
+  bool isTypeWplusKminusX(unsigned int numChildren, ASTNode* rightChild,
+      ASTNode* leftChild, ASTNodeType_t type, SubstitutionValues_t* value);
+
   // member variables populated during analysis
   pairODEs mODEs;
 
   Model* mModel;
 
   std::vector <SubstitutionValues_t*> mExpressions;
+   
+  // list of hidden species that are identified during the analysis
+  List* mHiddenSpecies;
 
   // variables to ensure unique new parameter name
 
diff --git a/src/sbml/conversion/SBMLRateRuleConverter.cpp b/src/sbml/conversion/SBMLRateRuleConverter.cpp
index 7f2744ed93..7d0a1b0792 100644
--- a/src/sbml/conversion/SBMLRateRuleConverter.cpp
+++ b/src/sbml/conversion/SBMLRateRuleConverter.cpp
@@ -241,6 +241,7 @@ SBMLRateRuleConverter::convert()
   }
 
   // Fages algo 3.6 Steps 1-2
+  populateInitialODEinfo();
   populateODEinfo();
 
   if (getMathNotSupportedFlag() == true)
@@ -835,37 +836,45 @@ SBMLRateRuleConverter::addToTerms(ASTNode* node, bool isToplevel)
 
 
 void 
-SBMLRateRuleConverter::populateODEinfo()
+SBMLRateRuleConverter::populateInitialODEinfo()
 {
-  Model* model = mDocument->getModel();
+    Model* model = mDocument->getModel();
 
-  // Fages algo 3.6 create set O
-  // create pairs of variables and their corresponding ODE
-  // eg. ODEs[0] = [S1, -k1*S1]
-  //     ODES[1] = [S2, k1*S1]
-  //     ODES[2] = [S3, k2*S3]
-  //
-  for (unsigned int n = 0; n < model->getNumSpecies(); n++)
-  {
-    Species *s = model->getSpecies(n);
-    if (s->getConstant() == false)
+    // Fages algo 3.6 create set O
+    // create pairs of variables and their corresponding ODE
+    // eg. ODEs[0] = [S1, -k1*S1]
+    //     ODES[1] = [S2, k1*S1]
+    //     ODES[2] = [S3, k2*S3]
+    //
+    for (unsigned int n = 0; n < model->getNumSpecies(); n++)
     {
-      addODEPair(s->getId(), model);
+        Species* s = model->getSpecies(n);
+        if (s->getConstant() == false)
+        {
+            addODEPair(s->getId(), model);
+        }
     }
-  }
-  for (unsigned int n = 0; n < model->getNumParameters(); n++)
-  {
-    Parameter *p = model->getParameter(n);
-    if (p->getConstant() == false)
+    for (unsigned int n = 0; n < model->getNumParameters(); n++)
     {
-      addODEPair(p->getId(), model);
+        Parameter* p = model->getParameter(n);
+        if (p->getConstant() == false)
+        {
+            addODEPair(p->getId(), model);
+        }
     }
-  }
 
-  for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++)
-  {
-    cout << mODEs[odeIndex].first << ": " << SBML_formulaToL3String(mODEs[odeIndex].second) << endl;
-  }
+    for (unsigned int odeIndex = 0; odeIndex < mODEs.size(); odeIndex++)
+    {
+        cout << mODEs[odeIndex].first << ": " << SBML_formulaToL3String(mODEs[odeIndex].second) << endl;
+    }
+}
+
+
+void
+SBMLRateRuleConverter::populateODEinfo()
+{
+    Model* model = mDocument->getModel();
+
   // implement Algo 3.1 here (hidden variables!)
   // check for hidden variables, and add an appropriate ODE if a hidden variable is found
   ExpressionAnalyser *ea = new ExpressionAnalyser(model, mODEs);
diff --git a/src/sbml/conversion/SBMLRateRuleConverter.h b/src/sbml/conversion/SBMLRateRuleConverter.h
index 67a0b9a77c..9805aab548 100644
--- a/src/sbml/conversion/SBMLRateRuleConverter.h
+++ b/src/sbml/conversion/SBMLRateRuleConverter.h
@@ -188,7 +188,7 @@ class LIBSBML_EXTERN SBMLRateRuleConverter : public SBMLConverter
   typedef std::vector< std::vector<double> > setRnCoeffs;
   typedef std::vector< std::pair< std::string, ASTNode*> >::iterator odeIt;
 
-  pairODEs getOde() { return mODEs; };
+  pairODEs getOdePairs() { return mODEs; };
 
   setCoeff getCoeff() { return mCoefficients; };
   setRnCoeffs getReactants() { return mReactants; };
@@ -200,6 +200,9 @@ class LIBSBML_EXTERN SBMLRateRuleConverter : public SBMLConverter
   std::vector< std::vector<bool> > getNegDer() { return mNegDerivative; };
   bool getMathNotSupportedFlag() const;
 
+  void populateInitialODEinfo();
+
+  void populateODEinfo();
 
 private:
   /** @cond doxygenLibsbmlInternal */
@@ -208,7 +211,6 @@ class LIBSBML_EXTERN SBMLRateRuleConverter : public SBMLConverter
 
   bool isDocumentAppropriate(OperationReturnValues_t& returnVal);
 
-  void populateODEinfo();
 
   void populateReactionCoefficients();
 
diff --git a/src/sbml/conversion/test/TestExpressionAnalyser.cpp b/src/sbml/conversion/test/TestExpressionAnalyser.cpp
new file mode 100644
index 0000000000..df2df3ee56
--- /dev/null
+++ b/src/sbml/conversion/test/TestExpressionAnalyser.cpp
@@ -0,0 +1,555 @@
+/**
+ * @file    TestExpressionAnalyser.cpp
+ * @brief   Tests for raterule to reaction converter
+ * @author  Sarah Keating
+ * @author  Alessandro Felder
+ *
+ * <!--------------------------------------------------------------------------
+ * This file is part of libSBML.  Please visit http://sbml.org for more
+ * information about SBML, and the latest version of libSBML.
+ *
+ * Copyright (C) 2013-2018 jointly by the following organizations:
+ *     1. California Institute of Technology, Pasadena, CA, USA
+ *     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
+ *     3. University of Heidelberg, Heidelberg, Germany
+ *
+ * Copyright (C) 2009-2013 jointly by the following organizations: 
+ *     1. California Institute of Technology, Pasadena, CA, USA
+ *     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
+ *  
+ * Copyright (C) 2006-2008 by the California Institute of Technology,
+ *     Pasadena, CA, USA 
+ *  
+ * Copyright (C) 2002-2005 jointly by the following organizations: 
+ *     1. California Institute of Technology, Pasadena, CA, USA
+ *     2. Japan Science and Technology Agency, Japan
+ * 
+ * This library is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation.  A copy of the license agreement is provided
+ * in the file named "LICENSE.txt" included with this software distribution
+ * and also available online as http://sbml.org/software/libsbml/license.html
+ * ---------------------------------------------------------------------- -->*/
+
+#include <sbml/common/common.h>
+
+#include <sbml/SBase.h>
+#include <sbml/SBMLTypes.h>
+
+#include <sbml/conversion/SBMLConverter.h>
+#include <sbml/conversion/SBMLConverterRegistry.h>
+#include <sbml/conversion/SBMLRateRuleConverter.h>
+#include <sbml/conversion/ExpressionAnalyser.h>
+
+#include <sbml/math/FormulaParser.h>
+
+#include <string>
+#include <iostream>
+#include <check.h>
+
+using namespace std;
+LIBSBML_CPP_NAMESPACE_USE
+
+CK_CPPSTART
+
+static Model* m;
+static SBMLDocument* d;
+static ConversionProperties props;
+static SBMLRateRuleConverter* converter;  
+
+
+static bool
+equals(const char* expected, const char* actual)
+{
+  if (!strcmp(expected, actual)) return true;
+
+  printf("\nStrings are not equal:\n");
+  printf("Expected:\n[%s]\n", expected);
+  printf("Actual:\n[%s]\n", actual);
+
+  return false;
+}
+
+static bool
+formulas_equal(const char* expected, ASTNode* actual)
+{
+	return equals(expected, SBML_formulaToL3String(actual));
+}
+
+extern char *TestDataDirectory;
+
+static Parameter* setupZeroParameter(Model* model, const char* name, bool is_constant)
+{
+	Parameter* parameter = model->createParameter();
+	parameter->setId(name);
+	parameter->setConstant(is_constant);
+	parameter->setValue(0);
+	return parameter;
+}
+
+Species* setupSpecies(Model* model, const char* name, const char* compartment) {
+    Species* species = model->createSpecies();
+    species->setId(name);
+    species->setCompartment(compartment);
+    species->setInitialAmount(0);
+	species->setConstant(false);
+    return species;
+}
+
+Model* setupModel(SBMLDocument* doc) {
+    Model* model = doc->createModel();
+    model->setId("m");
+
+	// create compartment
+    Compartment* compartment = model->createCompartment();
+    compartment->setId("c");
+    compartment->setSpatialDimensions(3.0);
+    compartment->setSize(1);
+    compartment->setConstant(true);
+
+	// create species
+    setupSpecies(model, "x", "c");
+    setupSpecies(model, "y", "c");
+    setupSpecies(model, "a", "c");
+	setupSpecies(model, "b", "c");
+
+	// create parameters
+    setupZeroParameter(model, "k", true);
+    setupZeroParameter(model, "v", false);
+	setupZeroParameter(model, "w", false);
+
+
+    return model;
+}
+
+
+void
+ExpressionAnalyser_setup(void)
+{
+	props.addOption("inferReactions", true);
+
+	converter = new SBMLRateRuleConverter();
+	converter->setProperties(&props);
+
+	d = new SBMLDocument();
+	m = setupModel(d);
+	converter->setDocument(d);
+}
+
+void
+ExpressionAnalyser_teardown(void)
+{
+	delete converter;
+	delete d;
+}
+
+extern char *TestDataDirectory;
+START_TEST(test_analyse)
+{
+	RateRule* rr = d->getModel()->createRateRule();
+	rr->setVariable("a");
+	rr->setMath(SBML_parseFormula("k-x-y"));
+	converter->populateInitialODEinfo();
+	ExpressionAnalyser* analyser = new ExpressionAnalyser(m, converter->getOdePairs());
+
+	fail_unless(analyser->getNumExpressions() == 0);
+    
+	analyser->analyse();
+
+	fail_unless(analyser->getNumExpressions() == 1);
+	SubstitutionValues_t *value = analyser->getExpression(0);
+	fail_unless(value->k_value == "k");
+	fail_unless(value->x_value == "x");
+	fail_unless(value->y_value == "y");
+	fail_unless(value->z_value.empty());
+	fail_unless(value->type == TYPE_K_MINUS_X_MINUS_Y);
+	fail_unless(formulas_equal("k - x - y", value->current));
+	fail_unless(formulas_equal( "0", value->dxdt_expression));
+	fail_unless(formulas_equal( "0", value->dydt_expression));
+	fail_unless(value->v_expression == NULL);
+	fail_unless(value->w_expression == NULL);
+	fail_unless(value->z_expression == NULL);
+	fail_unless(value->odeIndex == 0);
+	fail_unless(util_isNaN(value->k_real_value));
+}
+END_TEST
+
+START_TEST(test_analyse_1)
+{
+	RateRule* rr = d->getModel()->createRateRule();
+	rr->setVariable("a");
+	rr->setMath(SBML_parseFormula("k + v - x - y"));
+	converter->populateInitialODEinfo();
+	ExpressionAnalyser* analyser = new ExpressionAnalyser(m, converter->getOdePairs());
+
+	fail_unless(analyser->getNumExpressions() == 0);
+
+	analyser->analyse();
+
+	fail_unless(analyser->getNumExpressions() == 1);
+	SubstitutionValues_t* value = analyser->getExpression(0);
+	fail_unless(value->k_value == "k");
+	fail_unless(value->x_value == "x");
+	fail_unless(value->y_value == "y");
+	fail_unless(value->z_value.empty());
+	fail_unless(value->type == TYPE_K_PLUS_V_MINUS_X_MINUS_Y);
+	fail_unless(formulas_equal("k + v - x - y", value->current));
+	fail_unless(formulas_equal("0", value->dxdt_expression));
+	fail_unless(formulas_equal("0", value->dydt_expression));
+	fail_unless(formulas_equal("v", value->v_expression));
+	fail_unless(value->w_expression == NULL);
+	fail_unless(value->z_expression == NULL);
+	fail_unless(value->odeIndex == 0);
+	fail_unless(util_isNaN(value->k_real_value));
+}
+END_TEST
+START_TEST(test_analyse_2)
+{
+	RateRule* rr = d->getModel()->createRateRule();
+	rr->setVariable("a");
+	rr->setMath(SBML_parseFormula("k - x + w - y"));
+	converter->populateInitialODEinfo();
+	ExpressionAnalyser* analyser = new ExpressionAnalyser(m, converter->getOdePairs());
+
+	fail_unless(analyser->getNumExpressions() == 0);
+
+	analyser->analyse();
+
+	fail_unless(analyser->getNumExpressions() == 1);
+	SubstitutionValues_t* value = analyser->getExpression(0);
+	fail_unless(value->k_value == "k");
+	fail_unless(value->x_value == "x");
+	fail_unless(value->y_value == "y");
+	fail_unless(value->z_value.empty());
+	fail_unless(value->type == TYPE_K_MINUS_X_PLUS_W_MINUS_Y);
+	fail_unless(formulas_equal("w + (k - x) - y", value->current));
+	fail_unless(formulas_equal("0", value->dxdt_expression));
+	fail_unless(formulas_equal("0", value->dydt_expression));
+	fail_unless(formulas_equal("w", value->w_expression));
+	fail_unless(value->v_expression == NULL);
+	fail_unless(value->z_expression == NULL);
+	fail_unless(value->odeIndex == 0);
+	fail_unless(util_isNaN(value->k_real_value));
+}
+END_TEST
+
+START_TEST(test_analyse_3)
+{
+	RateRule* rr = d->getModel()->createRateRule();
+	rr->setVariable("a");
+	rr->setMath(SBML_parseFormula("k - x"));
+	converter->populateInitialODEinfo();
+	ExpressionAnalyser* analyser = new ExpressionAnalyser(m, converter->getOdePairs());
+
+	fail_unless(analyser->getNumExpressions() == 0);
+
+	analyser->analyse();
+
+	fail_unless(analyser->getNumExpressions() == 1);
+	SubstitutionValues_t* value = analyser->getExpression(0);
+	fail_unless(value->k_value == "k");
+	fail_unless(value->x_value == "x");
+	fail_unless(value->y_value.empty());
+	fail_unless(value->z_value.empty());
+	fail_unless(value->type == TYPE_K_MINUS_X);
+	fail_unless(formulas_equal("k - x", value->current));
+	fail_unless(formulas_equal("0", value->dxdt_expression));
+	fail_unless(value->dydt_expression == NULL);
+	fail_unless(value->v_expression == NULL);
+	fail_unless(value->w_expression == NULL);
+	fail_unless(value->z_expression == NULL);
+	fail_unless(value->odeIndex == 0);
+	fail_unless(util_isNaN(value->k_real_value));
+}
+END_TEST
+
+START_TEST(test_analyse_4)
+{
+	RateRule* rr = d->getModel()->createRateRule();
+	rr->setVariable("a");
+	rr->setMath(SBML_parseFormula("k + v - x"));
+	converter->populateInitialODEinfo();
+	ExpressionAnalyser* analyser = new ExpressionAnalyser(m, converter->getOdePairs());
+
+	fail_unless(analyser->getNumExpressions() == 0);
+
+	analyser->analyse();
+
+	fail_unless(analyser->getNumExpressions() == 1);
+	SubstitutionValues_t* value = analyser->getExpression(0);
+	fail_unless(value->k_value == "k");
+	fail_unless(value->x_value == "x");
+	fail_unless(value->y_value.empty());
+	fail_unless(value->z_value.empty());
+	fail_unless(value->type == TYPE_K_PLUS_V_MINUS_X);
+	fail_unless(formulas_equal("k + v - x", value->current));
+	fail_unless(formulas_equal("0", value->dxdt_expression));
+	fail_unless(value->dydt_expression == NULL);
+	fail_unless(formulas_equal("v", value->v_expression));
+	fail_unless(value->w_expression == NULL);
+	fail_unless(value->z_expression == NULL);
+	fail_unless(value->odeIndex == 0);
+	fail_unless(util_isNaN(value->k_real_value));
+}
+END_TEST
+
+START_TEST(test_reorder_minusXplusYIteratively_simple)
+{
+	RateRule* rr = d->getModel()->createRateRule();
+	rr->setVariable("a");
+	rr->setMath(SBML_parseFormula("-x + y"));
+	converter->populateInitialODEinfo();
+	ExpressionAnalyser* analyser = new ExpressionAnalyser(m, converter->getOdePairs());
+
+	fail_unless(analyser->getNumExpressions() == 0);
+
+	analyser->detect_minusXPlusYOnly();
+
+	fail_unless(analyser->getNumExpressions() == 1);
+
+	analyser->reorderMinusXPlusYIteratively();
+	SubstitutionValues_t* value = analyser->getExpression(0);
+	fail_unless(value->k_value.empty());
+	fail_unless(value->x_value == "x");
+	fail_unless(value->y_value == "y");
+	fail_unless(value->z_value.empty());
+	fail_unless(value->type == TYPE_MINUS_X_PLUS_Y);
+	fail_unless(formulas_equal("-x + y", value->current));
+	fail_unless(formulas_equal("0", value->dxdt_expression));
+	fail_unless(formulas_equal("0", value->dydt_expression));
+	fail_unless(value->v_expression == NULL);
+	fail_unless(value->w_expression == NULL);
+	fail_unless(value->z_expression == NULL);
+	fail_unless(value->odeIndex == 0);
+	fail_unless(util_isNaN(value->k_real_value));
+}
+END_TEST
+
+
+
+START_TEST(test_order_of_replacements)
+{
+	ConversionProperties props;
+	props.addOption("inferReactions", true);
+
+	SBMLRateRuleConverter* converter = new SBMLRateRuleConverter();
+	converter->setProperties(&props);
+
+	std::string filename(TestDataDirectory);
+	filename += "mraterules7.xml";
+
+
+	SBMLDocument* d = readSBMLFromFile(filename.c_str());
+	Model* model = d->getModel();
+	fail_unless(model != NULL);
+	fail_unless(model->getNumParameters() == 2);
+
+	converter->setDocument(d);
+	converter->populateInitialODEinfo();
+	converter->populateODEinfo();
+	fail_unless(model->getNumParameters() == 3);
+
+	delete converter;
+	delete d;
+}
+END_TEST
+
+START_TEST(test_order_of_replacements1)
+{
+	ConversionProperties props;
+	props.addOption("inferReactions", true);
+
+	SBMLRateRuleConverter* converter = new SBMLRateRuleConverter();
+	converter->setProperties(&props);
+
+	std::string filename(TestDataDirectory);
+	filename += "mraterules7.xml";
+
+
+	SBMLDocument* d = readSBMLFromFile(filename.c_str());
+	Model* model = d->getModel();
+	fail_unless(model != NULL);
+	fail_unless(model->getNumParameters() == 2);
+
+	converter->setDocument(d);
+	converter->populateInitialODEinfo();
+	ASTNode* currentNode = const_cast<ASTNode*>(model->getRule(0)->getMath())  ;
+
+	ExpressionAnalyser* ea = new ExpressionAnalyser(model, converter->getOdePairs());
+
+	SubstitutionValues_t* value = new SubstitutionValues_t;
+	value->type = TYPE_UNKNOWN;
+	value->k_real_value = util_NaN();
+	value->dxdt_expression = NULL;
+	value->dydt_expression = NULL;
+	value->v_expression = NULL;
+	value->w_expression = NULL;
+
+	fail_unless(ea->analyseNode(currentNode, value));
+
+	fail_unless(value->type == TYPE_K_PLUS_V_MINUS_X_MINUS_Y);
+	fail_unless(util_isNaN(value->k_real_value));
+	fail_unless(value->k_value == "k");
+	fail_unless(value->x_value == "a");
+    fail_unless(value->y_value == "b");
+
+	fail_unless(formulas_equal("0", value->dxdt_expression));
+	fail_unless(formulas_equal("0", value->dydt_expression));
+	fail_unless(formulas_equal("v", value->v_expression));
+    fail_unless(value->w_expression == NULL);
+
+	fail_unless(formulas_equal("k + v - a - b", value->current));
+	fail_unless(value->z_value == "");
+
+	ea->analyse();
+
+	fail_unless(ea->getNumExpressions() == 1);
+
+	SubstitutionValues_t* value1 = new SubstitutionValues_t;
+	value1 = ea->getExpression(0);
+	fail_unless(ea->areIdenticalSubstitutionValues(value, value1));
+	fail_unless(value1->odeIndex == 0);
+
+	delete converter;
+	delete d;
+}
+END_TEST
+
+
+START_TEST(test_order_of_replacements2)
+{
+	ConversionProperties props;
+	props.addOption("inferReactions", true);
+
+	SBMLRateRuleConverter* converter = new SBMLRateRuleConverter();
+	converter->setProperties(&props);
+
+	std::string filename(TestDataDirectory);
+	filename += "mraterules5.xml";
+
+
+	SBMLDocument* d = readSBMLFromFile(filename.c_str());
+	Model* model = d->getModel();
+	fail_unless(model != NULL);
+
+	converter->setDocument(d);
+	converter->populateInitialODEinfo();
+	ASTNode* currentNode = const_cast<ASTNode*>(model->getRule(1)->getMath()->getChild(1));
+	cout << SBML_formulaToL3String(currentNode);
+	ExpressionAnalyser* ea = new ExpressionAnalyser(model, converter->getOdePairs());
+	ea->analyse();
+	fail_unless(ea->getNumExpressions() == 2);
+
+	SubstitutionValues_t* value1 = new SubstitutionValues_t;
+	value1 = ea->getExpression(0);
+	fail_unless(value1->type == TYPE_K_PLUS_V_MINUS_X_MINUS_Y);
+	fail_unless(util_isNaN(value1->k_real_value));
+	fail_unless(value1->k_value == "k");
+	fail_unless(value1->x_value == "a");
+	fail_unless(value1->y_value == "b");
+
+	fail_unless(formulas_equal("-1 * (k - a - b)", value1->dxdt_expression));
+	fail_unless(formulas_equal("c", value1->dydt_expression));
+	fail_unless(formulas_equal("v", value1->v_expression ));
+	fail_unless(value1->w_expression == NULL);
+
+	fail_unless(formulas_equal("k + v - a - b", value1->current));
+	fail_unless(value1->z_value == "");
+	fail_unless(value1->odeIndex == 0);
+
+	SubstitutionValues_t* value = new SubstitutionValues_t;
+	value = ea->getExpression(1);
+
+	fail_unless(value->type == TYPE_K_MINUS_X_MINUS_Y);
+	fail_unless(util_isNaN(value->k_real_value));
+	fail_unless(value->k_value == "k");
+	fail_unless(value->x_value == "a");
+	fail_unless(value->y_value == "b");
+
+	fail_unless(formulas_equal("-1 * (k - a - b)", value->dxdt_expression));
+	fail_unless(formulas_equal("c", value->dydt_expression));
+	fail_unless(value->v_expression == NULL);
+	fail_unless(value->w_expression == NULL);
+
+	fail_unless(formulas_equal("k - a - b", value->current));
+	fail_unless(value->z_value == "");
+
+	fail_unless(value->odeIndex == 1);
+
+	ea->orderExpressions();
+	SubstitutionValues_t* value2 = new SubstitutionValues_t;
+	value2 = ea->getExpression(0);
+	fail_unless(ea->areIdenticalSubstitutionValues(value, value2));
+
+	delete converter;
+	delete d;
+}
+END_TEST
+
+START_TEST(test_variations)
+{
+	ConversionProperties props;
+	props.addOption("inferReactions", true);
+
+	SBMLRateRuleConverter* converter = new SBMLRateRuleConverter();
+	converter->setProperties(&props);
+
+	SBMLDocument* d = new SBMLDocument(3, 1);
+	Model* model = setupModel(d);
+
+	RateRule* rr = model->createRateRule();
+    rr->setVariable("a");
+    rr->setMath(SBML_parseFormula("k-x-y"));
+	RateRule* rr1 = model->createRateRule();
+	rr1->setVariable("b");
+	rr1->setMath(SBML_parseFormula("-1*(k-x-y)"));
+	cout << d->toSBML() << endl;
+
+
+	converter->setDocument(d);
+	//converter->populateInitialODEinfo();
+	//ExpressionAnalyser* ea = new ExpressionAnalyser(model, converter->getOde());
+	//ea->analyse();
+
+
+	fail_unless(converter->convert() == LIBSBML_OPERATION_SUCCESS);
+
+	cout << d->toSBML() << endl;
+    delete converter;
+    delete d;
+
+}
+END_TEST
+
+Suite *
+create_suite_TestExpressionAnalyser (void)
+{ 
+	bool testing = false;
+Suite *suite = suite_create("ExpressionAnalyser");
+  TCase *tcase = tcase_create("ExpressionAnalyser");
+  tcase_add_checked_fixture(tcase,
+	  ExpressionAnalyser_setup, ExpressionAnalyser_teardown);
+
+  if (testing)
+  {
+	  tcase_add_test(tcase, test_analyse_2);
+  }
+  else
+  {
+	  tcase_add_test(tcase, test_analyse); //k-x-y
+	  tcase_add_test(tcase, test_analyse_1); //k+v-x-y
+	  tcase_add_test(tcase, test_analyse_2); //k-x+w-y
+	  tcase_add_test(tcase, test_analyse_3); //k-x
+	  tcase_add_test(tcase, test_analyse_4); //k+v-x
+	  //tcase_add_test(tcase, test_order_of_replacements);
+	  //tcase_add_test(tcase, test_order_of_replacements1);
+	  //tcase_add_test(tcase, test_order_of_replacements2);
+
+  }
+  suite_add_tcase(suite, tcase);
+
+  return suite;
+
+}
+END_C_DECLS
+
diff --git a/src/sbml/conversion/test/TestSBMLRateRuleConverter.cpp b/src/sbml/conversion/test/TestSBMLRateRuleConverter.cpp
index b754c1961a..28da3512bc 100644
--- a/src/sbml/conversion/test/TestSBMLRateRuleConverter.cpp
+++ b/src/sbml/conversion/test/TestSBMLRateRuleConverter.cpp
@@ -39,6 +39,7 @@
 #include <sbml/conversion/SBMLConverter.h>
 #include <sbml/conversion/SBMLConverterRegistry.h>
 #include <sbml/conversion/SBMLRateRuleConverter.h>
+#include <sbml/conversion/ExpressionAnalyser.h>
 
 #include <sbml/math/FormulaParser.h>
 
@@ -62,12 +63,18 @@ equals(const char* expected, const char* actual)
   return false;
 }
 
+static bool
+formulas_equal(const char* expected, ASTNode* actual)
+{
+	return equals(expected, SBML_formulaToL3String(actual));
+}
+
 extern char *TestDataDirectory;
 
 
 
 // helper function to set up a parameter with 0 value
-Parameter* setupZeroParameter(Model* model, const char* name, bool is_constant)
+static Parameter* setupZeroParameter(Model* model, const char* name, bool is_constant)
 {
 	Parameter* parameter = model->createParameter();
 	parameter->setId(name);
@@ -77,7 +84,6 @@ Parameter* setupZeroParameter(Model* model, const char* name, bool is_constant)
 }
 
 
-
 extern char *TestDataDirectory;
 
 START_TEST(test_conversion_raterule_converter_invalid)
@@ -1029,35 +1035,43 @@ START_TEST(test_model_valid_55)
 END_TEST
 
 
+
+
+
 Suite *
 create_suite_TestSBMLRateRuleConverter (void)
 { 
-  Suite *suite = suite_create("SBMLRateRuleConverter");
+	bool testing = false;
+Suite *suite = suite_create("SBMLRateRuleConverter");
   TCase *tcase = tcase_create("SBMLRateRuleConverter");
-
-    tcase_add_test(tcase, test_conversion_raterule_converter_invalid);
-     tcase_add_test(tcase, test_conversion_raterule_converter);
-  tcase_add_test(tcase, test_conversion_raterule_converter_non_standard_stoichiometry);
-  tcase_add_test(tcase, test_conversion_raterule_converter_hidden_variable);
-   tcase_add_test(tcase, test_crash_converter);
-   tcase_add_test(tcase, test_model);
-    tcase_add_test(tcase, test_model1);
-   tcase_add_test(tcase, test_model2);
-   tcase_add_test(tcase, test_model3);
-   tcase_add_test(tcase, test_model4);
-   tcase_add_test(tcase, test_model5);
-   tcase_add_test(tcase, test_model6);
-   tcase_add_test(tcase, test_model_valid_01); // this one works
-  tcase_add_test(tcase, test_model_valid_02);	// this one works
- tcase_add_test(tcase, test_model_valid_03);	// this one replaces the kinetic law'b' with the assignment rule for b
-  tcase_add_test(tcase, test_model_valid_04);	// this one works
-   tcase_add_test(tcase, test_model_valid_05);	// this one works
-   tcase_add_test(tcase, test_model_valid_51);	// this one works
- tcase_add_test(tcase, test_model_valid_52);	// this one works although puts the kinetic law in a different order
-  tcase_add_test(tcase, test_model_valid_53);	// this one has hidden species that I'm not picked up
-tcase_add_test(tcase, test_model_valid_54);	//
-  tcase_add_test(tcase, test_model_valid_55);	// 
-  
+  if (testing)
+  {
+	  tcase_add_test(tcase, test_conversion_raterule_converter_invalid);
+  }
+  else
+  {
+	  tcase_add_test(tcase, test_conversion_raterule_converter_invalid); 
+	  tcase_add_test(tcase, test_conversion_raterule_converter); 
+	  tcase_add_test(tcase, test_conversion_raterule_converter_non_standard_stoichiometry); 
+	  tcase_add_test(tcase, test_crash_converter); 
+	  tcase_add_test(tcase, test_conversion_raterule_converter_hidden_variable);
+	  tcase_add_test(tcase, test_model); 
+	  tcase_add_test(tcase, test_model1); // ??? not sure these are accurate
+	  tcase_add_test(tcase, test_model2); 
+	  tcase_add_test(tcase, test_model3);
+	  tcase_add_test(tcase, test_model4); 
+	  tcase_add_test(tcase, test_model5); // not working again
+	  //tcase_add_test(tcase, test_model6); // not working
+	  tcase_add_test(tcase, test_model_valid_01); // fixed
+	  tcase_add_test(tcase, test_model_valid_02);	
+	  tcase_add_test(tcase, test_model_valid_04);	
+	  tcase_add_test(tcase, test_model_valid_05);	// have changed the output model but I think it was only valid if compartment volume is one
+	  tcase_add_test(tcase, test_model_valid_51);	
+	  tcase_add_test(tcase, test_model_valid_52);	 
+	  tcase_add_test(tcase, test_model_valid_53); //	
+	  tcase_add_test(tcase, test_model_valid_54); //	
+	  tcase_add_test(tcase, test_model_valid_55); //	 
+  }
   suite_add_tcase(suite, tcase);
 
   return suite;
diff --git a/src/sbml/conversion/test/test-data/mraterules5.xml b/src/sbml/conversion/test/test-data/mraterules5.xml
index 9ba3424f1b..19507152f6 100644
--- a/src/sbml/conversion/test/test-data/mraterules5.xml
+++ b/src/sbml/conversion/test/test-data/mraterules5.xml
@@ -6,7 +6,7 @@
         </listOfCompartments>
         <listOfSpecies>
             <species id="c" name="c" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false"/>
-            <species id="a" name="c" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false"/>
+            <species id="a" name="a" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false"/>
             <species id="b" name="b" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false"/>
         </listOfSpecies>
         <listOfParameters>
@@ -52,6 +52,11 @@
                     </apply>
                 </math>
             </rateRule>
+            <rateRule variable="b">
+                <math xmlns="http://www.w3.org/1998/Math/MathML">
+                    <ci> c </ci>
+                </math>
+            </rateRule>
         </listOfRules>
     </model>
 </sbml>
diff --git a/src/sbml/conversion/test/test-data/mraterules7.xml b/src/sbml/conversion/test/test-data/mraterules7.xml
new file mode 100644
index 0000000000..57e20dfd1a
--- /dev/null
+++ b/src/sbml/conversion/test/test-data/mraterules7.xml
@@ -0,0 +1,36 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<sbml xmlns="http://www.sbml.org/sbml/level2/version4" level="2" version="4">
+    <model metaid="COPASI0" id="New_Model" name="New Model">
+        <listOfCompartments>
+            <compartment metaid="COPASI1" id="compartment" name="compartment" spatialDimensions="3" size="1" constant="true"/>
+        </listOfCompartments>
+        <listOfSpecies>
+            <species id="c" name="c" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false"/>
+            <species id="a" name="c" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false"/>
+            <species id="b" name="b" compartment="compartment" initialConcentration="1" boundaryCondition="false" constant="false"/>
+        </listOfSpecies>
+        <listOfParameters>
+            <parameter id="k" value="2.1"/>
+            <parameter id="v" value="2.1" constant="false"/>
+        </listOfParameters>
+        <listOfRules>
+            <rateRule variable="c">
+                <math xmlns="http://www.w3.org/1998/Math/MathML">
+                     <apply>
+                    <minus/>
+                    <apply>
+                        <minus/>
+                        <apply>
+                            <plus/>
+                            <ci> k </ci>
+                            <ci> v </ci>
+                        </apply>
+                        <ci> a </ci>
+                    </apply>
+                    <ci> b </ci>
+                    </apply>
+                </math>
+            </rateRule>
+        </listOfRules>
+    </model>
+</sbml>
diff --git a/src/sbml/conversion/test/test-data/test b/src/sbml/conversion/test/test-data/test
new file mode 100644
index 0000000000..83e56beae9
--- /dev/null
+++ b/src/sbml/conversion/test/test-data/test
@@ -0,0 +1,122 @@
+
+Model 53
+// Reactions: mine
+  J1: newVar1 => SimData_1 + newVar1; newVar1*I*k_1/(k_7 + newVar1);
+  J2: SimData_1 => newVar1; SimData_1*k_2/(SimData_1 + k_8);
+  J3: SimData_1 + newVar1 => SimData_1 + SimData_2 + newVar1; newVar1*SimData_1*k_3/(k_9 + newVar1);
+  J4: SimData_2 => ; SimData_2*k_4/(SimData_2 + k_10);
+  J5: SimData_1 + newVar1 => SimData_1 + SimData_3 + newVar1; newVar1*SimData_1*k_5/(k_11 + newVar1);
+  J6: SimData_2 + SimData_3 => SimData_2; SimData_3*SimData_2*k_6/(SimData_3 + k_12);
+
+// once from moccasin
+  J1: SimData_1             => SimData_6;             SimData_1*(k_2/(SimData_1 + k_8));
+  J2: SimData_6             => SimData_1;             k_1*SimData_6*(I/(SimData_6 + k_7));
+  J3: SimData_2             => SimData_5;             SimData_2*(k_4/(SimData_2 + k_10));
+  J4: SimData_1 + SimData_5 => SimData_1 + SimData_2; SimData_1*SimData_5*(k_3/(SimData_5 + k_9));
+  J5: SimData_2 + SimData_3 => SimData_2 + SimData_4; SimData_2*SimData_3*(k_6/(SimData_3 + k_12));
+  J6: SimData_1 + SimData_4 => SimData_1 + SimData_3; SimData_1*SimData_4*(k_5/(SimData_4 + k_11));
+
+// new code
+J1: newVar1               => SimData_1 + newVar1;             (newVar1 * I * k_1)/(k_7 + newVar1)              = j2 with added product
+J2: SimData_1             => newVar1;                         (SimData_1 * k_2)/(SimData_1 + k_8)               = j1
+J3: SimData_1 + newVar2   => SimData_1 + SimData_2 + newVar2; (newVar2 * SimData_1 * k_3)/(k_9 + newVar2)      = j4 with added product
+j4: SimData_2             => newVar2;                         (SimData_2 * k_4)(SimData_2 + k_10)              = j3
+j5: SimData_1 + newVar3   => SimData_1 + SimData_3 + newVar3; (newVar3 * SimData_1 * k_5)/(k_11 + newVar3)     = j6 with added product
+j6: SimData_2 + SimData_3 => SimData_2 + newVar3;             (SimData_3 * SimData_2 * k_6)/(SimData_3 + k_12) = j5 with added reactant6*
+
+model five
+
+// from moccasin
+j1: newVar1   =>  c + newVar1;               newVar1
+j2:           =>  c;                         v
+
+// my code
+j1: newVar1      =>  c + newVar1;               newVar1                = J1
+j2:              =>  c;                         v                      = j2
+j3: newVar2      =>  2*newVar2;                 newVar2
+j4:              =>  newVar1;                   k
+j5: a + newVar1  =>  a;                         a
+j6: newVar1      =>                             b
+
+model 54
+
+// from moccasin
+j1: SimData_1               => SimData_6;                         SimData_1 * (k__2/(SimData_1 + k__8))
+j2: SimData_6               => SimData_1                          (k__1 * SimData_6 * I)/(SimData_6 + k__7)              
+j3: SimData_2               => SimData_5                          SimData_2 * (k__4/(SimData_2 + k__10))
+j4: SimData_1 + SimData_5   => SimData_1 + SimData_2              (SimData_1 * SimData_5 * k__3)/(SimData_5 + k__9)
+j5: SimData_2 + SimData_3   => SimData_2 + SimData_4              (SimData_2 * SimData_3 * k__6)/(SimData_3 + k__12)
+j6: SimData_1 + SimData_4   => SimData_1 + SimData_3              (SimData_1 * SimData_4 * k__5)/(SimData_4 + k__11)
+
+// my code
+j1: j2 newVar1 = SimData_6 with added product newVar1
+j2: j1 newVar1 = SimData_6
+j3: j4 newVar2 = SimData_5 with added product newVar2
+j4: j3 newVar2 = SimData_5 
+j5: j6 newVar3 = SimData_4 with added product newVar3
+j6: j5 newVar3 = SimData_4
+
+model 55
+
+// from moccasin
+
+convert function
+
+populateInitialodeInfo 
+populateodeInfo
+  ea-> detectHiddenSpecies
+        ->analyze True only add -x+y
+          ->analyzeNode
+            ->isVariableSpeciesOrParameter
+            ->isNumericalConstantOrConstantParameter
+          ->hasExpressionAlreadyBeanRecorded
+            ->expressionExists
+              ->matches*
+        ->reorderMinusExPlusYIteratively
+          ->replaceExpressioninNodeWithNode
+        ->analyze False
+          -> as above
+        ->substituteParameterForExpression
+          ->getUniqueNoParameterName
+  -> addodePair
+  -> populateTerms
+      ->createTerms
+  -> createAnalysisVectors
+      ->populateCoefficientVector
+      ->populateDerivativeVector
+populateReactionCoefficients
+  ->createInitialValues
+  ->analyzeCoefficient
+  ->analyzePosDerivative
+  ->analyzeNagDerivative
+reconstructModel
+  ->dealWithSpecies
+  ->createReactions
+  ->removeRules
+                    
+    TYPE_K_MINUS_X_MINUS_Y
+        analyzeNode k-x-y 
+        right: y 
+        left: k-x 
+          analyzeNode k-x 
+          right: x 
+          left: k 
+          values: k x dx type=k-x 
+        values: y dy type=k-x-y   
+  , TYPE_K_PLUS_V_MINUS_X_MINUS_Y
+        analyzeNode k+v-x-y 
+        right: y 
+        left: k+v-x 
+          analyzeNode: k+v-x 
+          right: x 
+          left: k+v 
+            analyzeNode: k+v 
+            right: v 
+            left: k 
+            return False
+        values: k x dx v type=k+v-x 
+      values: y dy type=k+v-x-y 
+  , TYPE_K_MINUS_X_PLUS_W_MINUS_Y
+  , TYPE_K_MINUS_X
+  , TYPE_K_PLUS_V_MINUS_X
+ne
\ No newline at end of file
diff --git a/src/sbml/conversion/test/test-data/valid_54_bio.c b/src/sbml/conversion/test/test-data/valid_54_bio.c
new file mode 100644
index 0000000000..daced134b6
--- /dev/null
+++ b/src/sbml/conversion/test/test-data/valid_54_bio.c
@@ -0,0 +1,103 @@
+#ifdef SIZE_DEFINITIONS
+#define N_METABS 6
+#define N_ODE_METABS 0
+#define N_INDEP_METABS 3
+#define N_COMPARTMENTS 1
+#define N_GLOBAL_PARAMS 14
+#define N_KIN_PARAMS 0
+#define N_REACTIONS 6
+
+#define N_ARRAY_SIZE_P  14	// number of parameters
+#define N_ARRAY_SIZE_X  3	// number of initials
+#define N_ARRAY_SIZE_Y  4	// number of assigned elements
+#define N_ARRAY_SIZE_XC 3	// number of x concentration
+#define N_ARRAY_SIZE_PC 0	// number of p concentration
+#define N_ARRAY_SIZE_YC 3	// number of y concentration
+#define N_ARRAY_SIZE_DX 3	// number of ODEs 
+#define N_ARRAY_SIZE_CT 3	// number of conserved totals
+
+#endif // SIZE_DEFINITIONS
+
+#ifdef TIME
+#define T  <set here a user name for the time variable> 
+#endif // TIME
+
+#ifdef NAME_ARRAYS
+const char* p_names[] = {"compartmentOne", "I", "k__1", "k__2", "k__3", "k__4", "k__5", "k__6", "k__7", "k__8", "k__9", "k__10", "k__11", "k__12",  "" };
+const char* x_names[] = {"SimData_1", "SimData_2", "SimData_3",  "" };
+const char* y_names[] = {"SimData_4", "SimData_5", "SimData_6", "t",  "" };
+const char* xc_names[] = {"SimData_1", "SimData_2", "SimData_3",  "" };
+const char* pc_names[] = { "" };
+const char* yc_names[] = {"SimData_4", "SimData_5", "SimData_6",  "" };
+const char* dx_names[] = {"ODE SimData_1", "ODE SimData_2", "ODE SimData_3",  "" };
+const char* ct_names[] = {"CT SimData_4", "CT SimData_5", "CT SimData_6",  "" };
+#endif // NAME_ARRAYS
+
+#ifdef INITIAL
+x[0] = 1;	//metabolite 'SimData_1': reactions
+x[1] = 1;	//metabolite 'SimData_2': reactions
+x[2] = 1;	//metabolite 'SimData_3': reactions
+#endif /* INITIAL */
+
+#ifdef FIXED
+ct[0] = 0.99999999999999978;	//ct[0] conserved total for 'SimData_4'
+ct[1] = 0.99999999999999978;	//ct[1] conserved total for 'SimData_5'
+ct[2] = 0.99999999999999978;	//ct[2] conserved total for 'SimData_6'
+p[0] = 1;	//compartment 'compartmentOne':fixed
+p[1] = 1;	//global quantity 'I':fixed
+p[2] = 4.2416;	//global quantity 'k__1':fixed
+p[3] = 5.9816;	//global quantity 'k__2':fixed
+p[4] = 0.1009;	//global quantity 'k__3':fixed
+p[5] = 1.1549;	//global quantity 'k__4':fixed
+p[6] = 1.3618;	//global quantity 'k__5':fixed
+p[7] = 1.4219;	//global quantity 'k__6':fixed
+p[8] = 0.0051;	//global quantity 'k__7':fixed
+p[9] = 0.0972;	//global quantity 'k__8':fixed
+p[10] = 0.0012;	//global quantity 'k__9':fixed
+p[11] = 56.8583;	//global quantity 'k__10':fixed
+p[12] = 0.0111;	//global quantity 'k__11':fixed
+p[13] = 0.0014;	//global quantity 'k__12':fixed
+#endif /* FIXED */
+
+#ifdef ASSIGNMENT
+y[0] = ct[0]-x[2];	//metabolite 'SimData_4': reactions
+y[1] = ct[1]-x[1];	//metabolite 'SimData_5': reactions
+y[2] = ct[2]-x[0];	//metabolite 'SimData_6': reactions
+y[3] = T;	//model entity 't':assignment
+x_c[0] = x[0]/p[0];	//concentration of metabolite 'SimData_1': reactions
+x_c[1] = x[1]/p[0];	//concentration of metabolite 'SimData_2': reactions
+x_c[2] = x[2]/p[0];	//concentration of metabolite 'SimData_3': reactions
+y_c[0] = y[0]/p[0];	//concentration of metabolite 'SimData_4': reactions
+y_c[1] = y[1]/p[0];	//concentration of metabolite 'SimData_5': reactions
+y_c[2] = y[2]/p[0];	//concentration of metabolite 'SimData_6': reactions
+#endif /* ASSIGNMENT */
+
+#ifdef FUNCTIONS_HEADERS
+double FunctionForJ1_1(double sub_0, double volume_0, double param_0, double param_1); 
+double FunctionForJ2_1(double param_0, double sub_0, double volume_0, double param_1, double param_2); 
+double FunctionForJ3(double sub_0, double volume_0, double param_0, double param_1); 
+double FunctionForJ4(double sub_0, double sub_1, double volume_0, double param_0, double param_1); 
+double FunctionForJ5(double sub_0, double sub_1, double volume_0, double param_0, double param_1); 
+double FunctionForJ6(double sub_0, double sub_1, double volume_0, double param_0, double param_1); 
+#endif /* FUNCTIONS_HEADERS */
+
+#ifdef FUNCTIONS
+double FunctionForJ1_1(double sub_0, double volume_0, double param_0, double param_1) 	//Function for J1_1
+{return  sub_0*(param_0/(sub_0+param_1))/volume_0;} 
+double FunctionForJ2_1(double param_0, double sub_0, double volume_0, double param_1, double param_2) 	//Function for J2_1
+{return  param_1*sub_0*(param_0/(sub_0+param_2))/volume_0;} 
+double FunctionForJ3(double sub_0, double volume_0, double param_0, double param_1) 	//Function for J3
+{return  sub_0*(param_1/(sub_0+param_0))/volume_0;} 
+double FunctionForJ4(double sub_0, double sub_1, double volume_0, double param_0, double param_1) 	//Function for J4
+{return  sub_0*sub_1*(param_0/(sub_1+param_1))/volume_0;} 
+double FunctionForJ5(double sub_0, double sub_1, double volume_0, double param_0, double param_1) 	//Function for J5
+{return  sub_0*sub_1*(param_1/(sub_1+param_0))/volume_0;} 
+double FunctionForJ6(double sub_0, double sub_1, double volume_0, double param_0, double param_1) 	//Function for J6
+{return  sub_0*sub_1*(param_1/(sub_1+param_0))/volume_0;} 
+#endif /* FUNCTIONS */
+
+#ifdef ODEs
+dx[0] = -FunctionForJ1_1(x_c[0], p[0], p[3], p[9])*p[0]+FunctionForJ2_1(p[1], y_c[2], p[0], p[2], p[8])*p[0];
+dx[1] = -FunctionForJ3(x_c[1], p[0], p[11], p[5])*p[0]+FunctionForJ4(x_c[0], y_c[1], p[0], p[4], p[10])*p[0];
+dx[2] = -FunctionForJ5(x_c[1], x_c[2], p[0], p[13], p[7])*p[0]+FunctionForJ6(x_c[0], y_c[0], p[0], p[12], p[6])*p[0];
+#endif /* ODEs */
diff --git a/src/sbml/conversion/test/test-data/valid_54_rr.c b/src/sbml/conversion/test/test-data/valid_54_rr.c
new file mode 100644
index 0000000000..d3ae5c036e
--- /dev/null
+++ b/src/sbml/conversion/test/test-data/valid_54_rr.c
@@ -0,0 +1,79 @@
+#ifdef SIZE_DEFINITIONS
+#define N_METABS 3
+#define N_ODE_METABS 3
+#define N_INDEP_METABS 0
+#define N_COMPARTMENTS 1
+#define N_GLOBAL_PARAMS 14
+#define N_KIN_PARAMS 0
+#define N_REACTIONS 0
+
+#define N_ARRAY_SIZE_P  14	// number of parameters
+#define N_ARRAY_SIZE_X  3	// number of initials
+#define N_ARRAY_SIZE_Y  1	// number of assigned elements
+#define N_ARRAY_SIZE_XC 3	// number of x concentration
+#define N_ARRAY_SIZE_PC 0	// number of p concentration
+#define N_ARRAY_SIZE_YC 0	// number of y concentration
+#define N_ARRAY_SIZE_DX 3	// number of ODEs 
+#define N_ARRAY_SIZE_CT 0	// number of conserved totals
+
+#endif // SIZE_DEFINITIONS
+
+#ifdef TIME
+#define T  <set here a user name for the time variable> 
+#endif // TIME
+
+#ifdef NAME_ARRAYS
+const char* p_names[] = {"compartmentOne", "I", "k__1", "k__2", "k__3", "k__4", "k__5", "k__6", "k__7", "k__8", "k__9", "k__10", "k__11", "k__12",  "" };
+const char* x_names[] = {"SimData__1", "SimData__2", "SimData__3",  "" };
+const char* y_names[] = {"t",  "" };
+const char* xc_names[] = {"SimData__1", "SimData__2", "SimData__3",  "" };
+const char* pc_names[] = { "" };
+const char* yc_names[] = { "" };
+const char* dx_names[] = {"ODE SimData__1", "ODE SimData__2", "ODE SimData__3",  "" };
+const char* ct_names[] = { "" };
+#endif // NAME_ARRAYS
+
+#ifdef INITIAL
+x[0] = 1;	//metabolite 'SimData__1': ode
+x[1] = 1;	//metabolite 'SimData__2': ode
+x[2] = 1;	//metabolite 'SimData__3': ode
+#endif /* INITIAL */
+
+#ifdef FIXED
+p[0] = 1;	//compartment 'compartmentOne':fixed
+p[1] = 1;	//global quantity 'I':fixed
+p[2] = 4.2416;	//global quantity 'k__1':fixed
+p[3] = 5.9816;	//global quantity 'k__2':fixed
+p[4] = 0.1009;	//global quantity 'k__3':fixed
+p[5] = 1.1549;	//global quantity 'k__4':fixed
+p[6] = 1.3618;	//global quantity 'k__5':fixed
+p[7] = 1.4219;	//global quantity 'k__6':fixed
+p[8] = 0.0051;	//global quantity 'k__7':fixed
+p[9] = 0.0972;	//global quantity 'k__8':fixed
+p[10] = 0.0012;	//global quantity 'k__9':fixed
+p[11] = 56.8583;	//global quantity 'k__10':fixed
+p[12] = 0.0111;	//global quantity 'k__11':fixed
+p[13] = 0.0014;	//global quantity 'k__12':fixed
+#endif /* FIXED */
+
+#ifdef ASSIGNMENT
+y[0] = T;	//model entity 't':assignment
+x_c[0] = x[0]/p[0];	//concentration of metabolite 'SimData__1': ode
+x_c[1] = x[1]/p[0];	//concentration of metabolite 'SimData__2': ode
+x_c[2] = x[2]/p[0];	//concentration of metabolite 'SimData__3': ode
+#endif /* ASSIGNMENT */
+
+#ifdef FUNCTIONS_HEADERS
+#endif /* FUNCTIONS_HEADERS */
+
+#ifdef FUNCTIONS
+#endif /* FUNCTIONS */
+
+#ifdef ODEs
+dx[0] = p[2]*p[1]*(1.00000000000000000-x_c[0])/(1.00000000000000000-x_c[0]+p[8])-p[3]*x_c[0]/(x_c[0]+p[9]) * p[0];	//model entity 'SimData__1':ode
+dx[0] = p[2]*p[1]*(1.00000000000000000-x_c[0])/(1.00000000000000000-x_c[0]+p[8])-p[3]*x_c[0]/(x_c[0]+p[9]) * p[0];	//model entity 'SimData__1':ode
+dx[1] = p[4]*x_c[0]*(1.00000000000000000-x_c[1])/(1.00000000000000000-x_c[1]+p[10])-p[5]*x_c[1]/(x_c[1]+p[11]) * p[0];	//model entity 'SimData__2':ode
+dx[1] = p[4]*x_c[0]*(1.00000000000000000-x_c[1])/(1.00000000000000000-x_c[1]+p[10])-p[5]*x_c[1]/(x_c[1]+p[11]) * p[0];	//model entity 'SimData__2':ode
+dx[2] = p[6]*x_c[0]*(1.00000000000000000-x_c[2])/(1.00000000000000000-x_c[2]+p[12])-p[7]*x_c[1]*x_c[2]/(x_c[2]+p[13]) * p[0];	//model entity 'SimData__3':ode
+dx[2] = p[6]*x_c[0]*(1.00000000000000000-x_c[2])/(1.00000000000000000-x_c[2]+p[12])-p[7]*x_c[1]*x_c[2]/(x_c[2]+p[13]) * p[0];	//model entity 'SimData__3':ode
+#endif /* ODEs */