diff --git a/HISTORY.md b/HISTORY.md
index d786a47bdd..02ac9a266e 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -24,6 +24,13 @@
   ```
   X's value will be `1.0` in the first vertex, but `0.0` in the remaining one (the system have 25 vertexes in total). SInce th parameters `p` and `d` are part of the non-spatial reaction network, their values are tied to vertexes. However, if the `D` parameter (which governs diffusion between vertexes) is given several values, these will instead correspond to the specific edges (and transportation along those edges.)
 
+- Update how compounds are created. E.g. use
+```julia
+@variables t C(t) O(t)
+@compound CO2 ~ C + 2O
+```
+to create a compound species `CO2` that consists of `C` and 2 `O`.
+- Added documentation for chemistry related functionality (compound creation and reaction balancing).
 - Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage:
 ```julia
 using Catalyst
diff --git a/docs/pages.jl b/docs/pages.jl
index 9b21e12b10..51f3c02863 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -7,6 +7,7 @@ pages = Any["Home" => "index.md",
                                             "catalyst_functionality/constraint_equations.md",
                                             "catalyst_functionality/parametric_stoichiometry.md",
                                             "catalyst_functionality/network_analysis.md",
+                                            "catalyst_functionality/chemistry_related_functionality.md",
                                             "Model creation examples" => Any["catalyst_functionality/example_networks/basic_CRN_examples.md",
                                                                              "catalyst_functionality/example_networks/hodgkin_huxley_equation.md",
                                                                              "catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md"]],
diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index 2aa153c3e7..a6c48d03cc 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -280,6 +280,17 @@ Base.convert
 ModelingToolkit.structural_simplify
 ```
 
+## Chemistry-related functionalities
+Various functionalities primarily relevant to modelling of chemical systems (but potentially also in biology). 
+```@docs
+@compound
+@compounds
+iscompound
+components
+coefficients
+component_coefficients
+```
+
 ## Unit validation
 ```@docs
 validate(rx::Reaction; info::String = "")
diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md
new file mode 100644
index 0000000000..01c44094ef
--- /dev/null
+++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md
@@ -0,0 +1,138 @@
+# [Chemistry-related functionality](@id chemistry_functionality)
+
+While Catalyst has primarily been designed around the modelling of biological systems, reaction network models are also common across chemistry. This section describes two types of functionality, that while of general interest, should be especially useful in the modelling of chemical systems.
+- The `@compound` option, which enables the user to designate that a specific species is composed of certain subspecies.
+- The `balance_reaction` function, which enables the user to balance a reaction so the same number of components occur on both sides.
+
+## Modelling with compound species
+
+### Creating compound species programmatically
+We will first show how to create compound species through [programmatic model construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and their stoichiometries). In this example, we will create a CO₂ molecule, consisting of one C atom and two O atoms. First, we create species corresponding to the components:
+```@example chem1
+using Catalyst
+@variables t
+@species C(t) O(t) 
+```
+Next, we create the `CO2` compound species:
+```@example chem1
+@compound CO2 ~ C + 2O
+```
+Here, the compound is the first argument to the macro, followed by its component (with the left-hand and right-hand sides separated by a `~` sign). While non-compound species (such as `C` and `O`) have their independent variable (in this case `t`) designated, independent variables are generally not designated for compounds (these are instead directly inferred from their components). Components with non-unit stoichiometries have this value written before the component (generally, the rules for designating the components of a compound are identical to those of designating the substrates or products of a reaction). The created compound, `CO2`, is also a species, and can be used wherever e.g. `C` can be used:
+```@example chem1
+isspecies(CO2)
+```
+In its metadata, however, is stored information of its components, which can be retrieved using the `components` (returning a vector of its component species) and `coefficients` (returning a vector with each component's stoichiometry) functions:
+```@example chem1
+components(CO2)
+```
+```@example chem1
+coefficients(CO2)
+```
+Alternatively, we can retrieve the components and their stoichiometric coefficients as a single vector using the `component_coefficients` function:
+```@example chem1
+component_coefficients(CO2)
+```
+Finally, it is possible to check whether a species is a compound or not using the `iscompound` function:
+```@example chem1
+iscompound(CO2)
+```
+
+Compound components that are also compounds are allowed, e.g. we can create a carbonic acid compound (H₂CO₃) that consists of CO₂ and H₂O:
+```@example chem1
+@species H(t)
+@compound H2O ~ 2H + O
+@compound H2CO3 ~ CO2 + H2O
+```
+
+When multiple compounds are created, they can be created simultaneously using the `@compounds` macro, e.g. the previous code-block can be re-written as:
+```@example chem1
+@species H(t)
+@compounds begin
+    H2O ~ 2H + O
+    H2CO3 ~ CO2 + H2O
+end
+```
+
+### Creating compound species within the DSL
+It is also possible to declare species as compound species within the `@reaction_network` DSL, using the `@compounds` options:
+```@example chem1
+rn = @reaction_network begin
+    @species C(t) H(t) O(t)
+    @compounds begin
+        C2O ~ C + 2O
+        H2O ~ 2H + O
+        H2CO3 ~ CO2 + H2O
+    end
+    (k1,k2), H2O + CO2 <--> H2CO3
+end
+```
+When creating compound species using the DSL, it is important to note that *every component must be known to the system as a species, either by being declared using the `@species` or `@compound` options, or by appearing in a reaction*. E.g. the following is not valid
+```julia 
+rn = @reaction_network begin``
+    @compounds begin
+        C2O ~ C + 2O
+        H2O ~ 2H + O
+        H2CO3 ~ CO2 + H2O
+    end
+    (k1,k2), H2O+ CO2 <--> H2CO3
+end
+```
+as the components `C`, `H`, and `O` are not declared as a species anywhere. Please also note that only `@compounds` can be used as an option in the DSL, not `@compound`.
+
+### Designating metadata and default values for compounds
+Just like for normal species, it is possible to designate metadata and default values for compounds. Metadata is provided after the compound name, but separated from it by a `,`:
+```@example chem1
+@compound (CO2, [unit="mol"]) ~ C + 2O
+```
+Default values are designated using `=`, and provided directly after the compound name.:
+```@example chem1
+@compound (CO2 = 2.0) ~ C + 2O
+```
+If both default values and meta data are provided, the metadata is provided after the default value:
+```@example chem1
+@compound (CO2 = 2.0, [unit="mol"]) ~ C + 2O
+```
+In all of these cases, the side to the left of the `~` must be enclosed within `()`.
+
+### Compounds with multiple independent variables
+While we generally do not need to specify independent variables for compound, if the components (together) have more than one independent variable, this have to be done:
+```@example chem1
+@variables t s
+@species N(s) O(t) 
+@compound NO2(t,s) ~ N + 2O
+```
+Here, `NO2` depend both on a spatial independent variable (`s`) and a time one (`t`). This is required since, while multiple independent variables can be inferred, their internal order cannot (and must hence be provided by the user).
+
+## Balancing chemical reactions
+One use of defining a species as a compound is that they can be used to balance reactions so that the number of components are the same on both sides. Catalyst provides the `balance_reaction` function, which takes a reaction, and returns a balanced version. E.g. let us consider a reaction when carbon dioxide is formed from carbon and oxide `C + O --> CO2`. Here, `balance_reaction` enables us to find coefficients creating a balanced reaction (in this case, where the number of carbon and oxygen atoms are the same on both sides). To demonstrate, we first created the unbalanced reactions:
+```@example chem1
+rx = @reaction k, C + O --> $CO2
+```
+Here, the reaction rate (`k`) is not involved in the reaction balancing. We use interpolation for `CO2`, ensuring that the `CO2` used in the reaction is the same one we previously defined as a compound of `C` and `O`. Next, we call the `balance_reaction` function
+```@example chem1
+balance_reaction(rx)
+```
+which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that the reaction balancing problem may have several solutions. Typically, there is only a single solution (in which case this is the vector's only element). No, or an infinite number of, solutions is also possible depending on the given reaction.
+
+Let us consider a more elaborate example, the reaction between ammonia (NH₃) and oxygen (O₂) to form nitrogen monoxide (NO) and water (H₂O). Let us first create the components and the unbalanced reaction:
+```@example chem2
+using Catalyst # hide
+@variables t
+@species N(t) H(t) O(t) 
+@compounds begin
+    NH3 ~ N + 3H
+    O2 ~ 2O
+    NO ~ N + O
+    H2O ~ 2H + O
+end
+unbalanced_reaction = @reaction k, $NH3 + $O2 --> $NO + $H2O
+```
+We can now create a balanced version (where the amount of H, N, and O is the same on both sides):
+```@example chem2
+balanced_reaction = balance_reaction(unbalanced_reaction)[1]
+```
+
+Reactions declared as a part of a `ReactionSystem` (e.g. using the DSL) can be retrieved for balancing using the [`reactions`](@ref) function. Please note that balancing these will not mutate the `ReactionSystem`, but a new reaction system will need to be created using the balanced reactions.
+
+!!! note
+    Reaction balancing is currently not supported for reactions involving compounds of compounds.
\ No newline at end of file
diff --git a/docs/src/catalyst_functionality/constraint_equations.md b/docs/src/catalyst_functionality/constraint_equations.md
index 29dcbb2fc1..c0fd74bf99 100644
--- a/docs/src/catalyst_functionality/constraint_equations.md
+++ b/docs/src/catalyst_functionality/constraint_equations.md
@@ -17,7 +17,7 @@ $\lambda$. For now we'll assume the cell can grow indefinitely. We'll also keep
 track of one protein $P(t)$, which is produced at a rate proportional $V$ and
 can be degraded.
 
-## Coupling ODE constraints via extending a system
+## [Coupling ODE constraints via extending a system](@id constraint_equations_coupling_constratins)
 
 There are several ways we can create our Catalyst model with the two reactions
 and ODE for $V(t)$. One approach is to use compositional modeling, create
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index fa2b3e51c2..42398a7c95 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -102,8 +102,8 @@ export Graph, savegraph, complexgraph
 
 # for creating compounds
 include("chemistry_functionality.jl")
-export @compound
-export components, iscompound, coefficients
+export @compound, @compounds
+export iscompound, components, coefficients, component_coefficients
 export balance_reaction
 
 ### Extensions ###
diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl
index 96b0118f39..a839385ab3 100644
--- a/src/chemistry_functionality.jl
+++ b/src/chemistry_functionality.jl
@@ -1,3 +1,4 @@
+### Declares Compound Related Metadata ###
 struct CompoundSpecies end
 struct CompoundComponents end
 struct CompoundCoefficients end
@@ -6,93 +7,195 @@ Symbolics.option_to_metadata_type(::Val{:iscompound}) = CompoundSpecies
 Symbolics.option_to_metadata_type(::Val{:components}) = CompoundComponents
 Symbolics.option_to_metadata_type(::Val{:coefficients}) = CompoundCoefficients
 
-macro compound(species_expr, arr_expr...)
-    # Ensure the species name is a valid expression
-    if !(species_expr isa Expr && species_expr.head == :call)
-        error("Invalid species name in @compound macro")
+### Create @compound Macro(s) ###
+
+"""
+    @compound
+
+Macro that creates a compound species, which is composed of smaller component species.
+
+Example:
+```julia
+@variables t
+@species C(t) O(t)
+@compound CO2(t) ~ C + 2O
+```
+
+Notes: 
+- The component species must be defined before using the `@compound` macro.
+"""
+macro compound(expr)
+    make_compound(MacroTools.striplines(expr))
+end
+
+# Declares compound error messages:
+const COMPOUND_CREATION_ERROR_BASE = "Malformed input to @compound. Should use form like e.g. \"@compound CO2 ~ C + 2O\"."
+const COMPOUND_CREATION_ERROR_BAD_SEPARATOR = "Malformed input to @compound. Left-hand side (the compound) and the right-hand side (the components) should be separated by a \"~\" (e.g. \"@compound CO2 ~ C + 2O\"). If the left hand side contains metadata and/or default values, this should be enclosed by \"()\" (e.g. \"@compound (CO2 = 1.0, [output=true]) ~ C + 2O\")."
+const COMPOUND_CREATION_ERROR_DEPENDENT_VAR_REQUIRED = "When the components (collectively) have more than 1 independent variable, independent variables have to be specified for the compound (e.g. `@compound CO2(t,s) ~ C + 2O`). This is required since the @compound macro cannot infer the correct order of the independent variables."
+
+# Function managing the @compound macro.
+function make_compound(expr)
+    # Error checks.
+    (expr isa Expr) || error(COMPOUND_CREATION_ERROR_BASE)
+    ((expr.head == :call) && (expr.args[1] == :~) && (length(expr.args) == 3)) || error(COMPOUND_CREATION_ERROR_BAD_SEPARATOR)
+
+    # Loops through all components, add the component and the coefficients to the corresponding vectors
+    # Cannot extract directly using e.g. "getfield.(composition, :reactant)" because then 
+    # we get something like :([:C, :O]), rather than :([C, O]).
+    composition = Catalyst.recursive_find_reactants!(expr.args[3], 1, Vector{ReactantStruct}(undef, 0))
+    components = :([])                                      # Becomes something like :([C, O]).                                         
+    coefficients = :([])                                    # Becomes something like :([1, 2]). 
+    for comp in composition
+        push!(components.args, comp.reactant)
+        push!(coefficients.args, comp.stoichiometry)
     end
 
-    # Parse the species name to extract the species name and argument
-    species_name = species_expr.args[1]
-    species_arg = species_expr.args[2]
-
-    # Construct the expressions that define the species
-    species_expr = Expr(:macrocall, Symbol("@species"), LineNumberNode(0),
-                        Expr(:call, species_name, species_arg))
-
-    # Construct the expression to set the iscompound metadata
-    setmetadata_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
-                                                                       Catalyst.CompoundSpecies,
-                                                                       true))
-
-    # Ensure the expressions are evaluated in the correct scope by escaping them
-    escaped_species_expr = esc(species_expr)
-    escaped_setmetadata_expr = esc(setmetadata_expr)
-
-    # Construct the array from the remaining arguments
-    arr = Expr(:vect, (arr_expr)...)
-    coeffs = []
-    species = []
-
-    for expr in arr_expr
-        if isa(expr, Expr) && expr.head == :call && expr.args[1] == :*
-            push!(coeffs, expr.args[2])
-            push!(species, expr.args[3])
-        else
-            push!(coeffs, 1)
-            push!(species, expr)
-        end
+    # Extracts:
+    # - The compound species name (species_name, e.g. `:CO2`).
+    # - Any ivs attached to it (ivs, e.g. `[]` or `[t,x]`).
+    # - The expression which creates the compound (species_expr, e.g. `CO2 = 1.0, [metadata=true]`).
+    species_expr = expr.args[2]
+    species_name, ivs, _, _ = find_varinfo_in_declaration(expr.args[2])
+
+    # If no ivs were given, inserts `(..)` (e.g. turning `CO` to `CO(..)`).
+    isempty(ivs) && (species_expr = insert_independent_variable(species_expr, :(..)))      
+
+    # Expression which when evaluated gives a vector with all the ivs of the components.
+    ivs_get_expr = :(unique(reduce(vcat,[arguments(ModelingToolkit.unwrap(comp)) for comp in $components])))    
+    
+    # Creates the found expressions that will create the compound species.
+    # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in 
+    # the scope the users use the macro in (to e.g. detect already exiting species).     
+    # Creates something like (where `compound_ivs` and `component_ivs` evaluates to all the compound's and components' ivs):
+    #   `@species CO2(..)`
+    #   `isempty([])` && length(component_ivs) && error("When ...)
+    #   `CO2 = CO2(component_ivs..)` 
+    #   `issetequal(compound_ivs, component_ivs) || error("The ...)` 
+    #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)`
+    #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [C, O])`
+    #   `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [1, 2])`
+    species_declaration_expr = Expr(:escape, :(@species $species_expr)) 
+    multiple_ivs_error_check_expr = Expr(:escape, :($(isempty(ivs)) && (length($ivs_get_expr) > 1) && error($COMPOUND_CREATION_ERROR_DEPENDENT_VAR_REQUIRED)))      
+    iv_designation_expr = Expr(:escape, :($(isempty(ivs)) && ($species_name = $(species_name)($(ivs_get_expr)...)))) 
+    iv_check_expr = Expr(:escape, :(issetequal(arguments(ModelingToolkit.unwrap($species_name)), $ivs_get_expr) || error("The independent variable(S) provided to the compound ($(arguments(ModelingToolkit.unwrap($species_name)))), and those of its components ($($ivs_get_expr)))), are not identical.")))
+    compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true)))
+    components_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundComponents, $components)))
+    coefficients_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundCoefficients, $coefficients))) 
+
+    # Returns the rephrased expression.
+    return quote
+        $species_declaration_expr
+        $multiple_ivs_error_check_expr
+        $iv_designation_expr
+        $iv_check_expr
+        $compound_designation_expr
+        $components_designation_expr
+        $coefficients_designation_expr
     end
+end
 
-    coeffs_expr = Expr(:vect, coeffs...)
-    species_expr = Expr(:vect, species...)
+"""
+    @compounds
 
-    # Construct the expression to set the components metadata
-    setcomponents_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
-                                                                         Catalyst.CompoundComponents,
-                                                                         $species_expr))
+Macro that creates several compound species, which each is composed of smaller component species. Uses the same syntax as `@compound`, but with one compound species one each line.
 
-    # Ensure the expression is evaluated in the correct scope by escaping it
-    escaped_setcomponents_expr = esc(setcomponents_expr)
+Example:
+```julia
+@variables t
+@species C(t) H(t) O(t) 
+@compounds
+    CH4(t) = C + 4H
+    O2(t) = 2O
+    CO2(t) = C + 2O
+    H2O(t) = 2H + O
+end
+```
+
+Notes: 
+- The component species must be defined before using the `@compound` macro.
+"""
+macro compounds(expr)
+    make_compounds(MacroTools.striplines(expr))
+end
+
+# Function managing the @compound macro.
+function make_compounds(expr)
+    # Creates an empty block containing the output call.
+    compound_declarations = Expr(:block)
 
-    # Construct the expression to set the coefficients metadata
-    setcoefficients_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name),
-                                                                           Catalyst.CompoundCoefficients,
-                                                                           $coeffs_expr))
+    # For each compound in `expr`, creates the set of 7 compound creation lines (using `make_compound`).
+    # Next, loops through all 7*[Number of compounds] lines and add them to compound_declarations.
+    compound_calls = [Catalyst.make_compound(line) for line in expr.args] 
+    for compound_call in compound_calls, line in MacroTools.striplines(compound_call).args
+        push!(compound_declarations.args, line)
+    end
+
+    # The output of the macros should be a vector with the compounds (same as for e.g. "@species begin ... end", also require for things to work in the DSL).
+    # Creates an output vector, and loops through all compounds, adding them to it.
+    push!(compound_declarations.args, :($(Expr(:escape, :([])))))
+    compound_syms = :([])
+    for arg in expr.args
+        push!(compound_syms.args, find_varinfo_in_declaration(arg.args[2])[1])
+    end
+    push!(compound_declarations.args, :($(Expr(:escape, :($(compound_syms))))))
 
-    escaped_setcoefficients_expr = esc(setcoefficients_expr)
+    # The output needs to be converted to Vector{Num} (from  Vector{SymbolicUtils.BasicSymbolic{Real}}) to be consistent with e.g. @variables.
+    compound_declarations.args[end] = :([ModelingToolkit.wrap(cmp) for cmp in $(compound_declarations.args[end])])
 
-    # Return a block that contains the escaped expressions
-    return Expr(:block, escaped_species_expr, escaped_setmetadata_expr,
-                escaped_setcomponents_expr, escaped_setcoefficients_expr)
+    # Returns output that.
+    return compound_declarations
 end
 
-# Check if a species is a compound
+## Compound Getters ###
+
+"""
+    iscompound(s)
+
+Returns `true` if the input is a compound species (else false).
+"""
 iscompound(s::Num) = iscompound(MT.value(s))
 function iscompound(s)
     MT.getmetadata(s, CompoundSpecies, false)
 end
 
-coefficients(s::Num) = coefficients(MT.value(s))
-function coefficients(s)
-    MT.getmetadata(s, CompoundCoefficients)
-end
+"""
+    components(s)
 
+Returns a vector with a list of all the components of a compound species (created using e.g. the @compound macro).
+"""
 components(s::Num) = components(MT.value(s))
 function components(s)
     MT.getmetadata(s, CompoundComponents)
 end
 
+"""
+    coefficients(s)
+
+Returns a vector with a list of all the stoichiometric coefficients of the components of a compound species (created using e.g. the @compound macro).
+"""
+coefficients(s::Num) = coefficients(MT.value(s))
+function coefficients(s)
+    MT.getmetadata(s, CompoundCoefficients)
+end
+
+"""
+    component_coefficients(s)
+
+Returns a Vector{Pari{Symbol,Int64}}, listing a compounds species (created using e.g. the @compound macro) all the coefficients and their stoichiometric coefficients.
+"""
 component_coefficients(s::Num) = component_coefficients(MT.value(s))
 function component_coefficients(s)
     return [c => co for (c, co) in zip(components(s), coefficients(s))]
 end
 
 
-### Balancing Code
+### Reaction Balancing Functionality ###
+
+# Reaction balancing error.
 const COMPOUND_OF_COMPOUND_ERROR = ErrorException("Reaction balancing does not currently work for reactions involving compounds of compounds.")
 
-# note this does not correctly handle compounds of compounds currently
+# Note this does not correctly handle compounds of compounds currently.
+# Internal function used by "balance_reaction" (via "get_balanced_stoich").
 function create_matrix(reaction::Catalyst.Reaction)
     @unpack substrates, products = reaction
     unique_atoms = [] # Array to store unique atoms
@@ -143,6 +246,7 @@ function create_matrix(reaction::Catalyst.Reaction)
     return A
 end
 
+# Internal function used by "balance_reaction".
 function get_balanced_stoich(reaction::Reaction)
     # Create the reaction matrix A that is m atoms by n compounds
     A = create_matrix(reaction)
@@ -171,6 +275,52 @@ function get_balanced_stoich(reaction::Reaction)
     return stoichvecs
 end
 
+"""
+    balance_reaction(reaction::Reaction)
+
+Returns a vector of all possible stoichiometrically balanced `Reaction` objects for the given `Reaction`.
+
+Example:
+```julia
+@variables t
+@species Si(t) Cl(t) H(t) O(t)
+@compound SiCl4 ~ Si + 4Cl
+@compound H2O ~ 2H + O
+@compound H4SiO4 ~ 4H + Si + 4O
+@compound HCl ~ H + Cl
+rx = @reaction 1.0, SiCl4 + H2O --> H4SiO4 HCl
+balance_reaction(rx) # Exactly one solution.
+```
+
+```julia
+@variables t
+@species C(t) H(t) O(t)
+@compound CO ~ C + O
+@compound CO2 ~ C + 2O
+@compound H2 ~ 2H
+@compound CH4 ~ C + 4H
+@compound H2O ~ 2H + O
+rx = @reaction 1.0, CO + CO2 + H2--> CH4 H2O
+balance_reaction(rx) # Multiple solutions.
+```
+
+```julia
+@variables t
+@species Fe(t) S(t) O(t) H(t) N(t)
+@compound FeS2 ~ Fe + 2S
+@compound HNO3 ~ H + N + 3O
+@compound Fe2S3O12 ~ 2Fe + 3S + 12O
+@compound NO ~ N + O
+@compound H2SO4 ~ 2H + S + 4O
+rx = @reaction 1.0, FeS2 + HNO3 --> Fe2S3O12 NO + H2SO4
+brxs = balance_reaction(rx) # No solution.
+```
+
+Notes:
+- Balancing reactions that contain compounds of compounds is currently not supported.
+- A reaction may not always yield a single solution; it could have an infinite number of solutions or none at all. When there are multiple solutions, a vector of all possible `Reaction` objects is returned. However, substrates and products may be interchanged as we currently do not solve for a linear combination that maintains the set of substrates and products. 
+- If the reaction cannot be balanced, an empty `Reaction` vector is returned.
+"""
 function balance_reaction(reaction::Reaction)
     # Calculate the stoichiometric coefficients for the balanced reaction.
     stoichiometries = get_balanced_stoich(reaction)
diff --git a/src/expression_utils.jl b/src/expression_utils.jl
index 2778772e33..48c0cfb5ed 100644
--- a/src/expression_utils.jl
+++ b/src/expression_utils.jl
@@ -1,4 +1,4 @@
-#Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a  Symbol/Numerical).
+# Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a  Symbol/Numerical).
 function tup_leng(ex::ExprValues)
     (typeof(ex) == Expr && ex.head == :tuple) && (return length(ex.args))
     return 1
@@ -10,6 +10,83 @@ function get_tup_arg(ex::ExprValues, i::Int)
     return ex.args[i]
 end
 
+# In variable/species/parameters on the forms like:
+# X
+# X = 1.0
+# X, [metadata=true]
+# X = 1.0, [metadata=true]
+# X(t)
+# X(t) = 1.0
+# X(t), [metadata=true]
+# X(t) = 1.0, [metadata=true]
+# Finds the: Variable name (X), Independent variable name(s) ([t]), default value (2.0), and metadata (:([metadata=true])).
+# If a field does not exist (e.g. independent variable in `X, [metadata=true]`), gives nothing.
+# The independent variables are given as a vector (empty if none given).
+# Does not support e.g. "X [metadata=true]" (when metadata does not have a comma before).
+function find_varinfo_in_declaration(expr)
+    # Case: X
+    (expr isa Symbol) && (return expr, [], nothing, nothing)                          
+    # Case: X(t)         
+    (expr.head == :call) && (return expr.args[1], expr.args[2:end], nothing, nothing)                   
+    if expr.head == :(=)
+        # Case: X = 1.0
+        (expr.args[1] isa Symbol) && (return expr.args[1], [], expr.args[2], nothing)  
+        # Case: X(t) = 1.0        
+        (expr.args[1].head == :call) && (return expr.args[1].args[1], expr.args[1].args[2:end], expr.args[2].args[1], nothing) 
+    end
+    if expr.head == :tuple
+        # Case: X, [metadata=true]
+        (expr.args[1] isa Symbol) && (return expr.args[1], [], nothing, expr.args[2])          
+        # Case: X(t), [metadata=true]
+        (expr.args[1].head == :call) && (return expr.args[1].args[1], expr.args[1].args[2:end], nothing, expr.args[2]) 
+        if (expr.args[1].head == :(=)) 
+            # Case: X = 1.0, [metadata=true]
+            (expr.args[1].args[1] isa Symbol) && (return expr.args[1].args[1], [], expr.args[1].args[2], expr.args[2]) 
+            # Case: X(t) = 1.0, [metadata=true]
+            (expr.args[1].args[1].head == :call) && (return expr.args[1].args[1].args[1], expr.args[1].args[1].args[2:end], expr.args[1].args[2].args[1], expr.args[2]) 
+        end
+    end
+    error("Unable to detect the variable declared in expression: $expr.")
+end
+
+# Converts an expression of the forms:
+# X
+# X = 1.0
+# X, [metadata=true]
+# X = 1.0, [metadata=true]
+# To the form:
+# X(t)
+# X(t) = 1.0
+# X(t), [metadata=true]
+# X(t) = 1.0, [metadata=true]
+# (In this example the independent variable :t was inserted).
+# Here, the iv is a iv_expr, which can be anything, which is inserted
+function insert_independent_variable(expr_in, iv_expr)
+    # If expr is a symbol, just attach the iv. If not we have to create a new expr and mutate it. 
+    # Because Symbols (a possible input) cannot be mutated, this function cannot mutate the input
+    # (would have been easier if Expr input was guaranteed).
+    (expr_in isa Symbol) && (return Expr(:call, expr_in, iv_expr))
+    expr = deepcopy(expr_in)
+
+    # Loops through possible cases.
+    if expr.head == :(=) 
+        # Case: :(X = 1.0)
+        expr.args[1] = Expr(:call, expr.args[1], iv_expr)
+    elseif expr.head == :tuple
+        if expr.args[1] isa Symbol 
+            # Case: :(X, [metadata=true])
+            expr.args[1] = Expr(:call, expr.args[1], iv_expr)
+        elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol)
+            # Case: :(X = 1.0, [metadata=true])
+            expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], iv_expr)
+        end
+    end
+    return expr
+end
+
+
+### Old Stuff ###
+
 #This will be called whenever a function stored in funcdict is called.
 # function replace_names(expr, old_names, new_names)
 #     mapping = Dict(zip(old_names, new_names))
diff --git a/src/reaction_network.jl b/src/reaction_network.jl
index 6a5fb6871c..68d62796ba 100644
--- a/src/reaction_network.jl
+++ b/src/reaction_network.jl
@@ -123,7 +123,7 @@ const forbidden_variables_error = let
 end
 
 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -358,9 +358,12 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
                        option_lines))
 
+    # Reads compounds options.
+    compound_expr, compound_species = read_compound_options(options)
+
     # Parses reactions, species, and parameters.
     reactions = get_reactions(reaction_lines)
-    species_declared = extract_syms(options, :species)
+    species_declared = [extract_syms(options, :species); compound_species]
     parameters_declared = extract_syms(options, :parameters)
     variables = extract_syms(options, :variables)
 
@@ -399,6 +402,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps")
     vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars")
     sps, spssym = scalarize_macro(!isempty(species), sexprs, "specs")
+    comps, compssym = scalarize_macro(!isempty(compound_species), compound_expr, "comps")
     rxexprs = :(Catalyst.CatalystEqType[])
     for reaction in reactions
         push!(rxexprs.args, get_rxexprs(reaction))
@@ -410,8 +414,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         $ivexpr
         $vars
         $sps
+        $comps
 
-        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym),
+        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym, $compssym),
                                               $pssym; name = $name,
                                               spatial_ivs = $sivs)
     end
@@ -461,6 +466,20 @@ function esc_dollars!(ex)
     ex
 end
 
+# When compound species are declared using the "@compound begin ... end" option, get a list of the compound species, and also the expression that crates them.
+function read_compound_options(opts)
+    # If the compound option is used retrive a list of compound species (need to be added to teh reaction system's species), and the option that creates them (used to declare them as compounds at the end).
+    if haskey(opts, :compounds)
+        compound_expr = opts[:compounds]
+        # Find compound species names, and append the independent variable.
+        compound_species = [find_varinfo_in_declaration(arg.args[2])[1] for arg in compound_expr.args[3].args]
+    else  # If option is not used, return empty vectors and expressions.
+        compound_expr = :()
+        compound_species = Union{Symbol, Expr}[]
+    end
+    return compound_expr, compound_species
+end
+
 # When the user have used the @species (or @parameters) options, extract species (or
 # parameters) from its input.
 function extract_syms(opts, vartype::Symbol)
diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl
index be3afbb035..a93d6dfcbb 100644
--- a/test/miscellaneous_tests/compound_macro.jl
+++ b/test/miscellaneous_tests/compound_macro.jl
@@ -1,10 +1,99 @@
 using Catalyst, Test
 
-# Test base funcationality in two cases.
+### Tests Main Macro Creation Forms ### 
+let
+    @variables t
+    @species C(t) H(t) O(t) 
+    @parameters p1 p2
+
+    # Basic cases that should pass:
+    @compound H2O_1 ~ 2H + O
+    @compound (H2O_2, [output=true]) ~ 2H + O
+    @compound (H2O_3 = 1.5) ~ 2H + O
+    @compound (H2O_4 = 4, [output=true]) ~ 2H + O
+    @compound (H2O_5 = p1, [output=true]) ~ 2H + p2*O
+    @test iscompound(H2O_1)
+    @test iscompound(H2O_2)
+    @test iscompound(H2O_3)
+    @test iscompound(H2O_4)
+    @test iscompound(H2O_5)
+
+    # Other errors.
+    @test_throws Exception @eval @compound H2O = 2H + O
+    @test_throws Exception @eval @compound (H2O, [output=true]) = 2H + O
+    @test_throws Exception @eval @compound (H2O = 1.5) ~ 2H + O
+    @test_throws Exception @eval @compound (H2O = 4, [output=true]) ~ 2H + O
+    @test_throws Exception @eval @compound (H2O = p1, [output=true]) ~ 2H + p2*O
+
+    # Compounds created in block notation.
+    @compounds begin
+        CO2_1 ~ 2H + O
+    end
+    @compounds begin
+        (CO2_2, [output=true]) ~ 2H + O
+        (CO2_3 = 1.5) ~ 2H + O
+        (CO2_4 = 4, [output=true]) ~ 2H + O
+        (CO2_5 = p1, [output=true]) ~ 2H + p2*O
+    end
+    @test iscompound(CO2_1)
+    @test iscompound(CO2_2)
+    @test iscompound(CO2_3)
+    @test iscompound(CO2_4)
+    @test iscompound(CO2_5)
+
+    # Declares stuff in the DSL.
+    rn = @reaction_network begin
+        @species N(t) H(t) 
+        @parameters p1 p2
+        @compounds begin
+            NH3_1 ~ N + 3H
+            (NH3_2, [output=true]) ~ N + 3H
+            (NH3_3 = 1.5) ~ N + 3H
+            (NH3_4 = 4, [output=true]) ~ N + 3H
+            (NH3_5 = p1, [output=true]) ~ N + p2*H
+        end
+    end
+    @test iscompound(rn.NH3_1)
+    @test iscompound(rn.NH3_2)
+    @test iscompound(rn.NH3_3)
+    @test iscompound(rn.NH3_4)
+    @test iscompound(rn.NH3_5)
+end
+
+### Test Various Independent Variables ###
+let
+    @variables t x y z
+    @species C(t) H(x) N(x) O(t) P(t,x) S(x,y)
+
+    # Checks that wrong (or absent) independent variable produces errors.
+    @test_throws Exception @eval @compound CO2(t,x) ~ C + 2O
+    @test_throws Exception @eval @compound (NH4(s), [output=true]) ~ N + 4H
+    @test_throws Exception @eval @compound (H2O = 2.0) ~ 2H + O
+    @test_throws Exception @eval @compound PH4(x) ~ P + 4H
+    @test_throws Exception @eval @compound SO2(t,y) ~ S + 2O
+    
+    # Creates compounds.
+    @compound CO2 ~ C + 2O
+    @compound (NH4, [output=true]) ~ N + 4H
+    @compound (H2O(t,x) = 2.0) ~ 2H + O
+    @compound PH4(t,x) ~ P + 4H
+    @compound SO2(t,x,y) ~ S + 2O
+    
+    # Checks they have the correct independent variables.
+    @test issetequal(arguments(ModelingToolkit.unwrap(CO2)), [t])
+    @test issetequal(arguments(ModelingToolkit.unwrap(NH4)), [x])
+    @test issetequal(arguments(ModelingToolkit.unwrap(H2O)), [t, x])
+    @test issetequal(arguments(ModelingToolkit.unwrap(PH4)), [t, x])
+    @test issetequal(arguments(ModelingToolkit.unwrap(SO2)), [t, x, y])
+end
+
+### Other Minor Tests ###
+
+# Test base functionality in two cases.
 let
     @variables t
     @species C(t) H(t) O(t)
-    @compound C6H12O2(t) 6C 12H 2O
+    @compound C6H12O2 ~ 6C + 12H + 2O
 
     @test iscompound(C6H12O2)
     @test isspecies(C6H12O2)
@@ -21,7 +110,7 @@ end
 let
     @variables t
     @species O(t)
-    @compound O2(t) 2O
+    @compound O2 ~ 2O
 
     @test iscompound(O2)
     @test isspecies(O2)
@@ -37,19 +126,19 @@ end
 let
     @variables t
     @species C(t) H(t)
-    @test_throws Exception @compound C6H12O2(t) 6C 12H 2O    
+    @test_throws Exception @compound C6H12O2 ~ 6C + 12H + 2O    
 end
 let
     @variables t
-    @test_throws Exception @compound O2(t) 2O    
+    @test_throws Exception @compound O2 ~ 2O    
 end
 
 # Checks that nested components works as expected.
 let
     @variables t
     @species C(t) H(t) O(t)
-    @compound OH(t) 1O 1H
-    @compound C3H5OH3(t) 3C 5H 3OH
+    @compound OH ~ 1O + 1H
+    @compound C3H5OH3 ~ 3C + 5H + 3OH
 
     @test !iscompound(O)
     @test !iscompound(H)
@@ -71,22 +160,23 @@ let
     @variables t
     @species C(t) H(t) O(t)
     s = C
-    @compound C6H12O2_1(t) 6s 12H 2O
-    @compound C6H12O2_2(t) 6C 12H 2O
+    @compound C6H12O2_1 ~ 6s + 12H + 2O
+    @compound C6H12O2_2 ~ 6C + 12H + 2O
 
     @test iscompound(C6H12O2_1)
     @test iscompound(C6H12O2_2)
 
     @test isequal(components(C6H12O2_1), components(C6H12O2_2))
     @test isequal(coefficients(C6H12O2_1), coefficients(C6H12O2_2))
+    @test isequal(component_coefficients(C6H12O2_1), component_coefficients(C6H12O2_2))
 end
 
 let
     @variables t
     @species C(t) H(t)
-    @compound Cyclopentadiene(t) 5C 6H
+    @compound Cyclopentadiene ~ 5C + 6H
     C5H6 = Cyclopentadiene
-    @compound C10H12(t) 2C5H6
+    @compound C10H12 ~ 2C5H6
 
     @test iscompound(C10H12)
     @test iscompound(components(C10H12)[1])
@@ -101,14 +191,23 @@ let
     @species H(t)
 
     alpha = 2
-    @compound H2_1(t) alpha*H
-    @compound H2_2(t) 2H
+    h = H
+    @compound H2_1 ~ 2*H
+    @compound H2_2 ~ alpha*H
+    @compound H2_3 ~ 2*h
+    @compound H2_4 ~ alpha*H
 
     @test iscompound(H2_1)
     @test iscompound(H2_2)
+    @test iscompound(H2_2)
+    @test iscompound(H2_4)
 
     @test isequal(components(H2_1),components(H2_2))
+    @test isequal(components(H2_2),components(H2_3))
+    @test isequal(components(H2_3),components(H2_4))
     @test isequal(coefficients(H2_1),coefficients(H2_2))
+    @test isequal(coefficients(H2_2),coefficients(H2_3))
+    @test isequal(coefficients(H2_3),coefficients(H2_4))
 end
 
 let
@@ -116,8 +215,8 @@ let
     @parameters alpha = 2
     @species H(t)
 
-    @compound H2_1(t) alpha*H
-    @compound H2_2(t) 2H
+    @compound H2_1 ~ alpha*H
+    @compound H2_2 ~ 2H
 
     @test iscompound(H2_1)
     @test iscompound(H2_2)
@@ -130,13 +229,133 @@ let
     @variables t
     @species A(t)
     B = A
-    @compound A2(t) 2A
-    @compound B2(t) 2B
+    @compound A2 ~ 2A
+    @compound B2 ~ 2B
 
     @test iscompound(A2)
     @test iscompound(B2)
 
     @test isequal(components(A2),components(B2))
     @test isequal(coefficients(A2), coefficients(B2))
+    @test isequal(component_coefficients(A2), component_coefficients(B2))
+end
+
+###  Check @compounds Macro ###
+
+# Basic syntax.
+let 
+    @variables t
+    @species C(t) H(t) O(t)
+    @compound OH ~ 1O + 1H
+    @compound C3H5OH3 ~ 3C + 5H + 3OH
+
+    @compounds begin
+        OH_alt ~ 1O + 1H
+        C3H5OH3_alt ~ 3C + 5H + 3OH
+    end
+
+    @test iscompound(OH_alt)
+    @test iscompound(C3H5OH3_alt)
+
+    @test isequal(components(OH),components(OH_alt))
+    @test isequal(coefficients(OH), coefficients(OH_alt))
+    @test isequal(component_coefficients(OH), component_coefficients(OH_alt))
+    @test isequal(components(C3H5OH3),components(C3H5OH3_alt))
+    @test isequal(coefficients(C3H5OH3), coefficients(C3H5OH3_alt))
+    @test isequal(component_coefficients(C3H5OH3), component_coefficients(C3H5OH3_alt))
 end
 
+# Interpolation
+let 
+    @variables t
+    @species s1(t) s2(t) s3(t)
+    s2_alt = s2
+    s3_alt = s3
+
+    @compounds begin
+        comp ~ s1 + s2 + 4s3
+        comp_alt ~ s1 + s2_alt + 4s3_alt
+    end
+
+    @test iscompound(comp)
+    @test iscompound(comp_alt)
+
+    @test isequal(components(comp),components(comp_alt))
+    @test isequal(coefficients(comp), coefficients(comp_alt))
+    @test isequal(component_coefficients(comp), component_coefficients(comp_alt))
+end
+
+### Compounds in DSL ###
+
+# Checks with a single compound.
+# Checks using @unpack.
+# Check where compounds and components does not occur in reactions.
+let 
+    rn = @reaction_network begin
+        @species C(t) O(t)
+        @compounds begin
+            CO2 ~ C + 2O
+        end
+    end
+    @unpack C, O, CO2 = rn
+    
+    @test length(species(rn)) == 3
+    @test iscompound(CO2)
+    @test isequal([C, O], components(CO2))
+    @test isequal([1, 2], coefficients(CO2))
+    @test isequal([C => 1, O => 2], component_coefficients(CO2)) 
+end
+
+# Test using multiple compounds.
+# Test using rn. notation to fetch species.
+let 
+    rn = complete(@reaction_network begin
+        @species C(t) O(t) H(t)
+        @compounds begin
+            CH4 ~ C + 4H
+            O2 ~ 2O
+            CO2 ~ C + 2O
+            H2O ~ 2H + O
+        end
+        k, CH4 + O2 --> CO2 + H2O
+    end)
+    species(rn)
+    
+    @test length(species(rn)) == 7
+    @test isequal([rn.C, rn.H], components(rn.CH4))
+    @test isequal([1, 4], coefficients(rn.CH4))
+    @test isequal([rn.C => 1, rn.H => 4], component_coefficients(rn.CH4)) 
+    @test isequal([rn.O], components(rn.O2))
+    @test isequal([2], coefficients(rn.O2))
+    @test isequal([rn.O => 2], component_coefficients(rn.O2)) 
+    @test isequal([rn.C, rn.O], components(rn.CO2))
+    @test isequal([1, 2], coefficients(rn.CO2))
+    @test isequal([rn.C => 1, rn.O => 2], component_coefficients(rn.CO2)) 
+    @test isequal([rn.H, rn.O], components(rn.H2O))
+    @test isequal([2, 1], coefficients(rn.H2O))
+    @test isequal([rn.H => 2, rn.O => 1], component_coefficients(rn.H2O)) 
+end
+
+# Tests using compounds of compounds.
+# Tests where species are part of reactions and not declared using "@species".
+let
+    rn = @reaction_network begin
+        @compounds begin
+            SO2 ~ S + 2O
+            S2O4 ~ 2SO2
+        end
+        dS, S --> 0
+        dO, O --> 0
+    end
+    species(rn)
+    @unpack S, O, SO2, S2O4 = rn
+    
+    @test length(species(rn)) == 4
+    
+    @test isequal([S, O], components(SO2))
+    @test isequal([1, 2], coefficients(SO2))
+    @test isequal([S => 1, O => 2], component_coefficients(SO2)) 
+    @test isequal([SO2], components(S2O4))
+    @test isequal([2], coefficients(S2O4))
+    @test isequal([SO2 => 2], component_coefficients(S2O4)) 
+end
\ No newline at end of file
diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl
index b858cf7aa9..2a5194a080 100644
--- a/test/miscellaneous_tests/reaction_balancing.jl
+++ b/test/miscellaneous_tests/reaction_balancing.jl
@@ -5,9 +5,11 @@ let
     @variables t
     @parameters k
     @species H(t) O(t)
-    @compound H2(t) 2H
-    @compound O2(t) 2O
-    @compound H2O(t) 2H 1O
+    @compounds begin
+        H2 ~ 2H
+        O2 ~ 2O
+        H2O ~ 2H + 1O
+    end
 
     rx = Reaction(k,[H2,O2],[H2O])
 
@@ -18,16 +20,20 @@ let
     brxs = balance_reaction(rx)
     @test length(brxs) == 1
     @test isequal(balanced_rx, first(brxs))
+
+    rx_macro = @reaction k, $H2 + $O2 --> $H2O
+    brxs_macro = balance_reaction(rx_macro)
+    @test isequal(brxs, brxs_macro)
 end
 
 let
     @variables t
     @parameters k
     @species C(t) H(t) O(t)
-    @compound O2(t) 2O
-    @compound CO2(t) 1C 2O
-    @compound H2O(t) 2H 1O
-    @compound C6H12O6(t) 6C 12H 6O
+    @compound O2 ~ 2O
+    @compound CO2 ~ 1C + 2O
+    @compound H2O ~ 2H + 1O
+    @compound C6H12O6 ~ 6C + 12H + 6O
 
     rx = Reaction(k,[CO2,H2O],[C6H12O6,O2])
 
@@ -38,13 +44,17 @@ let
     brxs = balance_reaction(rx)
     @test length(brxs) == 1
     @test isequal(balanced_rx, first(brxs))
+
+    rx_macro = @reaction k, $CO2 + $H2O --> $C6H12O6 + $O2
+    brxs_macro = balance_reaction(rx_macro)
+    @test isequal(brxs, brxs_macro)
 end
 
 # @reaction k, H2O --> H2O
 let
     @variables t
     @species H(t) O(t)
-    @compound H2O(t) 2H O
+    @compound H2O ~ 2H + O
 
     rx = Reaction(1.0, [H2O], [H2O], [2], [2])
 
@@ -58,7 +68,7 @@ end
 let
     @variables t
     @species H(t) O(t)
-    @compound H2O(t) 2H O
+    @compound H2O ~ 2H + O
 
     rx = Reaction(1.0, [H, O], [H2O], [23, 1], [7])
 
@@ -72,10 +82,12 @@ end
 let
     @variables t
     @species H(t) O(t) C(t)
-    @compound CH4(t) C 4H
-    @compound O2(t) 2O
-    @compound CO2(t) C 2O
-    @compound H2O(t) 2H O
+    @compounds begin
+        CH4 ~ C + 4H
+        O2 ~ 2O
+        CO2 ~ C + 2O
+        H2O ~ 2H + O
+    end
 
     rx = Reaction(1.0, [CH4, O2], [CO2, H2O])
 
@@ -89,9 +101,9 @@ end
 let
     @variables t
     @species H(t) N(t)
-    @compound N2(t) 2N
-    @compound H2(t) 2H
-    @compound NH3(t) N 3H
+    @compound N2 ~ 2N
+    @compound H2 ~ 2H
+    @compound NH3 ~ N + 3H
 
     rx = Reaction(1.0, [N2, H2], [NH3])
 
@@ -105,10 +117,10 @@ end
 let
     @variables t
     @species C(t) H(t) O(t)
-    @compound C2H5OH(t) 2C 6H O
-    @compound CH3COOH(t) 2C 4H 2O
-    @compound C4H8O2(t) 4C 8H 2O
-    @compound H2O(t) 2H O
+    @compound C2H5OH ~ 2C + 6H + O
+    @compound CH3COOH ~ 2C + 4H + 2O
+    @compound C4H8O2 ~ 4C + 8H + 2O
+    @compound H2O ~ 2H + O
 
     rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O])
 
@@ -122,9 +134,9 @@ end
 let
     @variables t
     @species Ca(t) P(t) O(t)
-    @compound Ca3PO42(t) 3Ca 2P 8O
-    @compound CaO(t) Ca O
-    @compound P4O10(t) 4P 10O
+    @compound Ca3PO42 ~ 3Ca + 2P + 8O
+    @compound CaO ~ Ca + O
+    @compound P4O10 ~ 4P + 10O
 
     rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10])
 
@@ -138,9 +150,9 @@ end
 let
     @variables t
     @species Fe(t) O(t) H(t)
-    @compound O2(t) 2O
-    @compound H2O(t) 2H O
-    @compound FeOH3(t) Fe 3H 3O
+    @compound O2 ~ 2O
+    @compound H2O ~ 2H + O
+    @compound FeOH3 ~ Fe + 3H + 3O
 
     rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3])
 
@@ -154,11 +166,11 @@ end
 let
     @variables t
     @species Na(t) O(t) H(t) S(t)
-    @compound SO4(t) S 4O
-    @compound NaOH(t) Na O H
-    @compound H2SO4(t) 2H 1S 4O
-    @compound Na2SO4(t) 2Na 1S 4O
-    @compound H2O(t) 2H O
+    @compound SO4 ~ S + 4O
+    @compound NaOH ~ Na + O + H
+    @compound H2SO4 ~ 2H + 1S + 4O
+    @compound Na2SO4 ~ 2Na + 1S + 4O
+    @compound H2O ~ 2H + O
 
     rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O])
 
@@ -172,8 +184,8 @@ end
 let
     @variables t
     @species N(t) O(t)
-    @compound NO2(t) N 2O
-    @compound N2O4(t) 2N 4O
+    @compound NO2 ~ N + 2O
+    @compound N2O4 ~ 2N + 4O
 
     rx = Reaction(1.0, [NO2], [N2O4])
 
@@ -187,11 +199,11 @@ end
 let
     @variables t
     @species C(t) H(t) O(t) Ca(t) Cl(t)
-    @compound H2O(t) 2H 1O
-    @compound CO2(t) 1C 2O
-    @compound CaCO3(t) 1Ca 1C 3O
-    @compound HCl(t) 1H 1Cl
-    @compound CaCl2(t) 1Ca 2Cl
+    @compound H2O ~ 2H + 1O
+    @compound CO2 ~ 1C + 2O
+    @compound CaCO3 ~ 1Ca + 1C + 3O
+    @compound HCl ~ 1H + 1Cl
+    @compound CaCl2 ~ 1Ca + 2Cl
 
     rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O])
     balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1])
@@ -204,10 +216,10 @@ end
 let
     @variables t
     @species Si(t) Cl(t) H(t) O(t)
-    @compound SiCl4(t) 1Si 4Cl
-    @compound H2O(t) 2H O
-    @compound H4SiO4(t) 4H Si 4O
-    @compound HCl(t) H Cl
+    @compound SiCl4 ~ 1Si + 4Cl
+    @compound H2O ~ 2H + O
+    @compound H4SiO4 ~ 4H + Si + 4O
+    @compound HCl ~ H + Cl
 
     rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl])
     balanced_rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl], [1,4], [1,4])
@@ -220,9 +232,9 @@ end
 let
     @variables t
     @species Al(t) Cl(t) H(t)
-    @compound HCl(t) H Cl
-    @compound AlCl3(t) Al 3Cl
-    @compound H2(t) 2H
+    @compound HCl ~ H + Cl
+    @compound AlCl3 ~ Al + 3Cl
+    @compound H2 ~ 2H
 
     rx = Reaction(1.0,[Al,HCl],[AlCl3,H2])
     balanced_rx = Reaction(1.0,[Al,HCl],[AlCl3,H2],[2,6], [2,3])
@@ -235,11 +247,11 @@ end
 let
     @variables t
     @species Na(t) C(t) O(t) H(t) Cl(t)
-    @compound Na2CO3(t) 2Na C 3O
-    @compound HCl(t) H Cl
-    @compound NaCl(t) Na Cl
-    @compound H2O(t) 2H O
-    @compound CO2(t) C 2O
+    @compound Na2CO3 ~ 2Na + C + 3O
+    @compound HCl ~ H + Cl
+    @compound NaCl ~ Na + Cl
+    @compound H2O ~ 2H + O
+    @compound CO2 ~ C + 2O
 
     rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2])
     balanced_rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2], [1,2], [2,1,1])
@@ -252,10 +264,10 @@ end
 let
     @variables t
     @species C(t) H(t) O(t)
-    @compound C7H6O2(t) 7C 6H 2O
-    @compound O2(t) 2O
-    @compound CO2(t) C 2O
-    @compound H2O(t) 2H O
+    @compound C7H6O2 ~ 7C + 6H + 2O
+    @compound O2 ~ 2O
+    @compound CO2 ~ C + 2O
+    @compound H2O ~ 2H + O
 
     rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O])
     balanced_rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O], [2,15], [14,6])
@@ -268,10 +280,10 @@ end
 let
     @variables t
     @species Fe(t) S(t) O(t) H(t) K(t)
-    @compound Fe2S3O12(t) 2Fe 3S 12O
-    @compound KOH(t) K O H
-    @compound K2SO4(t) 2K S 4O
-    @compound FeO3H3(t) Fe 3O 3H
+    @compound Fe2S3O12 ~ 2Fe + 3S + 12O
+    @compound KOH ~ K + O + H
+    @compound K2SO4 ~ 2K + S + 4O
+    @compound FeO3H3 ~ Fe + 3O + 3H
 
     rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3]) #5x4 matrix
     balanced_rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3], [1,6], [3,2])
@@ -284,10 +296,10 @@ end
 let
     @variables t
     @species Ca(t) P(t) O(t) Si(t)
-    @compound Ca3P2O8(t) 3Ca 2P 8O
-    @compound SiO2(t) Si 2O
-    @compound P4O10(t) 4P 10O
-    @compound CaSiO3(t) Ca Si 3O
+    @compound Ca3P2O8 ~ 3Ca + 2P + 8O
+    @compound SiO2 ~ Si + 2O
+    @compound P4O10 ~ 4P + 10O
+    @compound CaSiO3 ~ Ca + Si + 3O
 
     rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3]) #5x4 matrix
     balanced_rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3], [2,6] , [1,6])
@@ -300,9 +312,9 @@ end
 let
     @variables t
     @species K(t) Cl(t) O(t)
-    @compound KClO3(t) K Cl 3O
-    @compound KClO4(t) K Cl  4O
-    @compound KCl(t) K Cl
+    @compound KClO3 ~ K + Cl + 3O
+    @compound KClO4 ~ K + Cl + 4O
+    @compound KCl ~ K + Cl
 
     rx = Reaction(1.0,[KClO3],[KClO4,KCl])
     balanced_rx = Reaction(1.0,[KClO3],[KClO4,KCl], [4], [3,1])
@@ -314,11 +326,11 @@ end
 # @reaction k, Al2(SO4)3 + 3Ca(OH)2 → 2Al(OH)3 + 3CaSO4
 let
     @variables t
-    @species Al(t) S(t) O(t) Ca(t) O(t) (H)
-    @compound Al2S3O12(t) 2Al 3S 12O
-    @compound CaO2H2(t) Ca 2O 2H
-    @compound AlO3H3(t) Al 3O 3H
-    @compound CaSO4(t) Ca S 4O
+    @species Al(t) S(t) O(t) Ca(t) O(t) H(t)
+    @compound Al2S3O12 ~ 2Al + 3S + 12O
+    @compound CaO2H2 ~ Ca + 2O + 2H
+    @compound AlO3H3 ~ Al + 3O + 3H
+    @compound CaSO4 ~ Ca + S + 4O
 
     rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4])
     balanced_rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4], [1,3], [2,3])
@@ -331,11 +343,11 @@ end
 let
     @variables t
     @species H(t) S(t) O(t) I(t)
-    @compound H2SO4(t) 2H S 4O
-    @compound HI(t) H I
-    @compound H2S(t) 2H S
-    @compound I2(t) 2I
-    @compound H2O(t) 2H O
+    @compound H2SO4 ~ 2H + S + 4O
+    @compound HI ~ H + I
+    @compound H2S ~ 2H + S
+    @compound I2 ~ 2I
+    @compound H2O ~ 2H + O
 
     rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O])
     balanced_rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O], [1,8], [1,4,4])
@@ -348,10 +360,10 @@ end
 let
     @variables t
     @species C(t) H(t) O(t)
-    @compound C2H4(t) 2C 4H
-    @compound O2(t) 2O
-    @compound CO2(t) C 2O
-    @compound H2O(t) 2H O
+    @compound C2H4 ~ 2C + 4H
+    @compound O2 ~ 2O
+    @compound CO2 ~ C + 2O
+    @compound H2O ~ 2H + O
 
     rx = Reaction(1.0,[C2H4,O2],[CO2,H2O])
     balanced_rx = Reaction(1.0,[C2H4,O2],[CO2,H2O],[1,3],[2,2])
@@ -364,11 +376,11 @@ end
 let
     @variables t
     @species C(t) H(t) O(t)
-    @compound CO(t) C O
-    @compound CO2(t) C 2O
-    @compound H2(t) 2H
-    @compound CH4(t) C 4H
-    @compound H2O(t) 2H O
+    @compound CO ~ C + O
+    @compound CO2 ~ C + 2O
+    @compound H2 ~ 2H
+    @compound CH4 ~ C + 4H
+    @compound H2O ~ 2H + O
 
     rx = Reaction(1.0,[CO,CO2,H2],[CH4,H2O])
     brxs = balance_reaction(rx)
@@ -381,11 +393,11 @@ let
     @variables t
     @species Fe(t) S(t) O(t) H(t) N(t)
 
-    @compound FeS2(t) Fe 2S
-    @compound HNO3(t) H N 3O
-    @compound Fe2S3O12(t) 2Fe 3S 12O
-    @compound NO(t) N O
-    @compound H2SO4(t) 2H S 4O
+    @compound FeS2 ~ Fe + 2S
+    @compound HNO3 ~ H + N + 3O
+    @compound Fe2S3O12 ~ 2Fe + 3S + 12O
+    @compound NO ~ N + O
+    @compound H2SO4 ~ 2H + S + 4O
 
     rx = Reaction(1.0,[FeS2,HNO3],[Fe2S3O12,NO,H2SO4])
     brxs = balance_reaction(rx)
@@ -397,10 +409,32 @@ end
 let
     @variables t
     @species C(t) H(t) O(t)
-    @compound CO(t) C O
-    @compound H2(t) 2H
-    @compound COH2(t) CO H2
+    @compound CO ~ C + O
+    @compound H2 ~ 2H
+    @compound COH2 ~ CO + H2
 
     rx = Reaction(1.0, [CO, H2], [COH2])
     @test_throws Catalyst.COMPOUND_OF_COMPOUND_ERROR balance_reaction(rx)
 end
+
+# Checks that balancing work for a reaction from a reaction_network.
+let
+    rn = complete(@reaction_network begin
+        @species C(t) H(t) O(t)
+        @compounds begin
+            O2 ~ 2O
+            CO2 ~ 1C + 2O
+            H2O ~ 2H + 1O
+            C6H12O6 ~ 6C + 12H + 6O
+        end
+        k, CO2 + H2O --> C6H12O6 + O2
+    end)
+    
+    brxs = balance_reaction(reactions(rn)[1])[1]
+    
+    @test isequal(rn.k, brxs.rate)
+    @test isequal([rn.CO2,  rn.H2O], brxs.substrates)
+    @test isequal([rn.C6H12O6,  rn.O2], brxs.products)
+    @test isequal([6, 6], brxs.substoich)
+    @test isequal([1,  6], brxs.prodstoich)    
+end
\ No newline at end of file