diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 81208e1bf7..abab4f66ac 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -260,25 +260,19 @@ end
 Construct a directed simple graph where nodes correspond to reaction complexes and directed
 edges to reactions converting between two complexes.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-complexes,incidencemat = reactioncomplexes(sir)
 incidencematgraph(sir)
 ```
 """
 function incidencematgraph(rn::ReactionSystem)
     nps = get_networkproperties(rn)
     if Graphs.nv(nps.incidencegraph) == 0
-        isempty(nps.incidencemat) &&
-            error("Please call reactioncomplexes(rn) first to construct the incidence matrix.")
+        isempty(nps.incidencemat) && reactioncomplexes(rn)
         nps.incidencegraph = incidencematgraph(nps.incidencemat)
     end
     nps.incidencegraph
@@ -329,17 +323,12 @@ Given the incidence graph of a reaction network, return a vector of the
 connected components of the graph (i.e. sub-groups of reaction complexes that
 are connected in the incidence graph).
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-complexes,incidencemat = reactioncomplexes(sir)
 linkageclasses(sir)
 ```
 gives
@@ -371,17 +360,12 @@ Here the deficiency, ``\delta``, of a network with ``n`` reaction complexes,
 \delta = n - \ell - s
 ```
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-rcs,incidencemat = reactioncomplexes(sir)
 δ = deficiency(sir)
 ```
 """
@@ -419,17 +403,12 @@ end
 
 Find subnetworks corresponding to each linkage class of the reaction network.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-complexes,incidencemat = reactioncomplexes(sir)
 subnetworks(sir)
 ```
 """
@@ -456,17 +435,12 @@ end
 
 Calculates the deficiency of each sub-reaction network within `network`.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-rcs,incidencemat = reactioncomplexes(sir)
 linkage_deficiencies = linkagedeficiencies(sir)
 ```
 """
@@ -487,17 +461,12 @@ end
 
 Given a reaction network, returns if the network is reversible or not.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-rcs,incidencemat = reactioncomplexes(sir)
 isreversible(sir)
 ```
 """
@@ -511,30 +480,27 @@ end
 
 Determine if the reaction network with the given subnetworks is weakly reversible or not.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-rcs,incidencemat = reactioncomplexes(sir)
 subnets = subnetworks(rn)
 isweaklyreversible(rn, subnets)
 ```
 """
 function isweaklyreversible(rn::ReactionSystem, subnets)
-    im = get_networkproperties(rn).incidencemat
-    isempty(im) &&
-        error("Error, please call reactioncomplexes(rn::ReactionSystem) to ensure the incidence matrix has been cached.")
-    sparseig = issparse(im)
+    nps = get_networkproperties(rn)
+    isempty(nps.incidencemat) && reactioncomplexes(rn)
+    sparseig = issparse(nps.incidencemat)
+
     for subnet in subnets
-        nps = get_networkproperties(subnet)
-        isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig)
+        subnps = get_networkproperties(subnet)
+        isempty(subnps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig)
     end
+
+    # A network is weakly reversible if all of its subnetworks are strongly connected
     all(Graphs.is_strongly_connected ∘ incidencematgraph, subnets)
 end
 
diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl
index 5507a9f655..db5b9f5893 100644
--- a/test/network_analysis/network_properties.jl
+++ b/test/network_analysis/network_properties.jl
@@ -325,4 +325,3 @@ let
     rates = Dict(zip(parameters(rn), k))
     @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
-
diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl
index 514dbc136f..538ce0f174 100644
--- a/test/reactionsystem_core/reactionsystem.jl
+++ b/test/reactionsystem_core/reactionsystem.jl
@@ -816,3 +816,42 @@ end
 let
     @test_nowarn Catalyst.reactionsystem_uptodate_check()
 end
+
+# Test that functions using the incidence matrix properly cache it
+let
+    rn = @reaction_network begin
+        k1, A --> B
+        k2, B --> C
+        k3, C --> A
+    end
+
+    nps = Catalyst.get_networkproperties(rn)
+    @test isempty(nps.incidencemat) == true
+
+    img = incidencematgraph(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    lcs = linkageclasses(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    sns = subnetworks(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    δ = deficiency(rn)
+    @test size(nps.incidencemat) == (3,3)
+    
+    Catalyst.reset!(nps)
+    δ_l = linkagedeficiencies(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    rev = isreversible(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    weakrev = isweaklyreversible(rn, sns)
+    @test size(nps.incidencemat) == (3,3)
+end