diff --git a/Project.toml b/Project.toml
index 5216df5c..60761477 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,8 @@
 name = "PEtab"
 uuid = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
 authors = ["Viktor Hasselgren", "Sebastian Persson", "Rafael Arutjunjan", "Damiano Ognissanti"]
-version = "2.0.1"
+version = "2.0.2"
+
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -52,10 +53,10 @@ PEtabCatalystExtension = ["Catalyst"]
 PEtabFidesExtension = ["PyCall"]
 PEtabIpoptExtension = ["Ipopt"]
 PEtabOptimExtension = ["Optim"]
+PEtabPlotsExtension = ["Plots"]
 PEtabSciMLSensitivityExtension = ["SciMLSensitivity", "Zygote"]
 PEtabSelectExtension = ["PyCall"]
 PEtabZygoteExtension = ["SciMLSensitivity", "Zygote"]
-PEtabPlotsExtension = ["Plots"]
 
 [compat]
 CSV = "0.10"
diff --git a/src/SBML/Common.jl b/src/SBML/Common.jl
index 7f483736..24615099 100644
--- a/src/SBML/Common.jl
+++ b/src/SBML/Common.jl
@@ -51,10 +51,17 @@ function rewrite_piecewise_to_ifelse(rule_formula, variable, model_dict, base_fu
 
         if condition[1:2] == "lt" || condition[1:2] == "gt" || condition[1:2] == "eq" || condition[1:3] == "neq" || condition[1:3] == "geq" || condition[1:3] == "leq" 
             eq_syntax_dict[variable_change] = simple_piecewise_to_ifelse(condition, variable_change, value_active, value_inactive, model_dict, base_functions)
+
         elseif condition[1:3] == "and" || condition[1:2] == "if" || condition[1:2] == "or" || condition[1:3] == "xor" || condition[1:3] == "not"
             eq_syntax_dict[variable_change] = complex_piecewise_to_ifelse(condition, variable, value_active, value_inactive, model_dict, base_functions)
+
+        # Recursion to handle nested piecewise in condition             
+        elseif length(condition) ≥ 9 && condition[1:9] == "piecewise"
+            condition = rewrite_piecewise_to_ifelse(condition, "foo", model_dict, base_functions, model_SBML, ret_formula=true)
+            condition = rewrite_derivatives(condition, model_dict, base_functions, model_SBML)
+            eq_syntax_dict[variable_change] = simple_piecewise_to_ifelse(condition, variable_change, value_active, value_inactive, model_dict, base_functions)
         else
-            @error "Somehow we cannot process the piecewise expression"
+            @error "Somehow we cannot process the piecewise expression, condition = $condition"
         end
     end
 
@@ -124,6 +131,7 @@ end
 
 function simple_piecewise_to_ifelse(condition, variable, value_active, value_inactive, dicts, base_functions)
 
+    nested_condition::Bool = false
     if "leq" == condition[1:3]
         stripped_condition = condition[5:end-1]
         inEqUse = " <= "
@@ -146,13 +154,18 @@ function simple_piecewise_to_ifelse(condition, variable, value_active, value_ina
         return "true"
     elseif "false" == condition[1:5]
         return "false"
+    elseif "ifelse" == condition[1:6]
+        nested_condition = true
     else
         @error "Cannot recognize form of inequality, condition = $condition"
     end
 
-    parts = split_between(stripped_condition, ',')
-    # Trigger of event
-    expression = "ifelse(" * parts[1] * inEqUse * parts[2] * ", " * value_active * ", " * value_inactive * ")"
+    if nested_condition == false
+        parts = split_between(stripped_condition, ',')
+        expression = "ifelse(" * parts[1] * inEqUse * parts[2] * ", " * value_active * ", " * value_inactive * ")"
+    else
+        expression = "ifelse(" * condition * ", " * value_active * ", " * value_inactive * ")"
+    end
 
     return expression
 end
diff --git a/src/SBML/Process_rules.jl b/src/SBML/Process_rules.jl
index efca55a4..100ca286 100644
--- a/src/SBML/Process_rules.jl
+++ b/src/SBML/Process_rules.jl
@@ -18,6 +18,8 @@ function process_assignment_rule!(model_dict::Dict, rule_formula::String, rule_v
         return 
     end
 
+    rule_formula = replace_reactionid_with_math(rule_formula, model_SBML)
+
     #=
         If the rule does not involve a piecewise expression simply encode it as a function which downsteram
         is integrated into the equations when inserting "functions" into the reactions 
@@ -69,6 +71,8 @@ function process_rate_rule!(model_dict::Dict, rule_formula::String, rule_variabl
         end
     end
 
+    rule_formula = replace_reactionid_with_math(rule_formula, model_SBML)
+
     # Add rate rule as part of model derivatives and remove from parameters dict if rate rule variable
     # is a parameter
     if rule_variable ∈ keys(model_dict["parameters"])
@@ -153,7 +157,7 @@ function _time_dependent_ifelse_to_bool(formula_with_ifelse::String, model_dict:
         time_left = check_for_time(string(lhsRule))
         rewrite_ifelse = true
         if time_left == false && time_left == false
-            println("Have ifelse statements which does not contain time. Hence we do not rewrite as event, but rather keep it as an ifelse.")
+            @info "Have ifelse statements which does not contain time. Hence we do not rewrite as event, but rather keep it as an ifelse." maxlog=1
             rewrite_ifelse = false
             continue
         elseif time_left == true
diff --git a/src/SBML/SBML_to_ModellingToolkit.jl b/src/SBML/SBML_to_ModellingToolkit.jl
index cde1bb7c..7aae94ed 100644
--- a/src/SBML/SBML_to_ModellingToolkit.jl
+++ b/src/SBML/SBML_to_ModellingToolkit.jl
@@ -48,6 +48,10 @@ function as_trigger(trigger_formula, model_dict, model_SBML)
         parts[1] = replace_whole_word(parts[1], "time", "t")
     end
 
+    # Account for potential reaction-math kinetics making out the trigger 
+    parts[1] = replace_reactionid_with_math(parts[1], model_SBML)
+    parts[2] = replace_reactionid_with_math(parts[2], model_SBML)
+
     # States in ODE-system are typically in substance units, but formulas in 
     # concentratio. Thus, each state is divided with its corresponding 
     # compartment 
@@ -104,7 +108,10 @@ function process_initial_assignment(model_SBML, model_dict::Dict, base_functions
     initally_assigned_parameter = Dict{String, String}()
     for (assignId, initialAssignment) in model_SBML.initial_assignments
         
+        # The formula might refer to reaciton, in this case insert reaction dynamics 
         _formula = SBML_math_to_str(initialAssignment)
+        _formula = replace_reactionid_with_math(_formula, model_SBML)
+        
         formula = rewrite_derivatives(_formula, model_dict, base_functions, model_SBML)
         # Initial time i zero 
         formula = replace_whole_word(formula, "t", "0.0")
@@ -250,6 +257,13 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
         if model_dict["isBoundaryCondition"][state_id] == true
            model_dict["derivatives"][state_id] *= "0.0"
         end
+
+        # In case the conc. is given in initial conc, but the state should be in amounts this 
+        # must be acounted for with initial values
+        if model_dict["stateGivenInAmounts"][state_id][1] == false && model_dict["hasOnlySubstanceUnits"][state_id] == true
+            model_dict["stateGivenInAmounts"][state_id] = (true, state.compartment)
+            model_dict["states"][state_id] = string(state.initial_concentration) * " * " * state.compartment
+        end
     end
 
     # Extract model parameters and their default values. In case a parameter is non-constant 
@@ -306,6 +320,8 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
             event_assign_to[i] = event_assignment.variable
             event_formulas[i] = replace_function_with_formula(SBML_math_to_str(event_assignment.math), model_dict["modelFunctions"])
             event_formulas[i] = replace_whole_word(event_formulas[i], "t", "integrator.t")
+            event_formulas[i] = replace_reactionid_with_math(event_formulas[i], model_SBML)
+
             # Species typically given in substance units, but formulas in conc. Thus we must account for assignment 
             # formula being in conc., but we are changing something by amount 
             if event_assign_to[i] ∈ keys(model_SBML.species)
@@ -405,11 +421,19 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
             model_dict["derivatives"][product.species] *= "+" * stoichiometry * compartment_scaling * "(" * formula * ")"
         end
     end
-    # For states given in amount but model equations are in conc., multiply with compartment 
+    # For states given in amount but model equations are in conc., multiply with compartment, also handle potential 
+    # reaction identifayers in the derivative
     for (state_id, derivative) in model_dict["derivatives"]
+
+        model_dict["derivatives"][state_id] = replace_reactionid_with_math(model_dict["derivatives"][state_id], model_SBML)
+
         if model_dict["stateGivenInAmounts"][state_id][1] == false
             continue
         end
+        # Here equations should be given in amounts 
+        if model_dict["hasOnlySubstanceUnits"][state_id] == true
+            continue
+        end
         # Algebraic rule (see below)
         if replace(derivative, " " => "")[end] == '~' || replace(derivative, " " => "")[end] == '0'
             continue
@@ -420,7 +444,13 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
 
     # For states given by assignment rules 
     for (state, formula) in model_dict["assignmentRulesStates"]
-        model_dict["derivatives"][state] = state * " ~ " * formula
+        _formula = rewrite_derivatives(formula, model_dict, base_functions, model_SBML; check_scaling=true)
+        # Must track if species is given in amounts or conc.
+        if state ∈ keys(model_SBML.species) && model_SBML.species[state].only_substance_units == false
+            cmult = model_dict["stateGivenInAmounts"][state][1] == true ? " * " * model_SBML.species[state].compartment : ""
+            _formula = "(" * _formula * ")" * cmult
+        end
+        model_dict["derivatives"][state] = state * " ~ " * _formula
         if state ∈ non_constant_parameter_names
             delete!(model_dict["states"], state)
             delete!(model_dict["parameters"], state)
@@ -509,7 +539,17 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
         if divide_with_compartment == false
             continue
         end
-        model_dict["derivatives"][specie] = specie * " ~ (" * model_dict["states"][specie] * ") / " * c
+        if model_dict["stateGivenInAmounts"][specie][1] == true
+            model_dict["derivatives"][specie] = specie * " ~ (" * model_dict["states"][specie] * ") / " * c
+        else
+            # Must account for the fact that the compartment can change in size, and thus need to 
+            # account for its initial value 
+            if c ∈ keys(model_dict["parameters"])
+                model_dict["derivatives"][specie] = specie * " ~ (" * model_dict["states"][specie] * ")"
+            else
+                model_dict["derivatives"][specie] = specie * " ~ (" * model_dict["states"][specie] * ")" * " * " * string(model_dict["states"][c]) * " / " * c 
+            end
+        end
     end
 
     # In case the model has a conversion factor 
@@ -541,6 +581,35 @@ function build_model_dict(model_SBML, ifelse_to_event::Bool)
         non_constant_parameter_names = filter(x -> x != id, non_constant_parameter_names)
     end
 
+    #=
+        Sometimes the volume might change over time but the amount should stay constant, as we have a boundary condition
+        In this case it follows that amount n (amount), V (compartment) and conc. are related 
+        via the chain rule by - I need to change my ODE:s and add state
+        dn/dt = d(n/V)/dt*V + n*dV/dt 
+    =#
+    for specie in keys(model_SBML.species)
+
+        if !(model_dict["stateGivenInAmounts"][specie][1] == true && 
+             model_SBML.species[specie].compartment ∈  rate_rules_names && 
+             model_dict["isBoundaryCondition"][specie] == true &&
+             specie ∈ rate_rules_names && 
+             model_dict["hasOnlySubstanceUnits"][specie] == false)
+            continue
+        end
+
+        # Derivative and inital values for concentratin species
+        compartment = model_SBML.species[specie].compartment
+        specie_conc = "__" * specie * "__conc__"
+        model_dict["states"][specie_conc] = model_dict["states"][specie] * " / " * compartment
+        i_start = findfirst(x -> x == '~', model_dict["derivatives"][specie]) + 1
+        i_end = findlast(x -> x == ')', model_dict["derivatives"][specie])
+        model_dict["derivatives"][specie_conc] = "D(" * specie_conc  * ") ~ " * model_dict["derivatives"][specie][i_start:i_end]
+
+        # Rebuild derivative for amount specie
+        itmp1, itmp2 = findfirst(x -> x == '~', model_dict["derivatives"][specie_conc]) + 1, findfirst(x -> x == '~', model_dict["derivatives"][compartment]) + 1
+        model_dict["derivatives"][specie] = "D(" * specie * ") ~ " * model_dict["derivatives"][specie_conc][itmp1:end] * "*" *  compartment * " + " * specie * "*" * model_dict["derivatives"][compartment][itmp2:end] * " / " * compartment
+    end
+
     model_dict["numOfParameters"] = string(length(keys(model_dict["parameters"])))
     model_dict["numOfSpecies"] = string(length(keys(model_dict["states"])))
     model_dict["non_constant_parameter_names"] = non_constant_parameter_names
@@ -676,11 +745,12 @@ function create_ode_model(model_dict, path_jl_file, model_name, write_to_file::B
             s_index += 1
         end
     end
-    for key in keys(model_dict["assignmentRulesStates"])
+    for (key, formula) in model_dict["assignmentRulesStates"]
+        what_write = key ∈ keys(model_dict["derivatives"]) ? model_dict["derivatives"][key] : key * " ~ " * model_dict["assignmentRulesStates"][key]
         if s_index != 1
-            dict_model_str["derivatives"] *= ",\n    " * key * " ~ " * model_dict["assignmentRulesStates"][key]
+            dict_model_str["derivatives"] *= ",\n    " * what_write
         else
-            dict_model_str["derivatives"] *= "    " * key * " ~ " * model_dict["assignmentRulesStates"][key]
+            dict_model_str["derivatives"] *= "    " * what_write
             s_index += 1
         end
     end
@@ -923,3 +993,12 @@ function _SBML_math_to_str(math::SBML.MathConst)
         return math.id, false
     end
 end
+
+
+function replace_reactionid_with_math(formula::T, model_SBML)::T where T<:AbstractString
+    for (reaction_id, reaction) in model_SBML.reactions
+        reaction_math = SBML_math_to_str(reaction.kinetic_math)
+        formula = replace_whole_word(formula, reaction_id, reaction_math)
+    end
+    return formula
+end
diff --git a/test/Test_SBML.jl b/test/Test_SBML.jl
index 241d496b..dfe582c7 100644
--- a/test/Test_SBML.jl
+++ b/test/Test_SBML.jl
@@ -17,10 +17,11 @@ using DataFrames
     goes from false to true 
 =#
 
-
-testCase = "00001"
-testCase = "00607"
-testSBMLTestSuite(testCase, Rodas4P())
+# 01064 has stochiometry math 
+testCase = "00295"
+testCase = "00051"
+testCase = "01233"
+#testSBMLTestSuite(testCase, Rodas4P())
 # Next we must allow species to first be defined via an InitialAssignment, pretty stupied to me, but aja...
 function testSBMLTestSuite(testCase, solver)
     @info "Test case $testCase"
@@ -34,10 +35,12 @@ function testSBMLTestSuite(testCase, solver)
     if testCase == "00995" || testCase == "00996" || testCase == "00997"
         expected = expected[2:end, :]
     end
+    col_names =  Symbol.(replace.(string.(names(expected)), " " => ""))
+    rename!(expected, col_names)
 
-    t_save = Float64.(expected[!, :time])
+    t_save = "Time" in names(expected) ? Float64.(expected[!, :Time]) : Float64.(expected[!, :time])
     tmax = maximum(t_save)
-    whatCheck = Symbol.(filter(x -> x != "time", names(expected)))
+    whatCheck = filter(x -> x ∉ [:time, :Time], col_names)
     path_SBML = path_SBMLFiles[end]
 
     # Read settings file 
@@ -74,31 +77,41 @@ function testSBMLTestSuite(testCase, solver)
                 continue
             end
 
-            if toCheck ∈ speciesTest && toCheck ∈  speciesTestConc
-
+            if toCheck ∈ speciesTest && toCheck ∈ speciesTestConc && string(toCheck) ∈ keys(model_SBML.species)
                 compartmentName = model_SBML.species[string(toCheck)].compartment
-                if model_dict["stateGivenInAmounts"][string(toCheck)][1] == false
+                if model_dict["stateGivenInAmounts"][string(toCheck)][1] == false 
                     c = 1.0
                 elseif compartmentName in string.(model_parameters)
                     c = sol.prob.p[findfirst(x -> x == compartmentName, string.(model_parameters))]
                 else
                     c = sol[Symbol(compartmentName)]
                 end
+            elseif toCheck ∈ speciesTest && toCheck ∈ speciesTestAmount && string(toCheck) ∈ keys(model_SBML.species)
+                compartmentName = model_SBML.species[string(toCheck)].compartment
+                if model_dict["stateGivenInAmounts"][string(toCheck)][1] == false 
+                    if compartmentName in string.(model_parameters)
+                        c = 1 / (sol.prob.p[findfirst(x -> x == compartmentName, string.(model_parameters))])
+                    else
+                        c = 1 ./ sol[Symbol(compartmentName)]
+                    end
+                else
+                    c = 1.0
+                end
             else
                 c = 1.0
             end
-
+            
             @test all(abs.(sol[toCheck] ./ c .- expected[!, toCheck]) .< absTolTest .+ relTolTest .* abs.(expected[!, toCheck]))
         end
     end
 end
 
 
-# 00945 crash!!
+# 01014 current max 
 # 00369
 solver = Rodas4P()
 @testset "SBML test suite" begin
-    for i in 1:997
+    for i in 1:1236
         testCase = repeat("0", 5 - length(string(i))) *  string(i)
 
         if testCase == "00028"
@@ -106,11 +119,15 @@ solver = Rodas4P()
             continue
         end
 
-        # StoichiometryMath we do not aim to support 
+        # StoichiometryMath we do not yet support
         if testCase ∈ ["00068", "00069", "00070", "00129", "00130", "00131", "00388", "00391", "00394", "00516", 
                        "00517", "00518", "00519", "00520", "00521", "00522", "00561", "00562", "00563", 
                        "00564", "00731", "00827", "00828", "00829", "00898", "00899", "00900", "00609", 
-                       "00610", "00968", "00973", "00989", "00990", "00991", "00992", "00993", "00994"]
+                       "00610", "00968", "00973", "00989", "00990", "00991", "00992", "00993", "00994", 
+                       "01027", "01028", "01029", "01064", "01066", "01064", "01069", "01071", "01073", 
+                       "01084", "01085", "01086", "01088", "01095", "01096", "01097", "01100", "01101", 
+                       "01103", "01104", "01105", "01106", "01107", "01108", "01109", "01110", "01111", 
+                       "01121"]
             continue
         end
 
@@ -119,6 +136,18 @@ solver = Rodas4P()
             continue
         end
 
+        # We and SBML.jl do not currently support hierarchical models
+        not_test = ["011" * string(i) for i in 26:83]
+        if testCase ∈ not_test
+            continue
+        end
+
+        # We do not aim to support Flux-Balance-Analysis (FBA) models
+        not_test = ["01" * string(i) for i in 186:197]
+        if testCase ∈ not_test
+            continue
+        end
+
         # If user wants to add a random species, it must either be as a species, initialAssignment, assignmentRule
         # or by event, not by just random adding it to equations.
         if testCase ∈ ["00974"]
@@ -129,7 +158,7 @@ solver = Rodas4P()
         # As of yet we do not support events with priority, but could if there are interest. However should
         # be put up as an issue on GitHub 
         if testCase ∈ ["00931", "00934", "00935", "00962", "00963", "00964", "00965", "00966", "00967", 
-                       "00978", "00978"]
+                       "00978", "00978", "01229"]
             continue
         end
 
@@ -151,7 +180,22 @@ solver = Rodas4P()
 
         # Fast reactions can technically be handled via algebraic rules, will add support if wanted 
         if testCase ∈ ["00870", "00871", "00872", "00873", "00874", "00875", "00986", "00987", 
-                       "00988"]
+                       "00988", "01051", "01052", "01053"]
+            continue
+        end
+
+        # We do not support an event with multiple triggers 
+        if testCase ∈ ["01211"]
+            continue
+        end
+
+        # We do not support an event with piecewise in the activation 
+        if testCase ∈ ["01212", "01213", "01214", "01215"]
+            continue
+        end
+
+        # We do not lt etc... with multiple pair (more than two) parameters 
+        if testCase ∈ ["01216"]
             continue
         end
 
@@ -168,6 +212,12 @@ solver = Rodas4P()
             continue
         end
 
+        # Piecewise in initialAssignments we do not aim to support (can easily be 
+        # side-stepeed with assignmentrules)
+        if testCase ∈ ["01112", "01113", "01114", "01115", "01116", "01208", "01209", "01210"]
+            continue
+        end
+
         # We do not allow stochastic simulations 
         if testCase ∈ ["00952", "00953"]
             continue
@@ -186,53 +236,11 @@ solver = Rodas4P()
                        "00770", "00771", "00772", "00773", "00774", "00775", "00776", 
                        "00777", "00778", "00779", "00780", "00848", "00849", "00850", 
                        "00886", "00887", "00932", "00933", "00936", "00408", "00461", 
-                       "00655", "00656", "00657", "00980", ]) || testCase ∈ notTest
+                       "00655", "00656", "00657", "00980", "01000", "01048", "01049", 
+                       "01050", "01074", "01075", "01076", "01119", "01120", "01230"]) || testCase ∈ notTest
             continue
         end
 
         testSBMLTestSuite(testCase, solver)
     end
 end
-
-
-
-using OrdinaryDiffEq
-function f(du, u, p, t)
-    du[1] = -u[1]*p[1]
-end
-u0 = [10.0]
-const V = 1
-prob = ODEProblem(f, u0, (0.0, 10.0), [1.0])
-
-function condition_test(u, t, integrator)
-    
-    if integrator.p[1] > 0 && integrator.p[1] < 10
-        println("t = ", t)
-        return true
-    end
-    return false
-end
-affect!(integrator) = integrator.p[1] = 10
-function testinit(c,u,t,integrator)
-    cond = condition_test(u, t, integrator)
-    if cond == true
-        println("t in init = $t")
-        affect!(integrator)
-    end
-end
-
-cb = DiscreteCallback(condition_test, affect!, initialize=testinit)
-
-sol = solve(prob, Tsit5(), callback = cb)
-
-
-function hej(x, y)
-    println("y[1] = ", y[1])
-    y[1] += 1
-    return x * y[1]
-end
-
-
-_hej = let y=[1.0]
-    (x) -> hej(x, y)
-end
\ No newline at end of file