diff --git a/src/DiffEqBiological.jl b/src/DiffEqBiological.jl
index 513027c7e7..10d3f64ae3 100644
--- a/src/DiffEqBiological.jl
+++ b/src/DiffEqBiological.jl
@@ -25,8 +25,9 @@ export @reaction_network, @reaction_func, @min_reaction_network
 
 # functions to query network properties
 export speciesmap, paramsmap, numspecies, numreactions, numparams
-export odeexprs, jacobianexprs, noiseexprs, jumpexprs
-export get_substrate_stoich, get_net_stoich
+export oderhsfun, jacfun, paramjacfun, odefun, noisefun, sdefun, jumps, regularjumps
+export odeexprs, jacobianexprs, noiseexprs, jumpexprs, rateexpr, oderatelawexpr, ssaratelawexpr
+export substratestoich, productstoich, netstoich, ismassaction, dependants, dependents, substrates, products
 export rxtospecies_depgraph, speciestorx_depgraph, rxtorx_depgraph
 
 # functions to add mathematical equations to the network
diff --git a/src/massaction_jump_utils.jl b/src/massaction_jump_utils.jl
index 9bb1605a6d..393a442a46 100644
--- a/src/massaction_jump_utils.jl
+++ b/src/massaction_jump_utils.jl
@@ -3,8 +3,8 @@
 
 # given a ReactionStruct and a species map construct a MassActionJump
 function make_majump(rs, specmap, ratemap, params, param_context)
-    reactant_stoich = get_substrate_stoich(rs, specmap)
-    net_stoich      = get_net_stoich(rs, specmap)
+    reactant_stoich = substratestoich(rs, specmap)
+    net_stoich      = netstoich(rs, specmap)
     if isempty(net_stoich)
         error("Empty net stoichiometry vectors for mass action reactions are not allowed.")
     end
@@ -92,7 +92,7 @@ function jump_to_dep_specs_map(rn, rxidxs_jidxs)
     jtos_vec = Vector{Vector{valtype(specmap)}}(undef, numrxs)
     for rx in 1:numrxs
         jidx           = rxidxs_jidxs[rx]
-        jtos_vec[jidx] = sort!( [ns.first for ns in get_net_stoich(rn.reactions[rx], specmap)] )
+        jtos_vec[jidx] = sort!( [ns.first for ns in netstoich(rn.reactions[rx], specmap)] )
     end
 
     jtos_vec
@@ -125,7 +125,7 @@ function depgraph_from_network(rn, jset, rxidxs_to_jidxs, spec_to_dep_jumps)
         jidx = rxidxs_to_jidxs[rx]
 
         # get the net reaction stoichiometry
-        net_stoich = get_net_stoich(rn.reactions[rx], speciesmap(rn))
+        net_stoich = netstoich(rn.reactions[rx], speciesmap(rn))
 
         # rx changes spec, hence rxs depending on spec depend on rx
         for (spec,stoch) in net_stoich
diff --git a/src/network_properties.jl b/src/network_properties.jl
index d2eaf64f02..64176594be 100644
--- a/src/network_properties.jl
+++ b/src/network_properties.jl
@@ -1,39 +1,187 @@
-""" 
-Functions for querying network properties.
-"""
+# Functions for querying network properties.
 
 ######### Accessors: #########
 
 """
-Return a Dictionary mapping from species symbol to species index.
+    speciesmap(network)
+
+Given an `AbstractReactionNetwork`, return a Dictionary mapping from species
+symbol to species index. 
 """
 function speciesmap(network)
     network.syms_to_ints
 end
 
 """
-Return a Dictionary mapping from parameter symbol to parameter index.
+    paramsmap(network)
+
+Given an `AbstractReactionNetwork`, return a Dictionary mapping from parameter
+symbol to parameter index.
 """
 function paramsmap(network)
     network.params_to_ints
 end
 
+"""
+    numspecies(network)
+
+Return the number of species within the given `AbstractReactionNetwork`.
+"""
 function numspecies(network)
     length(speciesmap(network))
 end
 
+"""
+    numreactions(network)
+
+Return the number of reactions within the given `AbstractReactionNetwork`.
+"""
 function numreactions(network)
     length(network.reactions)
 end
 
+"""
+    numparams(network)
+
+Return the number of parameters within the given `AbstractReactionNetwork`.
+"""
 function numparams(network)
     length(paramsmap(network))
 end
 
+######### Generated Functions: #########
+
+"""
+    oderhsfun(network)
+
+Given an `AbstractReactionNetwork`, return a function, `f!(du,u,p,t)`, that
+evaluates the current value of the ODE model derivative functions, ``du/dt = f(u,t)``, 
+within `du`.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addodes!`
+must be called first.
+"""
+function oderhsfun(network)
+    isnothing(network.f) && error("Error, call addodes! first.")
+    network.f
+end
+
+"""
+    jacfun(network)
+
+Given an `AbstractReactionNetwork`, return a function, `jac!(J,u,p,t)`, that
+evaluates the current Jacobian matrix, `J`, of the ODE model, ``du/dt = f(u,t)``. 
+The Jacobian matrix has entries 
+
+``J_{i,j} = \\partial f_i(u,t) / \\partial u_j``.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addodes!` 
+must be called first.
+"""
+function jacfun(network)
+    isnothing(network.jac) && error("Error, call addodes! first.")
+    network.jac
+end
+
+"""
+    paramjacfun(network)
+
+Given an `AbstractReactionNetwork`, return a function, `pjac(pJ,u,p,t)`, that
+evaluates the current parameter Jacobian matrix, `pJ`, of the ODE model, ``du/dt = f(u,t)``. 
+The parameter Jacobian matrix has entries 
+
+``pJ_{i,j} = \\partial f_i(u,t) / \\partial p_j``.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addodes!` 
+must be called first.
+"""
+function paramjacfun(network)
+    isnothing(network.paramjac) && error("Error, call addodes! first.")
+    network.paramjac
+end
+
+"""
+    odefun(network)
+
+Given an `AbstractReactionNetwork`, return a `DiffEqBase.ODEFunction` encoding
+an ODE model for the reaction network. 
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addodes!` 
+must be called first.
+"""
+function odefun(network)
+    isnothing(network.odefun) && error("Error, call addodes! first.")
+    network.odefun
+end
+
+"""
+    noisefun(network)
+
+Given an `AbstractReactionNetwork`, return a function, `g(η,u,p,t)`, that
+evaluates the current noise coefficients for each reaction in the Chemical
+Langevin Equation representation within `η`.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addsdes!` 
+must be called first.
+"""
+function noisefun(network)
+    isnothing(network.g) && error("Error, call addsdes! first.")
+    network.g 
+end
+
+"""
+    sdefun(network)
+
+Given an `AbstractReactionNetwork`, return a `DiffEqBase.SDEFunction` encoding
+a Chemical Langevin Equation SDE model for the reaction network. 
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addsdes!` 
+must be called first.
+"""
+function sdefun(network)
+    isnothing(network.sdefun) && error("Error, call addsdes! first.")
+    network.sdefun
+end
+
+"""
+    jumps(network)
+
+Given an `AbstractReactionNetwork`, return a tuple of `AbstractJumps` encoding
+a stochastical chemical kinetics representation for the reaction network.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addjumps!` 
+must be called first.
+"""
+function jumps(network)
+    isnothing(network.jumps) && error("Error, call addjumps! first.")
+    network.jumps
+end
+
+"""
+    regularjumps(network)
+
+Given an `AbstractReactionNetwork`, return a `RegularJump` encoding a
+stochastical chemical kinetics representation of the reaction network for use in
+``\\tau``-leaping approximations.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addjumps!` 
+must be called first.
+"""
+function regularjumps(network)
+    isnothing(network.regular_jumps) && error("Error, call addjumps! first.")
+    network.regular_jumps
+end
+
+
 ######### Generated Expressions: #########
 
 """
-Return a vector of the ODE derivative expressions.
+    odeexprs(network)
+
+Given an `AbstractReactionNetwork`, return a vector of the ODE expressions.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addodes!` 
+must be called first.
 """
 function odeexprs(network)
     isnothing(network.f_func) && error("Error, call addodes! first.")
@@ -41,7 +189,13 @@ function odeexprs(network)
 end
 
 """
-Return a matrix with the ODE Jacobian expressions.
+    jacobianexprs(network)
+
+Given an `AbstractReactionNetwork`, return a matrix with the ODE Jacobian
+expressions.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addodes!` 
+must be called first.
 """
 function jacobianexprs(network)
     isnothing(network.symjac) && error("Error, call addodes! first.")
@@ -49,7 +203,13 @@ function jacobianexprs(network)
 end
 
 """
-Return a vector of the SDE noise expressions for each reaction.
+    noiseexprs(network)
+
+Given an `AbstractReactionNetwork`, return a vector of the SDE noise expressions
+for each reaction.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addsdes!` 
+must be called first.
 """
 function noiseexprs(network)
     isnothing(network.g_func) && error("Error, call addsdes! first.")
@@ -57,32 +217,108 @@ function noiseexprs(network)
 end
 
 """
-Return a tuple of the jump rates and affects expressions
+    jumpexprs(network)
+
+Given an `AbstractReactionNetwork`, return a tuple of the jump rates and affects
+expressions.
+
+*Note,* for a network generated with the `@min_reaction_network` macro `addjumps!` 
+must be called first.
 """
 function jumpexprs(network)
     (isnothing(network.jump_rate_expr) || isnothing(network.jump_affect_expr)) && error("Error, call addjumps! first.")
     network.jump_rate_expr,network.jump_affect_expr
 end
 
+"""
+    rateexpr(network, rxidx)
+
+Given an `AbstractReactionNetwork`, return the reaction rate expression for the
+reaction with index `rxidx`. Note, for a reaction defined by
 
-######### Reaction Properties: #########
+`k*X*Y, X+Z --> 2X + Y`
 
+the expression that is returned will be `:(k*X*Y)`, while the *rate law* used in
+ODEs, SDEs and jumps would be `k*X^2*Y*Z`.
 """
-Return a Vector of pairs, mapping ids of species that serve as substrates 
-in the reaction to the corresponding stoichiometric coefficient as a substrate.
+function rateexpr(network, rxidx)
+    network.reactions[rxidx].rate_org
+end
+
+"""
+    oderatelawexpr(network, rxidx)
+
+Given an `AbstractReactionNetwork`, return the reaction rate law expression used
+in generated ODEs for the reaction with index `rxidx`. Note, for a reaction
+defined by
+
+`k*X*Y, X+Z --> 2X + Y`
+
+the expression that is returned will be `:(k*X^2*Y*Z)`. For a reaction of the
+form 
+
+`k, 2X+3Y --> Z`
+
+the expression that is returned will be `:(k * (X^2/2) * (Y^3/6))`.
+""" 
+function oderatelawexpr(network, rxidx)
+    network.reactions[rxidx].rate_DE
+end
+
+"""
+    ssaratelawexpr(network, rxidx)
+
+Given an `AbstractReactionNetwork`, return the reaction rate law expression used
+in generated stochastic chemical kinetic model SSAs for the reaction with index
+`rxidx`. Note, for a reaction defined by
+
+`k*X*Y, X+Z --> 2X + Y`
+
+the expression that is returned will be `:(k*X^2*Y*Z)`. For a reaction of the
+form 
+
+`k, 2X+3Y --> Z`
+
+the expression that is returned will be `:(k * binomial(X,2) * binomial(Y,3))`.
+""" 
+function ssaratelawexpr(network, rxidx)
+    network.reactions[rxidx].rate_SSA
+end
+
+######### Reaction Properties: #########
+
+function substratestoich(rs::ReactionStruct, specmap)
+    sort!( [specmap[s.reactant] => s.stoichiometry for s in rs.substrates] )
+end
+
 """
-function get_substrate_stoich(rs::ReactionStruct, specmap)
-    reactant_stoich = [specmap[s.reactant] => s.stoichiometry for s in rs.substrates]
-    sort!(reactant_stoich)
-    reactant_stoich
+    substratestoich(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+vector of pairs, mapping ids of species that serve as substrates in the reaction
+to the corresponding stoichiometric coefficient as a substrate. 
+""" 
+function substratestoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    substratestoich(rn.reactions[rxidx], speciesmap(rn))
+end
+
+
+function productstoich(rs::ReactionStruct, specmap)
+    sort( [specmap[p.reactant] => p.stoichiometry for p in rs.products] )    
 end
 
 """
-Return a Vector of pairs, mapping ids of species that participate in the 
-reaction to the net stoichiometric coefficient of the species (i.e. net change
-in the species due to the reaction).
+    productstoich(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+vector of pairs, mapping ids of species that are products in the reaction to the
+corresponding stoichiometric coefficient as a product.
 """
-function get_net_stoich(rs::ReactionStruct, specmap)
+function productstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    productstoich(rn.reactions[rxidx], speciesmap(rn))
+end
+
+function netstoich(rs::ReactionStruct, specmap)
     nsdict = Dict{Int,Int}(specmap[s.reactant] => -s.stoichiometry for s in rs.substrates)
 
     for prod in rs.products
@@ -98,22 +334,111 @@ function get_net_stoich(rs::ReactionStruct, specmap)
     net_stoich
 end
 
+"""
+    netstoich(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+vector of pairs, mapping ids of species that participate in the reaction to the
+net stoichiometric coefficient of the species (i.e. net change in the species
+due to the reaction).
+"""
+function netstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    netstoich(rn.reactions[rxidx], speciesmap(rn))
+end
+
+"""
+    ismassaction(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+boolean indicating whether the given reaction is of mass action form. For
+example, the reaction
+
+`2*k, 2X + 3Y --> 5Z + W`
+
+would return true, while reactions with state-dependent rates like
+
+`k*X, X + Y --> Z`
+
+would return false.
+"""
+function ismassaction(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    rn.reactions[rxidx].is_pure_mass_action
+end
+
+
+"""
+    dependents(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+vector of symbols of species the *reaction rate law* depends on. i.e. for
+
+`k*W, 2X + 3Y --> 5Z + W`
+
+the returned vector would be `[:W,:X,:Y]`.
+"""
+function dependents(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    rn.reactions[rxidx].dependants
+end
+
+"""
+    dependants(network, rxidx)
+
+See documentation for [`dependents(network, rxidx)`](@ref).
+"""
+function dependants(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    dependents(rn, rxidx)
+end
+
+"""
+    substrates(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+vector of symbols of species that correspond to substrates in the reaction. 
+i.e. for
+
+`k*W, X + 3Y --> X + W`
+
+the returned vector would be `[:X,:Y]`.
+"""
+function substrates(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    rn.reactions[rxidx].substrates
+end
+
+"""
+    products(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+vector of symbols of species that correspond to products in the reaction. 
+i.e. for
+
+`k*W, X + 3Y --> X + W`
+
+the returned vector would be `[:X,:W]`.
+"""
+function products(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    rn.reactions[rxidx].products
+end
+
 ######### Network Properties: #########
 
 """
-Returns a Vector{Vector{Int}} mapping a reaction index 
-to the indices of species that depend on it.
+    rxtospecies_depgraph(network)
+
+Given an `AbstractReactionNetwork`, returns a Vector{Vector{Int}} mapping a
+reaction index to the indices of species that depend on it.
 """
 function rxtospecies_depgraph(network)
     specmap = speciesmap(network)
 
     # map from a reaction to vector of species that depend on it
-    [sort!( [ns.first for ns in get_net_stoich(rx, specmap)] ) for rx in network.reactions]
+    [sort!( [ns.first for ns in netstoich(rx, specmap)] ) for rx in network.reactions]
 end
 
 """
-Returns a Vector{Vector{Int}} mapping a species index 
-to the indices of reactions that depend on it.
+    speciestorx_depgraph(network)
+
+Given an `AbstractReactionNetwork`, returns a Vector{Vector{Int}} mapping a
+species index to the indices of reactions that depend on it.
 """
 function speciestorx_depgraph(network)
     numrxs  = numreactions(network)    
@@ -122,7 +447,7 @@ function speciestorx_depgraph(network)
     # map from a species to reactions that depend on it
     spectorxs = [Vector{Int}() for i=1:numspecies(network)]
     for rx in 1:numrxs
-        for specsym in network.reactions[rx].dependants
+        for specsym in dependents(network, rx) 
             push!(spectorxs[specmap[specsym]], rx)
         end
     end
@@ -132,17 +457,18 @@ function speciestorx_depgraph(network)
 end
 
 """
-Returns a Vector{Vector{Int}} mapping a reaction index
-to the indices of reactions that depend on it.
+    rxtorx_depgraph(network)
+
+Given an `AbstractReactionNetwork`, returns a Vector{Vector{Int}} mapping a
+reaction index to the indices of reactions that depend on it.
 """
-function rxtorx_depgraph(network, sptorxsmap=nothing)
-    sptorxs = isnothing(sptorxsmap) ? speciestorx_depgraph(network) : sptorxsmap
+function rxtorx_depgraph(network, sptorxs=speciestorx_depgraph(network))
     numrxs = numreactions(network)
     dep_sets = [SortedSet{Int}() for n = 1:numrxs]
 
     # dg as vector of sets
     for rx in 1:numrxs
-        net_stoich = get_net_stoich(network.reactions[rx], speciesmap(network))
+        net_stoich = netstoich(network, rx)
 
         for (spec,stoich) in net_stoich
             for dependent_rx in sptorxs[spec]