diff --git a/Project.toml b/Project.toml
index d298d37..6ba044c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SBMLToolkit"
 uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 authors = ["paulflang", "anandijain"]
-version = "0.1.25"
+version = "0.1.26"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
diff --git a/README.md b/README.md
index 4d1f83a..355b430 100644
--- a/README.md
+++ b/README.md
@@ -27,28 +27,43 @@ Pkg.add("SBMLToolkit")
 SBML models can be simulated with the following steps (note that `sol` is in absolute quantities rather than concentrations):
 
 ```julia
-using SBMLToolkit, ModelingToolkit, OrdinaryDiffEq
+using SBMLToolkit, OrdinaryDiffEq
 
-SBMLToolkit.checksupport_file("my_model.xml")
-mdl = readSBML("my_model.xml", doc -> begin
-    set_level_and_version(3, 2)(doc)
-    convert_promotelocals_expandfuns(doc)
-end)
-
-rs = ReactionSystem(mdl)  # If you want to create a reaction system
-odesys = convert(ODESystem, rs)  # Alternatively: ODESystem(mdl)
+odesys = readSBML("my_model.xml", ODESystemImporter())
 
 tspan = (0.0, 1.0)
 prob = ODEProblem(odesys, [], tspan, [])
 sol = solve(prob, Tsit5())
 ```
 
-Alternatively, SBMLToolkit also provides more concise methods to import `SBML.Models`, `Catalyst.ReactionSystems`, and `ModelingToolkit.ODESystems`.
+While this imports an `ODESystem` directly, you can also import a Catalyst.jl `ReactionSystem`:
 
 ```julia
-mdl = readSBML("my_model.xml", DefaultImporter())
+using SBMLToolkit
+
 rs = readSBML("my_model.xml", ReactionSystemImporter())
-odesys = readSBML("my_model.xml", ODESystemImporter())
+```
+
+One common case where this is useful is if you want to run stochastic instead of ODE simulations.
+
+In the very unlikely case that you need fine-grained control over the SBML import, you can create an SBML.jl `Model` (we strongly recommend manually running `SBMLToolkit.checksupport_file("my_model.xml")` before)
+
+```julia
+using SBML
+
+mdl = readSBML("my_model.xml", doc -> begin
+    set_level_and_version(3, 2)(doc)
+    convert_promotelocals_expandfuns(doc)
+end)
+```
+
+The conversion to SBML level 3 version 2 is necessary, because older versions are not well supported in SBMLToolkit. `convert_promotelocals_expandfuns` basically flattens the SBML before the import. Once you have obtained the `Model`, you can convert it to a `ReactionSystem` and `ODESystem`.
+
+```julia
+using SBMLToolkit
+
+rs = ReactionSystem(mdl)
+odesys = convert(ODESystem, rs)
 ```
 
 ## License
diff --git a/src/reactions.jl b/src/reactions.jl
index 86e3409..86ff18c 100644
--- a/src/reactions.jl
+++ b/src/reactions.jl
@@ -89,7 +89,7 @@ function get_reagents(reactant_references::Vector{SBML.SpeciesReference},
         sn = rr.species
         stoich = rr.stoichiometry
         if isnothing(stoich)
-            @warn "Stoichiometries of SpeciesReferences are not defined. Setting to 1." maxlog=1
+            @warn "SBML SpeciesReferences does not contain stoichiometries. Assuming stoichiometry of 1." maxlog=1
             stoich = 1.0
         end
         iszero(stoich) && @error("Stoichiometry of $sn must be non-zero")
@@ -183,6 +183,8 @@ function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing},
         rate_const = kl
     elseif isnothing(reactants) | isnothing(stoich)
         throw(ErrorException("`reactants` and `stoich` are inconsistent: `reactants` are $(reactants) and `stoich` is $(stoich)."))
+    elseif max(stoich...) > 100  # simplify_fractions might StackOverflow
+        rate_const = kl
     else
         rate_const = SymbolicUtils.simplify_fractions(kl / *((.^(reactants, stoich))...))
     end
diff --git a/test/reactions.jl b/test/reactions.jl
index 0dafb76..a39147e 100644
--- a/test/reactions.jl
+++ b/test/reactions.jl
@@ -135,6 +135,7 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r))
 @test isequal(SBMLToolkit.get_massaction(k1 + c1, nothing, nothing), k1 + c1)  # Case zero order kineticLaw
 @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2 / (c1 + s2), [s1], [1]))  # Case Michaelis-Menten kinetics
 @test isnan(SBMLToolkit.get_massaction(k1 * s1 * IV, [s1], [1]))  # Case kineticLaw with time
+@test isnan(SBMLToolkit.get_massaction(k1 * s1, [s1], [101]))  # Case no simplification performed due to large rstoich
 
 @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2, [s1], [1]))
 @test isnan(SBMLToolkit.get_massaction(k1 + c1, [s1], [1]))