diff --git a/modcc/expression.cpp b/modcc/expression.cpp
index 06a70e0c99bb7a9d9c38d472f2b2382e0f44371d..e2694ea6f36e5dc33622e9d4fbe725fbd13d24d4 100644
--- a/modcc/expression.cpp
+++ b/modcc/expression.cpp
@@ -110,16 +110,12 @@ expression_ptr IdentifierExpression::clone() const {
     return make_expression<IdentifierExpression>(location_, spelling_);
 }
 
-bool IdentifierExpression::is_lvalue() {
-    // check for global variable that is writeable
-    auto var = symbol_->is_variable();
-    if(var) return var->is_writeable();
-
-    // else look for local symbol
-    if( symbol_->kind() == symbolKind::local_variable ) {
-        return true;
-    }
+bool IdentifierExpression::is_lvalue() const {
+    return is_global_lvalue() || symbol_->kind() == symbolKind::local_variable;
+}
 
+bool IdentifierExpression::is_global_lvalue() const {
+    if(auto var = symbol_->is_variable()) return var->is_writeable();
     return false;
 }
 
@@ -131,6 +127,14 @@ expression_ptr NumberExpression::clone() const {
     return make_expression<NumberExpression>(location_, value_);
 }
 
+/*******************************************************************************
+  IntegerExpression
+********************************************************************************/
+
+expression_ptr IntegerExpression::clone() const {
+    return make_expression<IntegerExpression>(location_, integer_);
+}
+
 /*******************************************************************************
   LocalDeclaration
 *******************************************************************************/
@@ -250,6 +254,98 @@ std::string IndexedVariable::to_string() const {
         + ", ion" + (ion_channel()==ionKind::none ? red(ch) : green(ch)) + ") ";
 }
 
+/*******************************************************************************
+  ReactionExpression
+*******************************************************************************/
+
+std::string ReactionExpression::to_string() const {
+    return blue("reaction") +
+           pprintf(" % <-> % (%, %)",
+               lhs()->to_string(), rhs()->to_string(),
+                fwd_rate()->to_string(), rev_rate()->to_string());
+}
+
+expression_ptr ReactionExpression::clone() const {
+    return make_expression<ReactionExpression>(
+        location_, lhs()->clone(), rhs()->clone(), fwd_rate()->clone(), rev_rate()->clone());
+}
+
+void ReactionExpression::semantic(std::shared_ptr<scope_type> scp) {
+    scope_ = scp;
+    lhs()->semantic(scp);
+    rhs()->semantic(scp);
+
+    fwd_rate()->semantic(scp);
+    rev_rate()->semantic(scp);
+    if(fwd_rate_->is_procedure_call() || rev_rate_->is_procedure_call()) {
+        error("procedure calls can't be made in an expression");
+    }
+}
+
+/*******************************************************************************
+  StoichTermExpression
+*******************************************************************************/
+
+expression_ptr StoichTermExpression::clone() const {
+    return make_expression<StoichTermExpression>(
+        location_, coeff()->clone(), ident()->clone());
+}
+
+void StoichTermExpression::semantic(std::shared_ptr<scope_type> scp) {
+    scope_ = scp;
+    ident()->semantic(scp);
+}
+
+/*******************************************************************************
+  StoichExpression
+*******************************************************************************/
+
+expression_ptr StoichExpression::clone() const {
+    std::vector<expression_ptr> cloned_terms;
+    for(auto& e: terms()) {
+        cloned_terms.emplace_back(e->clone());
+    }
+
+    return make_expression<StoichExpression>(location_, std::move(cloned_terms));
+}
+
+std::string StoichExpression::to_string() const {
+    std::string s;
+    bool first = true;
+    for(auto& e: terms()) {
+        if (!first) s += "+";
+        s += e->to_string();
+        first = false;
+    }
+    return s;
+}
+
+void StoichExpression::semantic(std::shared_ptr<scope_type> scp) {
+    scope_ = scp;
+    for(auto& e: terms()) {
+        e->semantic(scp);
+    }
+}
+
+/*******************************************************************************
+  ConserveExpression
+*******************************************************************************/
+
+expression_ptr ConserveExpression::clone() const {
+    return make_expression<ConserveExpression>(
+        location_, lhs()->clone(), rhs()->clone());
+}
+
+void ConserveExpression::semantic(std::shared_ptr<scope_type> scp) {
+    scope_ = scp;
+    lhs_->semantic(scp);
+    rhs_->semantic(scp);
+
+    if(rhs_->is_procedure_call()) {
+        error("procedure calls can't be made in an expression");
+    }
+}
+
 /*******************************************************************************
   CallExpression
 *******************************************************************************/
@@ -709,7 +805,6 @@ expression_ptr IfExpression::clone() const {
 void Expression::accept(Visitor *v) {
     v->visit(this);
 }
-
 void Symbol::accept(Visitor *v) {
     v->visit(this);
 }
@@ -746,6 +841,9 @@ void IndexedVariable::accept(Visitor *v) {
 void NumberExpression::accept(Visitor *v) {
     v->visit(this);
 }
+void IntegerExpression::accept(Visitor *v) {
+    v->visit(this);
+}
 void LocalDeclaration::accept(Visitor *v) {
     v->visit(this);
 }
@@ -794,6 +892,18 @@ void BinaryExpression::accept(Visitor *v) {
 void AssignmentExpression::accept(Visitor *v) {
     v->visit(this);
 }
+void ConserveExpression::accept(Visitor *v) {
+    v->visit(this);
+}
+void ReactionExpression::accept(Visitor *v) {
+    v->visit(this);
+}
+void StoichExpression::accept(Visitor *v) {
+    v->visit(this);
+}
+void StoichTermExpression::accept(Visitor *v) {
+    v->visit(this);
+}
 void AddBinaryExpression::accept(Visitor *v) {
     v->visit(this);
 }
diff --git a/modcc/expression.hpp b/modcc/expression.hpp
index b1076b5b4a59e50f12aa936207ac0f3ddf9097bb..a7d84e61a6a72220e4bbc34e692de703e1681a94 100644
--- a/modcc/expression.hpp
+++ b/modcc/expression.hpp
@@ -23,6 +23,7 @@ class IfExpression;
 class VariableExpression;
 class IndexedVariable;
 class NumberExpression;
+class IntegerExpression;
 class LocalDeclaration;
 class ArgumentExpression;
 class DerivativeExpression;
@@ -40,6 +41,10 @@ class CosUnaryExpression;
 class SinUnaryExpression;
 class BinaryExpression;
 class AssignmentExpression;
+class ReactionExpression;
+class StoichExpression;
+class StoichTermExpression;
+class ConserveExpression;
 class AddBinaryExpression;
 class SubBinaryExpression;
 class MulBinaryExpression;
@@ -77,6 +82,7 @@ enum class procedureKind {
     initial,     ///< INITIAL
     net_receive, ///< NET_RECEIVE
     breakpoint,  ///< BREAKPOINT
+    kinetic,     ///< KINETIC
     derivative   ///< DERIVATIVE
 };
 std::string to_string(procedureKind k);
@@ -157,16 +163,22 @@ public:
     virtual PrototypeExpression*   is_prototype()         {return nullptr;}
     virtual IdentifierExpression*  is_identifier()        {return nullptr;}
     virtual NumberExpression*      is_number()            {return nullptr;}
+    virtual IntegerExpression*     is_integer()           {return nullptr;}
     virtual BinaryExpression*      is_binary()            {return nullptr;}
     virtual UnaryExpression*       is_unary()             {return nullptr;}
     virtual AssignmentExpression*  is_assignment()        {return nullptr;}
+    virtual ConserveExpression*    is_conserve()          {return nullptr;}
+    virtual ReactionExpression*    is_reaction()          {return nullptr;}
+    virtual StoichExpression*      is_stoich()            {return nullptr;}
+    virtual StoichTermExpression*  is_stoich_term()       {return nullptr;}
     virtual ConditionalExpression* is_conditional()       {return nullptr;}
     virtual InitialBlock*          is_initial_block()     {return nullptr;}
     virtual SolveExpression*       is_solve_statement()   {return nullptr;}
     virtual Symbol*                is_symbol()            {return nullptr;}
     virtual ConductanceExpression* is_conductance_statement() {return nullptr;}
 
-    virtual bool is_lvalue() {return false;}
+    virtual bool is_lvalue() const {return false;}
+    virtual bool is_global_lvalue() const {return false;}
 
     // force all derived classes to implement visitor
     // this might be a bad idea
@@ -259,7 +271,8 @@ public:
 
     IdentifierExpression* is_identifier() override {return this;}
 
-    bool is_lvalue() override;
+    bool is_lvalue() const override;
+    bool is_global_lvalue() const override;
 
     ~IdentifierExpression() {}
 
@@ -298,14 +311,14 @@ public:
 class NumberExpression : public Expression {
 public:
     NumberExpression(Location loc, std::string const& value)
-        : Expression(loc), value_(std::stod(value))
+        : Expression(loc), value_(std::stold(value))
     {}
 
     NumberExpression(Location loc, long double value)
         : Expression(loc), value_(value)
     {}
 
-    long double value() const {return value_;};
+    virtual long double value() const {return value_;};
 
     std::string to_string() const override {
         return purple(pprintf("%", value_));
@@ -324,6 +337,38 @@ private:
     long double value_;
 };
 
+// an integral number
+class IntegerExpression : public NumberExpression {
+public:
+    IntegerExpression(Location loc, std::string const& value)
+        : NumberExpression(loc, value), integer_(std::stoll(value))
+    {}
+
+    IntegerExpression(Location loc, long long integer)
+        : NumberExpression(loc, static_cast<long double>(integer)), integer_(integer)
+    {}
+
+    long long integer_value() const {return integer_;}
+
+    std::string to_string() const override {
+        return purple(pprintf("%", integer_));
+    }
+
+    // do nothing for number semantic analysis
+    void semantic(std::shared_ptr<scope_type> scp) override {};
+    expression_ptr clone() const override;
+
+    IntegerExpression* is_integer() override {return this;}
+
+    ~IntegerExpression() {}
+
+    void accept(Visitor *v) override;
+private:
+    long long integer_;
+};
+
+
+
 // declaration of a LOCAL variable
 class LocalDeclaration : public Expression {
 public:
@@ -788,6 +833,102 @@ private:
     std::vector<expression_ptr> args_;
 };
 
+class ReactionExpression : public Expression {
+public:
+    ReactionExpression(Location loc,
+                       expression_ptr&& lhs_terms,
+                       expression_ptr&& rhs_terms,
+                       expression_ptr&& fwd_rate_expr,
+                       expression_ptr&& rev_rate_expr)
+    : Expression(loc),
+      lhs_(std::move(lhs_terms)), rhs_(std::move(rhs_terms)),
+      fwd_rate_(std::move(fwd_rate_expr)), rev_rate_(std::move(rev_rate_expr))
+    {}
+
+    ReactionExpression* is_reaction() override {return this;}
+
+    std::string to_string() const override;
+    void semantic(std::shared_ptr<scope_type> scp) override;
+    expression_ptr clone() const override;
+    void accept(Visitor *v) override;
+
+    expression_ptr& lhs() { return lhs_; }
+    const expression_ptr& lhs() const { return lhs_; }
+
+    expression_ptr& rhs() { return rhs_; }
+    const expression_ptr& rhs() const { return rhs_; }
+
+    expression_ptr& fwd_rate() { return fwd_rate_; }
+    const expression_ptr& fwd_rate() const { return fwd_rate_; }
+
+    expression_ptr& rev_rate() { return rev_rate_; }
+    const expression_ptr& rev_rate() const { return rev_rate_; }
+
+private:
+    expression_ptr lhs_;
+    expression_ptr rhs_;
+    expression_ptr fwd_rate_;
+    expression_ptr rev_rate_;
+};
+
+class StoichTermExpression : public Expression {
+public:
+    StoichTermExpression(Location loc,
+                         expression_ptr&& coeff,
+                         expression_ptr&& ident)
+    : Expression(loc),
+      coeff_(std::move(coeff)), ident_(std::move(ident))
+    {}
+
+    StoichTermExpression* is_stoich_term() override {return this;}
+
+    std::string to_string() const override {
+        return pprintf("% %", coeff()->to_string(), ident()->to_string());
+    }
+    void semantic(std::shared_ptr<scope_type> scp) override;
+    expression_ptr clone() const override;
+    void accept(Visitor *v) override;
+
+    expression_ptr& coeff() { return coeff_; }
+    const expression_ptr& coeff() const { return coeff_; }
+
+    expression_ptr& ident() { return ident_; }
+    const expression_ptr& ident() const { return ident_; }
+
+    bool negative() const {
+        auto iexpr = coeff_->is_integer();
+        return iexpr && iexpr->integer_value()<0;
+    }
+
+private:
+    expression_ptr coeff_;
+    expression_ptr ident_;
+};
+
+class StoichExpression : public Expression {
+public:
+    StoichExpression(Location loc, std::vector<expression_ptr>&& terms)
+    : Expression(loc), terms_(std::move(terms))
+    {}
+
+    StoichExpression(Location loc)
+    : Expression(loc)
+    {}
+
+    StoichExpression* is_stoich() override {return this;}
+
+    std::string to_string() const override;
+    void semantic(std::shared_ptr<scope_type> scp) override;
+    expression_ptr clone() const override;
+    void accept(Visitor *v) override;
+
+    std::vector<expression_ptr>& terms() { return terms_; }
+    const std::vector<expression_ptr>& terms() const { return terms_; }
+
+private:
+    std::vector<expression_ptr> terms_;
+};
+
 // marks a call site in the AST
 // is used to mark both function and procedure calls
 class CallExpression : public Expression {
@@ -1100,6 +1241,20 @@ public:
     void accept(Visitor *v) override;
 };
 
+class ConserveExpression : public BinaryExpression {
+public:
+    ConserveExpression(Location loc, expression_ptr&& lhs, expression_ptr&& rhs)
+    :   BinaryExpression(loc, tok::eq, std::move(lhs), std::move(rhs))
+    {}
+
+    ConserveExpression* is_conserve() override {return this;}
+    expression_ptr clone() const override;
+
+    void semantic(std::shared_ptr<scope_type> scp) override;
+
+    void accept(Visitor *v) override;
+};
+
 class AddBinaryExpression : public BinaryExpression {
 public:
     AddBinaryExpression(Location loc, expression_ptr&& lhs, expression_ptr&& rhs)
diff --git a/modcc/lexer.cpp b/modcc/lexer.cpp
index 0a85de24500d1f846c1f0156f863908453cca961..b15513a93690006e4696839fd9e158fa356c6419 100644
--- a/modcc/lexer.cpp
+++ b/modcc/lexer.cpp
@@ -28,6 +28,9 @@ inline bool is_eof(char c) {
 inline bool is_operator(char c) {
     return (c=='+' || c=='-' || c=='*' || c=='/' || c=='^' || c=='\'');
 }
+inline bool is_plusminus(char c) {
+    return (c=='+' || c=='-');
+}
 
 //*********************
 // Lexer
@@ -90,11 +93,7 @@ Token Lexer::parse() {
             // number
             case '0' ... '9':
             case '.':
-                t.spelling = number();
-
-                // test for error when reading number
-                t.type = (status_==lexerStatus::error) ? tok::reserved : tok::number;
-                return t;
+                return number();
 
             // identifier or keyword
             case 'a' ... 'z':
@@ -123,6 +122,10 @@ Token Lexer::parse() {
                 t.type = tok::rbrace;
                 t.spelling += character();
                 return t;
+            case '~':
+                t.type = tok::tilde;
+                t.spelling += character();
+                return t;
             case '=': {
                 t.spelling += character();
                 if(*current_=='=') {
@@ -172,6 +175,11 @@ Token Lexer::parse() {
                     t.spelling += character();
                     t.type = tok::lte;
                 }
+                else if(*current_=='-' && current_[1]=='>') {
+                    t.spelling += character();
+                    t.spelling += character();
+                    t.type = tok::arrow;
+                }
                 else {
                     t.type = tok::lt;
                 }
@@ -228,7 +236,7 @@ Token Lexer::peek() {
 }
 
 // scan floating point number from stream
-std::string Lexer::number() {
+Token Lexer::number() {
     std::string str;
     char c = *current_;
 
@@ -253,13 +261,21 @@ std::string Lexer::number() {
                 incorrectly_formed_mantisa = true;
             }
         }
-        else if(c=='e' || c=='E') {
-            uses_scientific_notation++;
-            str += c;
-            current_++;
-            // Consume the next char if +/-
-            if (*current_ == '+' || *current_ == '-') {
-                str += *current_++;
+        else if(!uses_scientific_notation && (c=='e' || c=='E')) {
+            if(is_numeric(current_[1]) ||
+               is_plusminus(current_[1]) && is_numeric(current_[2]))
+            {
+                uses_scientific_notation++;
+                str += c;
+                current_++;
+                // Consume the next char if +/-
+                if (is_plusminus(*current_)) {
+                    str += *current_++;
+                }
+            }
+            else {
+                // the 'e' or 'E' is the beginning of a new token
+                break;
             }
         }
         else {
@@ -278,13 +294,19 @@ std::string Lexer::number() {
         error_string_ = pprintf("too many .'s when reading the number '%'", yellow(str));
         status_ = lexerStatus::error;
     }
-    // check that e or E is not used more than once in the number
-    if(uses_scientific_notation>1) {
-        error_string_ = pprintf("can't parse the number '%'", yellow(str));
-        status_ = lexerStatus::error;
+
+    tok type;
+    if(status_==lexerStatus::error) {
+        type = tok::reserved;
+    }
+    else if(num_point<1 && !uses_scientific_notation) {
+        type = tok::integer;
+    }
+    else {
+        type = tok::real;
     }
 
-    return str;
+    return Token(type, str, location_);
 }
 
 // scan identifier from stream
diff --git a/modcc/lexer.hpp b/modcc/lexer.hpp
index 52917e94477a09173d7e580ea067f9e7b691f2ed..142a5d3b1f7483bd548eb0a377e70ced5e9f2a19 100644
--- a/modcc/lexer.hpp
+++ b/modcc/lexer.hpp
@@ -73,7 +73,7 @@ public:
     Token peek();
 
     // scan a number from the stream
-    std::string number();
+    Token number();
 
     // scan an identifier string from the stream
     std::string identifier();
diff --git a/modcc/parser.cpp b/modcc/parser.cpp
index ffa498c08f50e4b0b1c2d7a6a3b8370c4e8c38de..43855a7ef1ad0915e7312aea048c23448015e06a 100644
--- a/modcc/parser.cpp
+++ b/modcc/parser.cpp
@@ -104,11 +104,12 @@ bool Parser::parse() {
             case tok::assigned :
                 parse_assigned_block();
                 break;
-            // INITIAL, DERIVATIVE, PROCEDURE, NET_RECEIVE and BREAKPOINT blocks
+            // INITIAL, KINETIC, DERIVATIVE, PROCEDURE, NET_RECEIVE and BREAKPOINT blocks
             // are all lowered to ProcedureExpression
             case tok::net_receive:
             case tok::breakpoint :
             case tok::initial    :
+            case tok::kinetic    :
             case tok::derivative :
             case tok::procedure  :
                 {
@@ -169,7 +170,7 @@ std::vector<Token> Parser::comma_separated_identifiers() {
             error(pprintf("found keyword '%', expected a variable name", token_.spelling));
             return tokens;
         }
-        else if(token_.type == tok::number) {
+        else if(token_.type == tok::real || token_.type == tok::integer) {
             error(pprintf("found number '%', expected a variable name", token_.spelling));
             return tokens;
         }
@@ -490,7 +491,7 @@ void Parser::parse_parameter_block() {
                 parm.value = "-";
                 get_token();
             }
-            if(token_.type != tok::number) {
+            if(token_.type != tok::integer && token_.type != tok::real) {
                 success = 0;
                 goto parm_exit;
             }
@@ -595,7 +596,7 @@ ass_exit:
 }
 
 std::vector<Token> Parser::unit_description() {
-    static const tok legal_tokens[] = {tok::identifier, tok::divide, tok::number};
+    static const tok legal_tokens[] = {tok::identifier, tok::divide, tok::real, tok::integer};
     int startline = location_.line;
     std::vector<Token> tokens;
 
@@ -726,6 +727,12 @@ symbol_ptr Parser::parse_procedure() {
             if( !expect( tok::identifier ) ) return nullptr;
             p = parse_prototype();
             break;
+        case tok::kinetic:
+            kind = procedureKind::kinetic;
+            get_token(); // consume keyword token
+            if( !expect( tok::identifier ) ) return nullptr;
+            p = parse_prototype();
+            break;
         case tok::procedure:
             kind = procedureKind::normal;
             get_token(); // consume keyword token
@@ -810,6 +817,10 @@ expression_ptr Parser::parse_statement() {
             return parse_local();
         case tok::identifier :
             return parse_line_expression();
+        case tok::conserve :
+            return parse_conserve_expression();
+        case tok::tilde :
+            return parse_reaction_expression();
         case tok::initial :
             // only used for INITIAL block in NET_RECEIVE
             return parse_initial();
@@ -938,6 +949,152 @@ expression_ptr Parser::parse_line_expression() {
     return lhs;
 }
 
+expression_ptr Parser::parse_stoich_term() {
+    expression_ptr coeff = make_expression<IntegerExpression>(location_, 1);
+    auto here = location_;
+    bool negative = false;
+
+    while(token_.type==tok::minus) {
+        negative = !negative;
+        get_token(); // consume '-'
+    }
+
+    if(token_.type==tok::integer) {
+        coeff = parse_integer();
+    }
+
+    if(token_.type!=tok::identifier) {
+        error(pprintf("expected an identifier, found '%'", yellow(token_.spelling)));
+        return nullptr;
+    }
+
+    if(negative) {
+        coeff = make_expression<IntegerExpression>(here, -coeff->is_integer()->integer_value());
+    }
+    return make_expression<StoichTermExpression>(here, std::move(coeff), parse_identifier());
+}
+
+expression_ptr Parser::parse_stoich_expression() {
+    std::vector<expression_ptr> terms;
+    auto here = location_;
+
+    if(token_.type==tok::integer || token_.type==tok::identifier || token_.type==tok::minus) {
+        auto term = parse_stoich_term();
+        if (!term) return nullptr;
+
+        terms.push_back(std::move(term));
+
+        while(token_.type==tok::plus || token_.type==tok::minus) {
+            if (token_.type==tok::plus) {
+                get_token(); // consume plus
+            }
+
+            auto term = parse_stoich_term();
+            if (!term) return nullptr;
+
+            terms.push_back(std::move(term));
+        }
+    }
+
+    return make_expression<StoichExpression>(here, std::move(terms));
+}
+
+expression_ptr Parser::parse_reaction_expression() {
+    auto here = location_;
+
+    if(token_.type!=tok::tilde) {
+        error(pprintf("expected '%', found '%'", yellow("~"), yellow(token_.spelling)));
+        return nullptr;
+    }
+
+    get_token(); // consume tilde
+    expression_ptr lhs = parse_stoich_expression();
+    if (!lhs) return nullptr;
+
+    // reaction halves must comprise non-negative terms
+    for (const auto& term: lhs->is_stoich()->terms()) {
+        // should always be true
+        if (auto sterm = term->is_stoich_term()) {
+            if (sterm->negative()) {
+                error(pprintf("expected only non-negative terms in reaction lhs, found '%'",
+                    yellow(term->to_string())));
+                return nullptr;
+            }
+        }
+    }
+
+    if(token_.type != tok::arrow) {
+        error(pprintf("expected '%', found '%'", yellow("<->"), yellow(token_.spelling)));
+        return nullptr;
+    }
+
+    get_token(); // consume arrow
+    expression_ptr rhs = parse_stoich_expression();
+    if (!rhs) return nullptr;
+
+    for (const auto& term: rhs->is_stoich()->terms()) {
+        // should always be true
+        if (auto sterm = term->is_stoich_term()) {
+            if (sterm->negative()) {
+                error(pprintf("expected only non-negative terms in reaction rhs, found '%'",
+                    yellow(term->to_string())));
+                return nullptr;
+            }
+        }
+    }
+
+    if(token_.type != tok::lparen) {
+        error(pprintf("expected '%', found '%'", yellow("("), yellow(token_.spelling)));
+        return nullptr;
+    }
+
+    get_token(); // consume lparen
+    expression_ptr fwd = parse_expression();
+    if (!fwd) return nullptr;
+
+    if(token_.type != tok::comma) {
+        error(pprintf("expected '%', found '%'", yellow(","), yellow(token_.spelling)));
+        return nullptr;
+    }
+
+    get_token(); // consume comma
+    expression_ptr rev = parse_expression();
+    if (!rev) return nullptr;
+
+    if(token_.type != tok::rparen) {
+        error(pprintf("expected '%', found '%'", yellow(")"), yellow(token_.spelling)));
+        return nullptr;
+    }
+
+    get_token(); // consume rparen
+    return make_expression<ReactionExpression>(here, std::move(lhs), std::move(rhs),
+        std::move(fwd), std::move(rev));
+}
+
+expression_ptr Parser::parse_conserve_expression() {
+    auto here = location_;
+
+    if(token_.type!=tok::conserve) {
+        error(pprintf("expected '%', found '%'", yellow("CONSERVE"), yellow(token_.spelling)));
+        return nullptr;
+    }
+
+    get_token(); // consume 'CONSERVE'
+    auto lhs = parse_stoich_expression();
+    if (!lhs) return nullptr;
+
+    if(token_.type != tok::eq) {
+        error(pprintf("expected '%', found '%'", yellow("="), yellow(token_.spelling)));
+        return nullptr;
+    }
+
+    get_token(); // consume '='
+    auto rhs = parse_expression();
+    if (!rhs) return nullptr;
+
+    return make_expression<ConserveExpression>(here, std::move(lhs), std::move(rhs));
+}
+
 expression_ptr Parser::parse_expression() {
     auto lhs = parse_unaryop();
 
@@ -1005,8 +1162,10 @@ expression_ptr Parser::parse_unaryop() {
 ///  ::  parenthesis expression (parsed recursively)
 expression_ptr Parser::parse_primary() {
     switch(token_.type) {
-        case tok::number:
-            return parse_number();
+        case tok::real:
+            return parse_real();
+        case tok::integer:
+            return parse_integer();
         case tok::identifier:
             if( peek().type == tok::lparen ) {
                 return parse_call();
@@ -1043,11 +1202,15 @@ expression_ptr Parser::parse_parenthesis_expression() {
     return e;
 }
 
-expression_ptr Parser::parse_number() {
+expression_ptr Parser::parse_real() {
     auto e = make_expression<NumberExpression>(token_.location, token_.spelling);
-
     get_token(); // consume the number
+    return e;
+}
 
+expression_ptr Parser::parse_integer() {
+    auto e = make_expression<IntegerExpression>(token_.location, token_.spelling);
+    get_token(); // consume the number
     return e;
 }
 
@@ -1280,6 +1443,10 @@ expression_ptr Parser::parse_block(bool is_nested) {
                 error("LOCAL variable declarations are not allowed inside a nested scope");
                 return nullptr;
             }
+            if(e->is_reaction()) {
+                error("reaction expressions are not allowed inside a nested scope");
+                return nullptr;
+            }
         }
 
         body.emplace_back(std::move(e));
diff --git a/modcc/parser.hpp b/modcc/parser.hpp
index b55bede17e698aa9d790211ddbee5a54f90ba56a..dc673e9d7b9f9d1084f8d024801ec4eaf87b9d88 100644
--- a/modcc/parser.hpp
+++ b/modcc/parser.hpp
@@ -17,12 +17,17 @@ public:
     expression_ptr parse_prototype(std::string);
     expression_ptr parse_statement();
     expression_ptr parse_identifier();
-    expression_ptr parse_number();
+    expression_ptr parse_integer();
+    expression_ptr parse_real();
     expression_ptr parse_call();
     expression_ptr parse_expression();
     expression_ptr parse_primary();
     expression_ptr parse_parenthesis_expression();
     expression_ptr parse_line_expression();
+    expression_ptr parse_stoich_expression();
+    expression_ptr parse_stoich_term();
+    expression_ptr parse_reaction_expression();
+    expression_ptr parse_conserve_expression();
     expression_ptr parse_binop(expression_ptr&&, Token);
     expression_ptr parse_unaryop();
     expression_ptr parse_local();
diff --git a/modcc/token.cpp b/modcc/token.cpp
index 1216c91287a7be1df4a0311c65c8436bdb4fb1cb..7383ca9fc684b4f47122a98e8f1af035724c8d9a 100644
--- a/modcc/token.cpp
+++ b/modcc/token.cpp
@@ -29,6 +29,7 @@ static Keyword keywords[] = {
     {"STATE",       tok::state},
     {"BREAKPOINT",  tok::breakpoint},
     {"DERIVATIVE",  tok::derivative},
+    {"KINETIC",     tok::kinetic},
     {"PROCEDURE",   tok::procedure},
     {"FUNCTION",    tok::function},
     {"INITIAL",     tok::initial},
@@ -42,6 +43,7 @@ static Keyword keywords[] = {
     {"WRITE",       tok::write},
     {"RANGE",       tok::range},
     {"LOCAL",       tok::local},
+    {"CONSERVE",    tok::conserve},
     {"SOLVE",       tok::solve},
     {"THREADSAFE",  tok::threadsafe},
     {"GLOBAL",      tok::global},
@@ -72,6 +74,8 @@ static TokenString token_strings[] = {
     {">=",          tok::gte},
     {"==",          tok::equality},
     {"!=",          tok::ne},
+    {"<->",         tok::arrow},
+    {"~",           tok::tilde},
     {",",           tok::comma},
     {"'",           tok::prime},
     {"{",           tok::lbrace},
@@ -79,7 +83,8 @@ static TokenString token_strings[] = {
     {"(",           tok::lparen},
     {")",           tok::rparen},
     {"identifier",  tok::identifier},
-    {"number",      tok::number},
+    {"real",        tok::real},
+    {"integer",     tok::integer},
     {"TITLE",       tok::title},
     {"NEURON",      tok::neuron},
     {"UNITS",       tok::units},
@@ -88,6 +93,7 @@ static TokenString token_strings[] = {
     {"STATE",       tok::state},
     {"BREAKPOINT",  tok::breakpoint},
     {"DERIVATIVE",  tok::derivative},
+    {"KINETIC",     tok::kinetic},
     {"PROCEDURE",   tok::procedure},
     {"FUNCTION",    tok::function},
     {"INITIAL",     tok::initial},
diff --git a/modcc/token.hpp b/modcc/token.hpp
index 8e97e47f9fa8a46bed23913a4762a90324bf60e2..c31cdbc409a2dd411cc5c074bd351fd994d4d029 100644
--- a/modcc/token.hpp
+++ b/modcc/token.hpp
@@ -23,6 +23,12 @@ enum class tok {
     equality,// ==
     ne,      // !=
 
+    // <->
+    arrow,
+
+    // ~
+    tilde,
+
     // , '
     comma, prime,
 
@@ -35,7 +41,7 @@ enum class tok {
     identifier,
 
     // numbers
-    number,
+    real, integer,
 
     /////////////////////////////
     // keywords
@@ -44,14 +50,14 @@ enum class tok {
     title,
     neuron, units, parameter,
     assigned, state, breakpoint,
-    derivative, procedure, initial, function,
+    derivative, kinetic, procedure, initial, function,
     net_receive,
 
     // keywoards inside blocks
     unitsoff, unitson,
     suffix, nonspecific_current, useion,
     read, write,
-    range, local,
+    range, local, conserve,
     solve, method,
     threadsafe, global,
     point_process,
@@ -76,7 +82,7 @@ enum class tok {
 struct Token {
     // the spelling string contains the text of the token as it was written
     // in the input file
-    //   type = tok::number     : spelling = "3.1415"  (e.g.)
+    //   type = tok::real       : spelling = "3.1415"  (e.g.)
     //   type = tok::identifier : spelling = "foo_bar" (e.g.)
     //   type = tok::plus       : spelling = "+"       (always)
     //   type = tok::if         : spelling = "if"      (always)
diff --git a/modcc/visitor.hpp b/modcc/visitor.hpp
index faab5cca2500e7a144263a3402364284860bc227..ae5c42dcd29f7760a0b8dcc533fdb80e9561ed0f 100644
--- a/modcc/visitor.hpp
+++ b/modcc/visitor.hpp
@@ -20,6 +20,7 @@ public:
     virtual void visit(LocalVariable *e)        { visit((Expression*) e); }
     virtual void visit(IdentifierExpression *e) { visit((Expression*) e); }
     virtual void visit(NumberExpression *e)     { visit((Expression*) e); }
+    virtual void visit(IntegerExpression *e)    { visit((NumberExpression*) e); }
     virtual void visit(LocalDeclaration *e)     { visit((Expression*) e); }
     virtual void visit(ArgumentExpression *e)   { visit((Expression*) e); }
     virtual void visit(PrototypeExpression *e)  { visit((Expression*) e); }
diff --git a/tests/modcc/driver.cpp b/tests/modcc/driver.cpp
index f32ac59d1471f94c6fb3d7d6333eeaf86d2e1c28..67505a0d20d3406d1ed8142715fb2243e74c6b0b 100644
--- a/tests/modcc/driver.cpp
+++ b/tests/modcc/driver.cpp
@@ -2,10 +2,17 @@
  * unit test driver
  **************************************************************/
 
+#include <cstring>
+
 #include "test.hpp"
 
+bool g_verbose_flag = false;
+
 int main(int argc, char **argv) {
     ::testing::InitGoogleTest(&argc, argv);
+    if (argc>1 && (!std::strcmp(argv[1],"-v") || !std::strcmp(argv[1],"--verbose"))) {
+        g_verbose_flag = true;
+    }
     return RUN_ALL_TESTS();
 }
 
diff --git a/tests/modcc/test.hpp b/tests/modcc/test.hpp
index 588b75a36afcb09107adc65b862fd65eceb5a9c6..d3de72506d2ef679227fcc71f0d76aab141b3a9c 100644
--- a/tests/modcc/test.hpp
+++ b/tests/modcc/test.hpp
@@ -5,12 +5,9 @@
 #include "parser.hpp"
 #include "modccutil.hpp"
 
-//#define VERBOSE_TEST
-#ifdef VERBOSE_TEST
-#define VERBOSE_PRINT(x) std::cout << (x) << std::endl;
-#else
-#define VERBOSE_PRINT(x)
-#endif
+extern bool g_verbose_flag;
+
+#define VERBOSE_PRINT(x) (g_verbose_flag && std::cout << (x) << "\n")
 
 inline expression_ptr parse_line_expression(std::string const& s) {
     return Parser(s).parse_line_expression();
diff --git a/tests/modcc/test_lexer.cpp b/tests/modcc/test_lexer.cpp
index b7d94515798b966b2e6d39f26ab68f1787f46e4b..e148b12b0631e1622bd52c2a4771f9b96802b23b 100644
--- a/tests/modcc/test_lexer.cpp
+++ b/tests/modcc/test_lexer.cpp
@@ -1,11 +1,55 @@
+#include <cctype>
 #include <cmath>
+#include <cstdio>
 #include <iterator>
+#include <utility>
 
 #include "test.hpp"
 #include "lexer.hpp"
 
-//#define PRINT_LEX_STRING std::cout << "________________\n" << string << "\n________________\n";
-#define PRINT_LEX_STRING
+void verbose_print(const char* string) {
+    if (!g_verbose_flag) return;
+    std::cout << "________________\n" << string << "\n________________\n";
+}
+
+void verbose_print(const Token& token) {
+    if (!g_verbose_flag) return;
+    std::cout << "tok: " << token << "\n";
+}
+
+class VerboseLexer: public Lexer {
+public:
+    template <typename... Args>
+    VerboseLexer(Args&&... args): Lexer(std::forward<Args>(args)...) {
+        if (g_verbose_flag) {
+            std::cout << "________________\n" << std::string(begin_, end_) << "\n________________\n";
+        }
+    }
+
+    Token parse() {
+        auto tok = Lexer::parse();
+        if (g_verbose_flag) {
+            std::cout << "token: " << tok << "\n";
+        }
+        return tok;
+    }
+
+    char character() {
+        char c = Lexer::character();
+        if (g_verbose_flag) {
+            std::cout << "character: ";
+            if (!std::isprint(c)) {
+                char buf[5] = "XXXX";
+                snprintf(buf, sizeof buf, "0x%02x", (unsigned)c);
+                std::cout << buf << '\n';
+            }
+            else {
+                std::cout << c << '\n';
+            }
+        }
+        return c;
+    }
+};
 
 /**************************************************************
  * lexer tests
@@ -13,8 +57,7 @@
 // test identifiers
 TEST(Lexer, identifiers) {
     char string[] = "_foo:\nbar, buzz f_zz";
-    PRINT_LEX_STRING
-    Lexer lexer(string, string+sizeof(string));
+    VerboseLexer lexer(string, string+sizeof(string));
 
     auto t1 = lexer.parse();
     EXPECT_EQ(t1.type, tok::identifier);
@@ -43,9 +86,8 @@ TEST(Lexer, identifiers) {
 
 // test keywords
 TEST(Lexer, keywords) {
-    char string[] = "NEURON UNITS SOLVE else TITLE CONDUCTANCE";
-    PRINT_LEX_STRING
-    Lexer lexer(string, string+sizeof(string));
+    char string[] = "NEURON UNITS SOLVE else TITLE CONDUCTANCE KINETIC CONSERVE LOCAL";
+    VerboseLexer lexer(string, string+sizeof(string));
 
     // should skip all white space and go straight to eof
     auto t1 = lexer.parse();
@@ -69,19 +111,30 @@ TEST(Lexer, keywords) {
     EXPECT_NE(t5.type, tok::identifier);
     EXPECT_EQ(t5.spelling, "TITLE");
 
+    auto t6 = lexer.parse();
+    EXPECT_EQ(t6.type, tok::conductance);
+    EXPECT_EQ(t6.spelling, "CONDUCTANCE");
+
     auto t7 = lexer.parse();
-    EXPECT_EQ(t7.type, tok::conductance);
-    EXPECT_EQ(t7.spelling, "CONDUCTANCE");
+    EXPECT_EQ(t7.type, tok::kinetic);
+    EXPECT_EQ(t7.spelling, "KINETIC");
 
-    auto t6 = lexer.parse();
-    EXPECT_EQ(t6.type, tok::eof);
+    auto t8 = lexer.parse();
+    EXPECT_EQ(t8.type, tok::conserve);
+    EXPECT_EQ(t8.spelling, "CONSERVE");
+
+    auto t9 = lexer.parse();
+    EXPECT_EQ(t9.type, tok::local);
+    EXPECT_EQ(t9.spelling, "LOCAL");
+
+    auto tlast = lexer.parse();
+    EXPECT_EQ(tlast.type, tok::eof);
 }
 
 // test white space
 TEST(Lexer, whitespace) {
     char string[] = " \t\v\f";
-    PRINT_LEX_STRING
-    Lexer lexer(string, string+sizeof(string));
+    VerboseLexer lexer(string, string+sizeof(string));
 
     // should skip all white space and go straight to eof
     auto t1 = lexer.parse();
@@ -91,8 +144,7 @@ TEST(Lexer, whitespace) {
 // test new line
 TEST(Lexer, newline) {
     char string[] = "foo \n    bar \n +\r\n-";
-    PRINT_LEX_STRING
-    Lexer lexer(string, string+sizeof(string));
+    VerboseLexer lexer(string, string+sizeof(string));
 
     // get foo
     auto t1 = lexer.parse();
@@ -123,9 +175,8 @@ TEST(Lexer, newline) {
 
 // test operators
 TEST(Lexer, symbols) {
-    char string[] = "+-/*, t= ^ h'";
-    PRINT_LEX_STRING
-    Lexer lexer(string, string+sizeof(string));
+    char string[] = "+-/*, t= ^ h'<->~";
+    VerboseLexer lexer(string, string+sizeof(string));
 
     auto t1 = lexer.parse();
     EXPECT_EQ(t1.type, tok::plus);
@@ -161,13 +212,18 @@ TEST(Lexer, symbols) {
     EXPECT_EQ(t10.type, tok::prime);
 
     auto t11 = lexer.parse();
-    EXPECT_EQ(t11.type, tok::eof);
+    EXPECT_EQ(t11.type, tok::arrow);
+
+    auto t12 = lexer.parse();
+    EXPECT_EQ(t12.type, tok::tilde);
+
+    auto tlast = lexer.parse();
+    EXPECT_EQ(tlast.type, tok::eof);
 }
 
 TEST(Lexer, comparison_operators) {
     char string[] = "< <= > >= == != !";
-    PRINT_LEX_STRING
-    Lexer lexer(string, string+sizeof(string));
+    VerboseLexer lexer(string, string+sizeof(string));
 
     auto t1 = lexer.parse();
     EXPECT_EQ(t1.type, tok::lt);
@@ -191,8 +247,7 @@ TEST(Lexer, comparison_operators) {
 // test braces
 TEST(Lexer, braces) {
     char string[] = "foo}";
-    PRINT_LEX_STRING
-    Lexer lexer(string, string+sizeof(string));
+    VerboseLexer lexer(string, string+sizeof(string));
 
     auto t1 = lexer.parse();
     EXPECT_EQ(t1.type, tok::identifier);
@@ -209,8 +264,7 @@ TEST(Lexer, comments) {
     char string[] = "foo:this is one line\n"
                     "bar : another comment\n"
                     "foobar ? another comment\n";
-    PRINT_LEX_STRING
-    Lexer lexer(string, string+sizeof(string));
+    VerboseLexer lexer(string, string+sizeof(string));
 
     auto t1 = lexer.parse();
     EXPECT_EQ(t1.type, tok::identifier);
@@ -231,6 +285,7 @@ TEST(Lexer, comments) {
 
 // test numbers
 TEST(Lexer, numbers) {
+    auto numeric = [](tok t) { return t==tok::real || t==tok::integer; };
     std::istringstream floats_stream("1 23 .3 87.99 12. 1.e3 1.2e+2 23e-3 -3");
 
     std::vector<double> floats;
@@ -238,7 +293,11 @@ TEST(Lexer, numbers) {
               std::istream_iterator<double>(),
               std::back_inserter(floats));
 
-    Lexer lexer(floats_stream.str());
+    // hand-parse these ...
+    std::vector<long long> check_ints = {1, 23, 3};
+    std::vector<long long> ints;
+
+    VerboseLexer lexer(floats_stream.str());
     auto t = lexer.parse();
     auto iter = floats.cbegin();
     while (t.type != tok::eof && iter != floats.cend()) {
@@ -249,11 +308,13 @@ TEST(Lexer, numbers) {
             // decide if the minus is a binary or unary expression
             EXPECT_EQ(tok::minus, t.type);
             t = lexer.parse();
-            EXPECT_EQ(tok::number, t.type);
+            EXPECT_TRUE(numeric(t.type));
+            if (t.type==tok::integer) ints.push_back(std::stoll(t.spelling));
             EXPECT_EQ(-(*iter), std::stod(t.spelling));
         }
         else {
-            EXPECT_EQ(t.type, tok::number);
+            EXPECT_TRUE(numeric(t.type));
+            if (t.type==tok::integer) ints.push_back(std::stoll(t.spelling));
             EXPECT_EQ(*iter, std::stod(t.spelling));
         }
 
@@ -263,4 +324,30 @@ TEST(Lexer, numbers) {
 
     EXPECT_EQ(floats.cend(), iter);
     EXPECT_EQ(tok::eof, t.type);
+    EXPECT_EQ(check_ints, ints);
+
+    // check case where 'E' is not followed by +, -, or a digit explicitly
+    lexer = VerboseLexer("7.2E");
+    t = lexer.parse();
+    EXPECT_EQ(lexerStatus::happy, lexer.status());
+    EXPECT_EQ(tok::real, t.type);
+    EXPECT_EQ(t.spelling, "7.2");
+    EXPECT_EQ(lexer.character(), 'E');
+
+    lexer = VerboseLexer("3E+E2");
+    t = lexer.parse();
+    EXPECT_EQ(lexerStatus::happy, lexer.status());
+    EXPECT_EQ(tok::integer, t.type);
+    EXPECT_EQ(t.spelling, "3");
+    EXPECT_EQ(lexer.character(), 'E');
+    EXPECT_EQ(lexer.character(), '+');
+
+    // 'bad' numbers should give errors
+    lexer = VerboseLexer("1.2.3");
+    lexer.parse();
+    EXPECT_EQ(lexerStatus::error, lexer.status());
+
+    lexer = VerboseLexer("1.2E4.3");
+    lexer.parse();
+    EXPECT_EQ(lexerStatus::error, lexer.status());
 }
diff --git a/tests/modcc/test_optimization.cpp b/tests/modcc/test_optimization.cpp
index 9b69cee7d2aa5c16af82e2a143cc9632c9b359b5..01f0a69c7e9cc8ec918a307bb77dd2a00d5e961b 100644
--- a/tests/modcc/test_optimization.cpp
+++ b/tests/modcc/test_optimization.cpp
@@ -9,23 +9,23 @@ TEST(Optimizer, constant_folding) {
     auto v = make_unique<ConstantFolderVisitor>();
     {
         auto e = parse_line_expression("x = 2*3");
-        VERBOSE_PRINT( e->to_string() )
+        VERBOSE_PRINT( e->to_string() );
         e->accept(v.get());
         EXPECT_EQ(e->is_assignment()->rhs()->is_number()->value(), 6);
-        VERBOSE_PRINT( e->to_string() )
-        VERBOSE_PRINT( "" )
+        VERBOSE_PRINT( e->to_string() );
+        VERBOSE_PRINT( "" );
     }
     {
         auto e = parse_line_expression("x = 1 + 2 + 3");
-        VERBOSE_PRINT( e->to_string() )
+        VERBOSE_PRINT( e->to_string() );
         e->accept(v.get());
         EXPECT_EQ(e->is_assignment()->rhs()->is_number()->value(), 6);
-        VERBOSE_PRINT( e->to_string() )
-        VERBOSE_PRINT( "" )
+        VERBOSE_PRINT( e->to_string() );
+        VERBOSE_PRINT( "" );
     }
     {
         auto e = parse_line_expression("x = exp(2)");
-        VERBOSE_PRINT( e->to_string() )
+        VERBOSE_PRINT( e->to_string() );
         e->accept(v.get());
         // The tolerance has to be loosend to 1e-15, because the optimizer performs
         // all intermediate calculations in 80 bit precision, which disagrees in
@@ -33,24 +33,24 @@ TEST(Optimizer, constant_folding) {
         // This is a good thing: by using the constant folder we increase accuracy
         // over the unoptimized code!
         EXPECT_EQ(std::fabs(e->is_assignment()->rhs()->is_number()->value()-std::exp(2.0))<1e-15, true);
-        VERBOSE_PRINT( e->to_string() )
-        VERBOSE_PRINT( "" )
+        VERBOSE_PRINT( e->to_string() );
+        VERBOSE_PRINT( "" );
     }
     {
         auto e = parse_line_expression("x= 2*2 + 3");
-        VERBOSE_PRINT( e->to_string() )
+        VERBOSE_PRINT( e->to_string() );
         e->accept(v.get());
         EXPECT_EQ(e->is_assignment()->rhs()->is_number()->value(), 7);
-        VERBOSE_PRINT( e->to_string() )
-        VERBOSE_PRINT( "" )
+        VERBOSE_PRINT( e->to_string() );
+        VERBOSE_PRINT( "" );
     }
     {
         auto e = parse_line_expression("x= 3 + 2*2");
-        VERBOSE_PRINT( e->to_string() )
+        VERBOSE_PRINT( e->to_string() );
         e->accept(v.get());
         EXPECT_EQ(e->is_assignment()->rhs()->is_number()->value(), 7);
-        VERBOSE_PRINT( e->to_string() )
-        VERBOSE_PRINT( "" )
+        VERBOSE_PRINT( e->to_string() );
+        VERBOSE_PRINT( "" );
     }
     {
         // this doesn't work: the (y+2) expression is not a constant, so folding stops.
@@ -58,23 +58,23 @@ TEST(Optimizer, constant_folding) {
         // one approach would be try sorting communtative operations so that numbers
         // are adjacent to one another in the tree
         auto e = parse_line_expression("x= y + 2 + 3");
-        VERBOSE_PRINT( e->to_string() )
+        VERBOSE_PRINT( e->to_string() );
         e->accept(v.get());
-        VERBOSE_PRINT( e->to_string() )
-        VERBOSE_PRINT( "" )
+        VERBOSE_PRINT( e->to_string() );
+        VERBOSE_PRINT( "" );
     }
     {
         auto e = parse_line_expression("x= 2 + 3 + y");
-        VERBOSE_PRINT( e->to_string() )
+        VERBOSE_PRINT( e->to_string() );
         e->accept(v.get());
-        VERBOSE_PRINT( e->to_string() )
-        VERBOSE_PRINT("");
+        VERBOSE_PRINT( e->to_string() );
+        VERBOSE_PRINT("");;
     }
     {
         auto e = parse_line_expression("foo(2+3, log(32), 2*3 + x)");
-        VERBOSE_PRINT( e->to_string() )
+        VERBOSE_PRINT( e->to_string() );
         e->accept(v.get());
-        VERBOSE_PRINT( e->to_string() )
-        VERBOSE_PRINT("");
+        VERBOSE_PRINT( e->to_string() );
+        VERBOSE_PRINT("");;
     }
 }
diff --git a/tests/modcc/test_parser.cpp b/tests/modcc/test_parser.cpp
index 1edfe1a9ad61a0e3837dbddba2877037fb30aee2..ab25d942dd5c6165b61bfd4940454ac840b9f993 100644
--- a/tests/modcc/test_parser.cpp
+++ b/tests/modcc/test_parser.cpp
@@ -1,12 +1,76 @@
 #include <cmath>
+#include <memory>
 
 #include "test.hpp"
 #include "module.hpp"
+#include "modccutil.hpp"
 #include "parser.hpp"
 
+template <typename EPtr>
+void verbose_print(const EPtr& e, Parser& p, const char* text) {
+    if (!g_verbose_flag) return;
+
+    if (e) std::cout << e->to_string() << "\n";
+    if (p.status()==lexerStatus::error)
+        std::cout << "in " << red(text) << "\t" << p.error_message() << "\n";
+}
+
+template <typename Derived, typename RetUniqPtr>
+::testing::AssertionResult check_parse(
+    std::unique_ptr<Derived>& derived,
+    RetUniqPtr (Parser::*pmemfn)(),
+    const char* text)
+{
+    Parser p(text);
+    auto e = (p.*pmemfn)();
+    verbose_print(e, p, text);
+
+    if (e==nullptr) {
+        return ::testing::AssertionFailure() << "failed to parse '" << text << "'";
+    }
+
+    if (p.status()!=lexerStatus::happy) {
+        return ::testing::AssertionFailure() << "parser status is not happy";
+    }
+
+    Derived *ptr = e? dynamic_cast<Derived*>(e.get()): nullptr;
+    if (ptr==nullptr) {
+        return ::testing::AssertionFailure() << "failed to cast to derived type";
+    }
+    else {
+        e.release();
+        derived.reset(ptr);
+    }
+
+    return ::testing::AssertionSuccess();
+}
+
+template <typename RetUniqPtr>
+::testing::AssertionResult check_parse(RetUniqPtr (Parser::*pmemfn)(), const char* text) {
+    std::unique_ptr<Expression> e;
+    return check_parse(e, pmemfn, text);
+}
+
+template <typename RetUniqPtr>
+::testing::AssertionResult check_parse_fail(RetUniqPtr (Parser::*pmemfn)(), const char* text) {
+    Parser p(text);
+    auto e = (p.*pmemfn)();
+    verbose_print(e, p, text);
+
+    if (p.status()!=lexerStatus::error) {
+        return ::testing::AssertionFailure() << "parser status is not error";
+    }
+
+    if (e!=nullptr) {
+        return ::testing::AssertionFailure() << "parser returned non-null expression";
+    }
+
+    return ::testing::AssertionSuccess();
+}
+
 TEST(Parser, full_file) {
     Module m(DATADIR "/test.mod");
-    if(m.buffer().size()==0) {
+    if (m.buffer().size()==0) {
         std::cout << "skipping Parser.full_file test because unable to open input file" << std::endl;
         return;
     }
@@ -15,481 +79,415 @@ TEST(Parser, full_file) {
 }
 
 TEST(Parser, procedure) {
-    std::vector<const char*> calls =
-{
-"PROCEDURE foo(x, y) {"
-"  LOCAL a\n"
-"  LOCAL b\n"
-"  LOCAL c\n"
-"  a = 3\n"
-"  b = x * y + 2\n"
-"  y = x + y * 2\n"
-"  y = a + b +c + a + b\n"
-"  y = a + b *c + a + b\n"
-"}"
-,
-"PROCEDURE trates(v) {\n"
-"    LOCAL qt\n"
-"    qt=q10^((celsius-22)/10)\n"
-"    minf=1-1/(1+exp((v-vhalfm)/km))\n"
-"    hinf=1/(1+exp((v-vhalfh)/kh))\n"
-"    mtau = 0.6\n"
-"    htau = 1500\n"
-"}"
-};
-    for(auto const& str : calls) {
-        Parser p(str);
-        auto e = p.parse_procedure();
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
-        if(p.status()==lexerStatus::error) {
-            std::cout << str << std::endl;
-            std::cout << red("error ") << p.error_message() << std::endl;
-        }
+    std::vector<const char*> calls = {
+        "PROCEDURE foo(x, y) {\n"
+        "  LOCAL a\n"
+        "  LOCAL b\n"
+        "  LOCAL c\n"
+        "  a = 3\n"
+        "  b = x * y + 2\n"
+        "  y = x + y * 2\n"
+        "  y = a + b +c + a + b\n"
+        "  y = a + b *c + a + b\n"
+        "}"
+        ,
+        "PROCEDURE trates(v) {\n"
+        "    LOCAL qt\n"
+        "    qt=q10^((celsius-22)/10)\n"
+        "    minf=1-1/(1+exp((v-vhalfm)/km))\n"
+        "    hinf=1/(1+exp((v-vhalfh)/kh))\n"
+        "    mtau = 0.6\n"
+        "    htau = 1500\n"
+        "}"
+    };
+
+    for (const auto& str: calls) {
+        EXPECT_TRUE(check_parse(&Parser::parse_procedure, str));
     }
 }
 
 TEST(Parser, net_receive) {
     char str[] =
-    "NET_RECEIVE (x, y) {   \n"
-    "  LOCAL a              \n"
-    "  a = 3                \n"
-    "  x = a+3              \n"
-    "  y = x+a              \n"
-    "}";
-    Parser p(str);
-    auto e = p.parse_procedure();
-    #ifdef VERBOSE_TEST
-    if(e) std::cout << e->to_string() << std::endl;
-    #endif
-
-    EXPECT_NE(e, nullptr);
-    EXPECT_EQ(p.status(), lexerStatus::happy);
-
-    auto nr = e->is_symbol()->is_net_receive();
-    EXPECT_NE(nr, nullptr);
-    if(nr) {
-        EXPECT_EQ(nr->args().size(), (unsigned)2);
-    }
-    if(p.status()==lexerStatus::error) {
-        std::cout << str << std::endl;
-        std::cout << red("error ") << p.error_message() << std::endl;
+        "NET_RECEIVE (x, y) {   \n"
+        "  LOCAL a              \n"
+        "  a = 3                \n"
+        "  x = a+3              \n"
+        "  y = x+a              \n"
+        "}";
+
+    std::unique_ptr<Symbol> sym;
+
+    EXPECT_TRUE(check_parse(sym, &Parser::parse_procedure, str));
+    if (sym) {
+        auto nr = sym->is_net_receive();
+        EXPECT_NE(nullptr, nr);
+        if (nr) {
+            EXPECT_EQ(2u, nr->args().size());
+        }
     }
 }
 
 TEST(Parser, function) {
-    std::vector< const char*> calls =
-{
-"FUNCTION foo(x, y) {"
-"  LOCAL a\n"
-"  a = 3\n"
-"  b = x * y + 2\n"
-"  y = x + y * 2\n"
-"  foo = a * x + y\n"
-"}"
-};
-    for(auto const& str : calls) {
-        Parser p(str);
-        auto e = p.parse_function();
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
-        if(p.status()==lexerStatus::error) {
-            std::cout << str << std::endl;
-            std::cout << red("error ") << p.error_message() << std::endl;
-        }
-    }
+    char str[] =
+        "FUNCTION foo(x, y) {"
+        "  LOCAL a\n"
+        "  a = 3\n"
+        "  b = x * y + 2\n"
+        "  y = x + y * 2\n"
+        "  foo = a * x + y\n"
+        "}";
+
+    std::unique_ptr<Symbol> sym;
+    EXPECT_TRUE(check_parse(sym, &Parser::parse_function, str));
 }
 
 TEST(Parser, parse_solve) {
-    {
-        Parser p("SOLVE states METHOD cnexp");
-        auto e = p.parse_solve();
-
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
-
-        if(e) {
-            SolveExpression* s = dynamic_cast<SolveExpression*>(e.get());
-            EXPECT_EQ(s->method(), solverMethod::cnexp);
-            EXPECT_EQ(s->name(), "states");
-        }
+    std::unique_ptr<SolveExpression> s;
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error) {
-            std::cout << red("error") << p.error_message() << std::endl;
-        }
+    EXPECT_TRUE(check_parse(s, &Parser::parse_solve, "SOLVE states METHOD cnexp"));
+    if (s) {
+        EXPECT_EQ(s->method(), solverMethod::cnexp);
+        EXPECT_EQ(s->name(), "states");
     }
-    {
-        Parser p("SOLVE states");
-        auto e = p.parse_solve();
-
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
-
-        if(e) {
-            SolveExpression* s = dynamic_cast<SolveExpression*>(e.get());
-            EXPECT_EQ(s->method(), solverMethod::none);
-            EXPECT_EQ(s->name(), "states");
-        }
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error) {
-            std::cout << red("error") << p.error_message() << std::endl;
-        }
+    EXPECT_TRUE(check_parse(s, &Parser::parse_solve, "SOLVE states"));
+    if (s) {
+        EXPECT_EQ(s->method(), solverMethod::none);
+        EXPECT_EQ(s->name(), "states");
     }
 }
 
 TEST(Parser, parse_conductance) {
-    {
-        Parser p("CONDUCTANCE g USEION na");
-        auto e = p.parse_conductance();
-
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
-
-        if(e) {
-            ConductanceExpression* s = dynamic_cast<ConductanceExpression*>(e.get());
-            EXPECT_EQ(s->ion_channel(), ionKind::Na);
-            EXPECT_EQ(s->name(), "g");
-        }
+    std::unique_ptr<ConductanceExpression> s;
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error) {
-            std::cout << red("error") << p.error_message() << std::endl;
-        }
+    EXPECT_TRUE(check_parse(s, &Parser::parse_conductance, "CONDUCTANCE g USEION na"));
+    if (s) {
+        EXPECT_EQ(s->ion_channel(), ionKind::Na);
+        EXPECT_EQ(s->name(), "g");
     }
-    {
-        Parser p("CONDUCTANCE gnda");
-        auto e = p.parse_conductance();
-
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
-
-        if(e) {
-            ConductanceExpression* s = dynamic_cast<ConductanceExpression*>(e.get());
-            EXPECT_EQ(s->ion_channel(), ionKind::nonspecific);
-            EXPECT_EQ(s->name(), "gnda");
-        }
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error) {
-            std::cout << red("error") << p.error_message() << std::endl;
-        }
+    EXPECT_TRUE(check_parse(s, &Parser::parse_conductance, "CONDUCTANCE gnda"));
+    if (s) {
+        EXPECT_EQ(s->ion_channel(), ionKind::nonspecific);
+        EXPECT_EQ(s->name(), "gnda");
     }
 }
 
 TEST(Parser, parse_if) {
-    {
-        char expression[] =
+    std::unique_ptr<IfExpression> s;
+
+    EXPECT_TRUE(check_parse(s, &Parser::parse_if,
         "   if(a<b) {      \n"
         "       a = 2+b    \n"
         "       b = 4^b    \n"
-        "   }              \n";
-        Parser p(expression);
-        auto e = p.parse_if();
-        EXPECT_NE(e, nullptr);
-        if(e) {
-            auto ife = e->is_if();
-            EXPECT_NE(e->is_if(), nullptr);
-            if(ife) {
-                EXPECT_NE(ife->condition()->is_binary(), nullptr);
-                EXPECT_NE(ife->true_branch()->is_block(), nullptr);
-                EXPECT_EQ(ife->false_branch(), nullptr);
-            }
-            //std::cout << e->to_string() << std::endl;
-        }
-        else {
-            std::cout << p.error_message() << std::endl;
-        }
+        "   }              \n"
+    ));
+    if (s) {
+        EXPECT_NE(s->condition()->is_binary(), nullptr);
+        EXPECT_NE(s->true_branch()->is_block(), nullptr);
+        EXPECT_EQ(s->false_branch(), nullptr);
     }
-    {
-        char expression[] =
+
+    EXPECT_TRUE(check_parse(s, &Parser::parse_if,
         "   if(a<b) {      \n"
         "       a = 2+b    \n"
         "   } else {       \n"
         "       a = 2+b    \n"
-        "   }                ";
-        Parser p(expression);
-        auto e = p.parse_if();
-        EXPECT_NE(e, nullptr);
-        if(e) {
-            auto ife = e->is_if();
-            EXPECT_NE(ife, nullptr);
-            if(ife) {
-                EXPECT_NE(ife->condition()->is_binary(), nullptr);
-                EXPECT_NE(ife->true_branch()->is_block(), nullptr);
-                EXPECT_NE(ife->false_branch(), nullptr);
-            }
-            //std::cout << std::endl << e->to_string() << std::endl;
-        }
-        else {
-            std::cout << p.error_message() << std::endl;
-        }
+        "   }                "
+    ));
+    if (s) {
+        EXPECT_NE(s->condition()->is_binary(), nullptr);
+        EXPECT_NE(s->true_branch()->is_block(), nullptr);
+        EXPECT_NE(s->false_branch(), nullptr);
     }
-    {
-        char expression[] =
+
+    EXPECT_TRUE(check_parse(s, &Parser::parse_if,
         "   if(a<b) {      \n"
         "       a = 2+b    \n"
         "   } else if(b>a){\n"
         "       a = 2+b    \n"
-        "   }              ";
-        Parser p(expression);
-        auto e = p.parse_if();
-        EXPECT_NE(e, nullptr);
-        if(e) {
-            auto ife = e->is_if();
-            EXPECT_NE(ife, nullptr);
-            if(ife) {
-                EXPECT_NE(ife->condition()->is_binary(), nullptr);
-                EXPECT_NE(ife->true_branch()->is_block(), nullptr);
-                EXPECT_NE(ife->false_branch(), nullptr);
-                EXPECT_NE(ife->false_branch()->is_if(), nullptr);
-                EXPECT_EQ(ife->false_branch()->is_if()->false_branch(), nullptr);
-            }
-            //std::cout << std::endl << e->to_string() << std::endl;
-        }
-        else {
-            std::cout << p.error_message() << std::endl;
-        }
+        "   }              "
+    ));
+    if (s) {
+        EXPECT_NE(s->condition()->is_binary(), nullptr);
+        EXPECT_NE(s->true_branch()->is_block(), nullptr);
+        ASSERT_NE(s->false_branch(), nullptr);
+        ASSERT_NE(s->false_branch()->is_if(), nullptr);
+        EXPECT_EQ(s->false_branch()->is_if()->false_branch(), nullptr);
     }
 }
 
 TEST(Parser, parse_local) {
-    ////////////////////// test for valid expressions //////////////////////
-    {
-        Parser p("LOCAL xyz");
-        auto e = p.parse_local();
-
-        #ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-        #endif
-        EXPECT_NE(e, nullptr);
-        if(e) {
-            EXPECT_NE(e->is_local_declaration(), nullptr);
-            EXPECT_EQ(p.status(), lexerStatus::happy);
-        }
+    std::unique_ptr<LocalDeclaration> s;
+    EXPECT_TRUE(check_parse(s, &Parser::parse_local, "LOCAL xyz"));
+    if (s) {
+        ASSERT_EQ(1u, s->variables().size());
+    }
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error)
-            std::cout << red("error") << p.error_message() << std::endl;
+    EXPECT_TRUE(check_parse(s, &Parser::parse_local, "LOCAL x, y, z"));
+    if (s) {
+        auto vars = s->variables();
+        ASSERT_EQ(3u, vars.size());
+        ASSERT_TRUE(vars.count("x"));
+        ASSERT_TRUE(vars.count("y"));
+        ASSERT_TRUE(vars.count("z"));
     }
 
-    {
-        Parser p("LOCAL x, y, z");
-        auto e = p.parse_local();
-
-        #ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-        #endif
-        EXPECT_NE(e, nullptr);
-        if(e) {
-            EXPECT_NE(e->is_local_declaration(), nullptr);
-            EXPECT_EQ(p.status(), lexerStatus::happy);
-            auto vars = e->is_local_declaration()->variables();
-            EXPECT_EQ(vars.size(), (unsigned)3);
-            EXPECT_NE(vars.find("x"), vars.end());
-            EXPECT_NE(vars.find("y"), vars.end());
-            EXPECT_NE(vars.find("z"), vars.end());
-        }
+    EXPECT_TRUE(check_parse_fail(&Parser::parse_local, "LOCAL x,"));
+}
+
+TEST(Parser, parse_unary_expression) {
+    const char* good_expr[] = {
+        "+x             ",
+        "-x             ",
+        "(x + -y)       ",
+        "-(x - + -y)    ",
+        "exp(x + y)     ",
+        "-exp(x + -y)   "
+    };
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error)
-            std::cout << red("error") << p.error_message() << std::endl;
+    for (auto& text: good_expr) {
+        EXPECT_TRUE(check_parse(&Parser::parse_unaryop, text));
     }
+}
 
-    ////////////////////// test for invalid expressions //////////////////////
-    {
-        Parser p("LOCAL 2");
-        auto e = p.parse_local();
+// test parsing of parenthesis expressions
+TEST(Parser, parse_parenthesis_expression) {
+    const char* good_expr[] = {
+        "((celsius-22)/10)      ",
+        "((celsius-22)+10)      ",
+        "(x+2)                  ",
+        "((x))                  ",
+        "(((x)))                ",
+        "(x + (x * (y*(2)) + 4))",
+    };
 
-        EXPECT_EQ(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::error);
+    for (auto& text: good_expr) {
+        EXPECT_TRUE(check_parse(&Parser::parse_parenthesis_expression, text));
+    }
 
-        #ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-        if(p.status()==lexerStatus::error)
-            std::cout << "in " << cyan(bad_expression) << "\t" << p.error_message() << std::endl;
-        #endif
+    const char* bad_expr[] = {
+        "(x             ",
+        "((x+3)         ",
+        "(x+ +)         ",
+        "(x=3)          ",  // assignment inside parenthesis isn't allowed
+        "(a + (b*2^(x)) ",  // missing closing parenthesis
+    };
+
+    for (auto& text: bad_expr) {
+        EXPECT_TRUE(check_parse_fail(&Parser::parse_parenthesis_expression, text));
     }
+}
 
-    {
-        Parser p("LOCAL x, ");
-        auto e = p.parse_local();
+// test parsing of line expressions
+TEST(Parser, parse_line_expression) {
+    const char* good_expr[] = {
+        "qt=q10^((celsius-22)/10)",
+        "x=2        ",
+        "x=2        ",
+        "x = -y\n   "
+        "x=2*y      ",
+        "x=y + 2 * z",
+        "x=(y + 2) * z      ",
+        "x=(y + 2) * z ^ 3  ",
+        "x=(y + 2 * z ^ 3)  ",
+        "foo(x+3, y, bar(21.4))",
+        "y=exp(x+3) + log(exp(x/y))",
+        "a=x^y^z",
+        "a=x/y/z"
+    };
 
-        EXPECT_EQ(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::error);
+    for (auto& text: good_expr) {
+        EXPECT_TRUE(check_parse(&Parser::parse_line_expression, text));
+    }
 
-        #ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-        if(p.status()==lexerStatus::error)
-            std::cout << "in " << cyan(bad_expression) << "\t" << p.error_message() << std::endl;
-        #endif
+    const char* bad_expr[] = {
+        "x=2+        ",      // incomplete binary expression on rhs
+        "x=          ",      // missing rhs of assignment
+        "x=)y + 2 * z",
+        "x=(y + 2    ",
+        "x=(y ++ z   ",
+        "x/=3        ",      // compound binary expressions not supported
+        "foo+8       ",      // missing assignment
+        "foo()=8     ",      // lhs of assingment must be an lvalue
+    };
+
+    for (auto& text: bad_expr) {
+        EXPECT_TRUE(check_parse_fail(&Parser::parse_line_expression, text));
     }
 }
 
-TEST(Parser, parse_unary_expression) {
-    std::vector<const char*> good_expressions =
-    {
-"+x             ",
-"-x             ",
-"(x + -y)       ",
-"-(x - + -y)    ",
-"exp(x + y)     ",
-"-exp(x + -y)   ",
+TEST(Parser, parse_stoich_term) {
+    const char* good_pos_expr[] = {
+        "B", "B3", "3B3", "0A", "12A", "4E"
     };
 
-    for(auto const& expression : good_expressions) {
-        Parser p(expression);
-        auto e = p.parse_unaryop();
+    for (auto& text: good_pos_expr) {
+        std::unique_ptr<StoichTermExpression> s;
+        EXPECT_TRUE(check_parse(s, &Parser::parse_stoich_term, text));
+        EXPECT_TRUE((s && !s->negative()));
+    }
 
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
+    const char* good_neg_expr[] = {
+        "-3B3", "-A", "-12A"
+    };
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error)
-            std::cout << red("error") << p.error_message() << std::endl;
+    for (auto& text: good_neg_expr) {
+        std::unique_ptr<StoichTermExpression> s;
+        EXPECT_TRUE(check_parse(s, &Parser::parse_stoich_term, text));
+        EXPECT_TRUE((s && s->negative()));
+    }
+    const char* bad_expr[] = {
+        "0.2A", "5", "3e2" // "3e2" should lex as real number 300.0
+    };
+
+    for (auto& text: bad_expr) {
+        EXPECT_TRUE(check_parse_fail(&Parser::parse_stoich_term, text));
     }
 }
 
-// test parsing of parenthesis expressions
-TEST(Parser, parse_parenthesis_expression) {
-    std::vector<const char*> good_expressions =
-    {
-"((celsius-22)/10)      ",
-"((celsius-22)+10)      ",
-"(x+2)                  ",
-"((x))                  ",
-"(((x)))                ",
-"(x + (x * (y*(2)) + 4))",
+TEST(Parser, parse_stoich_expression) {
+    const char* single_expr[] = {
+        "B", "B3", "3xy"
     };
 
-    for(auto const& expression : good_expressions) {
-        Parser p(expression);
-        auto e = p.parse_parenthesis_expression();
+    for (auto& text: single_expr) {
+        std::unique_ptr<StoichExpression> s;
+        EXPECT_TRUE(check_parse(s, &Parser::parse_stoich_expression, text));
+        EXPECT_EQ(1, s->terms().size());
+    }
 
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
+    const char* double_expr[] = {
+        "B+A", "a1 + 2bn", "4c+d"
+    };
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error)
-            std::cout << cyan(expression) << "\t"
-                      << red("error") << p.error_message() << std::endl;
+    for (auto& text: double_expr) {
+        std::unique_ptr<StoichExpression> s;
+        EXPECT_TRUE(check_parse(s, &Parser::parse_stoich_expression, text));
+        EXPECT_EQ(2, s->terms().size());
     }
 
-    std::vector<const char*> bad_expressions =
-    {
-"(x             ",
-"((x+3)         ",
-"(x+ +)         ",
-"(x=3)          ",  // assignment inside parenthesis isn't allowed
-"(a + (b*2^(x)) ",  // missing closing parenthesis
+    const char* other_good_expr[] = {
+        "", "a+b+c", "1a-2b+3c+4d"
     };
 
-    for(auto const& expression : bad_expressions) {
-        Parser p(expression);
-        auto e = p.parse_parenthesis_expression();
+    for (auto& text: other_good_expr) {
+        std::unique_ptr<StoichExpression> s;
+        EXPECT_TRUE(check_parse(s, &Parser::parse_stoich_expression, text));
+    }
 
-        EXPECT_EQ(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::error);
+    const char* check_coeff = "-3a+2b-c+d";
+    {
+        std::unique_ptr<StoichExpression> s;
+        EXPECT_TRUE(check_parse(s, &Parser::parse_stoich_expression, check_coeff));
+        EXPECT_EQ(4, s->terms().size());
+        std::vector<int> confirm = {-3,2,-1,1};
+        for (unsigned i = 0; i<4; ++i) {
+            auto term = s->terms()[i]->is_stoich_term();
+            EXPECT_EQ(confirm[i], term->coeff()->is_integer()->integer_value());
+        }
+    }
 
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-        if(p.status()==lexerStatus::error)
-            std::cout << "in " << cyan(expression) << "\t" << p.error_message() << std::endl;
-#endif
+    const char* bad_expr[] = {
+        "A+B+", "A+5+B"
+    };
+
+    for (auto& text: bad_expr) {
+        EXPECT_TRUE(check_parse_fail(&Parser::parse_stoich_expression, text));
     }
 }
 
-// test parsing of line expressions
-TEST(Parser, parse_line_expression) {
-    std::vector<const char*> good_expressions =
-    {
-"qt=q10^((celsius-22)/10)"
-"x=2        ",
-"x=2        ",
-"x = -y\n   "
-"x=2*y      ",
-"x=y + 2 * z",
-"x=(y + 2) * z      ",
-"x=(y + 2) * z ^ 3  ",
-"x=(y + 2 * z ^ 3)  ",
-"foo(x+3, y, bar(21.4))",
-"y=exp(x+3) + log(exp(x/y))",
-"a=x^y^z",
-"a=x/y/z"
+// test parsing of stoich and reaction expressions
+TEST(Parser, parse_reaction_expression) {
+    const char* good_expr[] = {
+        "~ A + B <-> C + D (k1, k2)",
+        "~ 2B <-> C + D + E (k1(3,v), k2)",
+        "~ <-> C + D + 7 E (k1, f(a,b)-2)",
+        "~ <-> C + D + 7E+F (k1, f(a,b)-2)",
+        "~ <-> (f,g)",
+        "~ A + 3B + C<-> (f,g)"
     };
 
-    for(auto const& expression : good_expressions) {
-        Parser p(expression);
-        auto e = p.parse_line_expression();
+    for (auto& text: good_expr) {
+        std::unique_ptr<ReactionExpression> s;
+        EXPECT_TRUE(check_parse(s, &Parser::parse_reaction_expression, text));
+    }
 
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
+    const char* bad_expr[] = {
+        "~ A + B <-> C + D (k1, k2, k3)",
+        "~ A + B <-> C + (k1, k2)",
+        "~ 2.3B <-> C + D + E (k1(3,v), k2)",
+        "~ <-> C + D + 7E",
+        "~ <-> C + D + 7E+2F (k1, f(a,b)-2)", // "7E+2" will lex as real number
+        "~ <-> (,g)",
+        "~ A - 3B + C<-> (f,g)",
+        "  A <-> B (k1, k2)",
+        "~ A <- B (k1)",
+        "~ A -> B (k2)",
+    };
 
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error)
-            std::cout << red("error") << p.error_message() << std::endl;
+    for (auto& text: bad_expr) {
+        EXPECT_TRUE(check_parse_fail(&Parser::parse_reaction_expression, text));
     }
+}
 
-    std::vector<const char*> bad_expressions =
+TEST(Parser, parse_conserve) {
+    std::unique_ptr<ConserveExpression> s;
+    const char* text;
+
+    text = "CONSERVE a + b = 1";
+    ASSERT_TRUE(check_parse(s, &Parser::parse_conserve_expression, text));
+    EXPECT_TRUE(s->rhs()->is_number());
+    ASSERT_TRUE(s->lhs()->is_stoich());
+    EXPECT_EQ(2, s->lhs()->is_stoich()->terms().size());
+
+    text = "CONSERVE a = 1.23e-2";
+    ASSERT_TRUE(check_parse(s, &Parser::parse_conserve_expression, text));
+    EXPECT_TRUE(s->rhs()->is_number());
+    ASSERT_TRUE(s->lhs()->is_stoich());
+    EXPECT_EQ(1, s->lhs()->is_stoich()->terms().size());
+
+    text = "CONSERVE = 0";
+    ASSERT_TRUE(check_parse(s, &Parser::parse_conserve_expression, text));
+    EXPECT_TRUE(s->rhs()->is_number());
+    ASSERT_TRUE(s->lhs()->is_stoich());
+    EXPECT_EQ(0, s->lhs()->is_stoich()->terms().size());
+
+    text = "CONSERVE -2a + b -c = foo*2.3-bar";
+    ASSERT_TRUE(check_parse(s, &Parser::parse_conserve_expression, text));
+    EXPECT_TRUE(s->rhs()->is_binary());
+    ASSERT_TRUE(s->lhs()->is_stoich());
     {
-"x=2+       ",      // incomplete binary expression on rhs
-"x=         ",      // missing rhs of assignment
-"x=)y + 2 * z",
-"x=(y + 2   ",
-"x=(y ++ z  ",
-"x/=3       ",      // compound binary expressions not supported
-"foo+8      ",      // missing assignment
-"foo()=8    ",      // lhs of assingment must be an lvalue
-    };
-
-    for(auto const& expression : bad_expressions) {
-        Parser p(expression);
-        auto e = p.parse_line_expression();
+        auto& terms = s->lhs()->is_stoich()->terms();
+        ASSERT_EQ(3, terms.size());
+        auto coeff = terms[0]->is_stoich_term()->coeff()->is_integer();
+        ASSERT_TRUE(coeff);
+        EXPECT_EQ(-2, coeff->integer_value());
+        coeff = terms[1]->is_stoich_term()->coeff()->is_integer();
+        ASSERT_TRUE(coeff);
+        EXPECT_EQ(1, coeff->integer_value());
+        coeff = terms[2]->is_stoich_term()->coeff()->is_integer();
+        ASSERT_TRUE(coeff);
+        EXPECT_EQ(-1, coeff->integer_value());
+    }
 
-        EXPECT_EQ(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::error);
+    const char* bad_expr[] = {
+        "CONSERVE a + 3*b -c = 1",
+        "CONSERVE a + 3b -c = ",
+        "a+b+c = 2",
+        "CONSERVE a + 3b +c"
+    };
 
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-        if(p.status()==lexerStatus::error)
-            std::cout << "in " << cyan(expression) << "\t" << p.error_message() << std::endl;
-#endif
+    for (auto& text: bad_expr) {
+        EXPECT_TRUE(check_parse_fail(&Parser::parse_conserve_expression, text));
     }
 }
 
 long double eval(Expression *e) {
-    if(auto n = e->is_number()) {
+    if (auto n = e->is_number()) {
         return n->value();
     }
-    if(auto b = e->is_binary()) {
+    if (auto b = e->is_binary()) {
         auto lhs = eval(b->lhs());
         auto rhs = eval(b->rhs());
         switch(b->op()) {
@@ -501,7 +499,7 @@ long double eval(Expression *e) {
             default:;
         }
     }
-    if(auto u = e->is_unary()) {
+    if (auto u = e->is_unary()) {
         auto val = eval(u->expression());
         switch(u->op()) {
             case tok::plus  : return  val;
@@ -515,8 +513,7 @@ long double eval(Expression *e) {
 // test parsing of expressions for correctness
 // by parsing rvalue expressions with numeric atoms, which can be evalutated using eval
 TEST(Parser, parse_binop) {
-    std::vector<std::pair<const char*, double>> tests =
-    {
+    std::pair<const char*, double> tests[] = {
         // simple
         {"2+3", 2.+3.},
         {"2-3", 2.-3.},
@@ -540,15 +537,9 @@ TEST(Parser, parse_binop) {
         {"3^2*5.", std::pow(3.,2.)*5.},
     };
 
-    for(auto const& test_case : tests) {
-        Parser p(test_case.first);
-        auto e = p.parse_expression();
-
-#ifdef VERBOSE_TEST
-        if(e) std::cout << e->to_string() << std::endl;
-#endif
-        EXPECT_NE(e, nullptr);
-        EXPECT_EQ(p.status(), lexerStatus::happy);
+    for (const auto& test_case: tests) {
+        std::unique_ptr<Expression> e;
+        EXPECT_TRUE(check_parse(e, &Parser::parse_expression, test_case.first));
 
         // A loose tolerance of 1d-10 is required here because the eval()
         // function uses long double for intermediate results (like constant
@@ -556,15 +547,11 @@ TEST(Parser, parse_binop) {
         // operations this can see relatively large divergence between the
         // double and long double results.
         EXPECT_NEAR(eval(e.get()), test_case.second, 1e-10);
-
-        // always print the compiler errors, because they are unexpected
-        if(p.status()==lexerStatus::error)
-            std::cout << red("error") << p.error_message() << std::endl;
     }
 }
 
 TEST(Parser, parse_state_block) {
-    std::vector<const char*> state_blocks = {
+    const char* state_blocks[] = {
         "STATE {\n"
         "    h\n"
         "    m r\n"
@@ -587,14 +574,26 @@ TEST(Parser, parse_state_block) {
         "}"
     };
 
-    for (auto const& str: state_blocks) {
-        Module m(str, sizeof(str));
+    expression_ptr null;
+    for (auto& text: state_blocks) {
+        Module m(text, sizeof(text));
         Parser p(m, false);
         p.parse_state_block();
         EXPECT_EQ(lexerStatus::happy, p.status());
-        if (p.status() == lexerStatus::error) {
-            std::cout << str << "\n"
-                      << red("error") << p.error_message() << "\n";
-        }
+        verbose_print(null, p, text);
     }
 }
+
+TEST(Parser, parse_kinetic) {
+    char str[] =
+        "KINETIC kin {\n"
+        "    rates(v)             \n"
+        "    ~ s1 <-> s2 (f1, r1) \n"
+        "    ~ s2 <-> s3 (f2, r2) \n"
+        "    ~ s2 <-> s4 (f3, r3) \n"
+        "    CONSERVE s1 + s3 + s4 - s2 = 2.3\n"
+        "}";
+
+    std::unique_ptr<Symbol> sym;
+    EXPECT_TRUE(check_parse(sym, &Parser::parse_procedure, str));
+}