diff --git a/docs/Project.toml b/docs/Project.toml
index 39cdb81ad0..08c3ed624d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -35,7 +35,6 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
-#StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
@@ -58,7 +57,7 @@ JumpProcesses = "9.13.2"
 Latexify = "0.16.5"
 LinearSolve = "2.30"
 ModelingToolkit = "9.32"
-NonlinearSolve = "3.12"
+NonlinearSolve = "3.12, 4"
 Optim = "1.9"
 Optimization = "4"
 OptimizationBBO = "0.4"
@@ -74,5 +73,4 @@ SpecialFunctions = "2.4"
 StaticArrays = "1.9"
 SteadyStateDiffEq = "2.2"
 StochasticDiffEq = "6.65"
-#StructuralIdentifiability = "0.5.8"
 Symbolics = "5.30.1, 6"
diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 3e42e064c1..93bebe1888 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -398,6 +398,19 @@ function isterminal(lc::Vector, rn::ReactionSystem)
     true
 end
 
+function isforestlike(rn::ReactionSystem)
+    subnets = subnetworks(rn)
+    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)
+    end
+    all(Graphs.is_tree ∘ SimpleGraph ∘ incidencematgraph, subnets)
+end
+
 @doc raw"""
     deficiency(rn::ReactionSystem)
 
@@ -724,6 +737,86 @@ function conservationlaw_errorcheck(rs, pre_varmap)
         error("The system has conservation laws but initial conditions were not provided for some species.")
 end
 
+"""
+    isdetailedbalanced(rs::ReactionSystem, parametermap; reltol=1e-9, abstol)
+
+Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium
+solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+"""
+function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, reltol=1e-9)
+    if length(parametermap) != numparams(rs)
+        error("Incorrect number of parameters specified.")
+    elseif !isreversible(rs)
+        return false
+    elseif !all(r -> ismassaction(r, rs), reactions(rs))
+        error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.")
+    end
+
+    isforestlike(rs) && deficiency(rs) == 0 && return true
+
+    pmap = symmap_to_varmap(rs, parametermap)
+    pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)
+
+    # Construct reaction-complex graph 
+    complexes, D = reactioncomplexes(rs)
+    img = incidencematgraph(rs)
+    undir_img = SimpleGraph(incidencematgraph(rs))
+    K = ratematrix(rs, pmap)
+
+    spanning_forest = Graphs.kruskal_mst(undir_img)
+    outofforest_edges = setdiff(collect(edges(undir_img)), spanning_forest)
+
+    # Independent Cycle Conditions: for any cycle we create by adding in an out-of-forest reaction, the product of forward reaction rates over the cycle must equal the product of reverse reaction rates over the cycle.  
+    for edge in outofforest_edges
+        g = SimpleGraph([spanning_forest..., edge])
+        ic = Graphs.cycle_basis(g)[1]
+        fwd = prod([K[ic[r], ic[r + 1]] for r in 1:(length(ic) - 1)]) * K[ic[end], ic[1]]
+        rev = prod([K[ic[r + 1], ic[r]] for r in 1:(length(ic) - 1)]) * K[ic[1], ic[end]]
+        isapprox(fwd, rev; atol = abstol, rtol = reltol) ? continue : return false
+    end
+
+    # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (we will take the one given by default from kruskal_mst).  
+
+    if deficiency(rs) > 0
+        rxn_idxs = [edgeindex(D, Graphs.src(e), Graphs.dst(e)) for e in spanning_forest]
+        S_F = netstoichmat(rs)[:, rxn_idxs]
+        sols = positive_nullspace(S_F)
+
+        for i in 1:size(sols, 2)
+            α = sols[:, i]
+            fwd = prod([K[Graphs.src(e), Graphs.dst(e)]^α[i]
+                        for (e, i) in zip(spanning_forest, 1:length(α))])
+            rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i]
+                        for (e, i) in zip(spanning_forest, 1:length(α))])
+            isapprox(fwd, rev; atol = abstol, rtol = reltol) ? continue : return false
+        end
+    end
+
+    true
+end
+
+# Helper to find the index of the reaction with a given reactant and product complex.
+function edgeindex(imat, src::T, dst::T) where T <: Int
+    for i in 1:size(imat, 2)
+        (imat[src, i] == -1) && (imat[dst, i] == 1) && return i
+    end
+    error("This edge does not exist in this reaction graph.")
+end
+
+function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{<:Pair})
+    pdict = Dict(parametermap)
+    isdetailedbalanced(rs, pdict)
+end
+
+function isdetailedbalanced(rs::ReactionSystem, parametermap::Tuple{<:Pair})
+    pdict = Dict(parametermap)
+    isdetailedbalanced(rs, pdict)
+end
+
+function isdetailedbalanced(rs::ReactionSystem, parametermap)
+    error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
+end
+
 """
     iscomplexbalanced(rs::ReactionSystem, parametermap)
 
@@ -731,7 +824,6 @@ Constructively compute whether a network will have complex-balanced equilibrium
 solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). 
 Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
 """
-
 function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict)
     if length(parametermap) != numparams(rs)
         error("Incorrect number of parameters specified.")
@@ -808,7 +900,6 @@ end
     constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple 
     of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
 """
-
 function ratematrix(rs::ReactionSystem, rates::Vector{Float64})
     complexes, D = reactioncomplexes(rs)
     n = length(complexes)
@@ -903,7 +994,6 @@ end
     Returns the matrix of a basis of cycles (or flux vectors), or a basis for reaction fluxes for which the system is at steady state. 
     These correspond to right eigenvectors of the stoichiometric matrix. Equivalent to [`fluxmodebasis`](@ref). 
 """
-
 function cycles(rs::ReactionSystem)
     nps = get_networkproperties(rs)
     nsm = netstoichmat(rs)
@@ -938,7 +1028,6 @@ end
 
     See documentation for [`cycles`](@ref). 
 """
-
 function fluxvectors(rs::ReactionSystem)
     cycles(rs)
 end
@@ -950,7 +1039,6 @@ end
 
     Check if a reaction network obeys the conditions of the deficiency one theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class.
 """
-
 function satisfiesdeficiencyone(rn::ReactionSystem)
     all(r -> ismassaction(r, rn), reactions(rn)) ||
         error("The deficiency one theorem is only valid for reaction networks that are mass action.")
@@ -973,7 +1061,6 @@ end
 
     Check if a reaction network obeys the conditions of the deficiency zero theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class, this equilibrium is asymptotically stable, and this equilibrium is complex balanced.
 """
-
 function satisfiesdeficiencyzero(rn::ReactionSystem)
     all(r -> ismassaction(r, rn), reactions(rn)) ||
         error("The deficiency zero theorem is only valid for reaction networks that are mass action.")
@@ -988,7 +1075,6 @@ end
 
     Note: This function currently only works for networks of deficiency one, and is not currently guaranteed to return *all* the concentration-robust species in the network. Any species returned by the function will be robust, but this may not include all of them. Use with caution. Support for higher deficiency networks and necessary conditions for robustness will be coming in the future.  
 """
-
 function robustspecies(rn::ReactionSystem)
     complexes, D = reactioncomplexes(rn)
     nps = get_networkproperties(rn)
diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl
index decdbaa0ac..21f7cb402e 100644
--- a/test/network_analysis/network_properties.jl
+++ b/test/network_analysis/network_properties.jl
@@ -355,6 +355,7 @@ let
     k = rand(rng, numparams(rn))
     rates = Dict(zip(parameters(rn), k))
     @test Catalyst.iscomplexbalanced(rn, rates) == true 
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
 end
 
 ### STRONG LINKAGE CLASS TESTS
@@ -498,6 +499,107 @@ let
     @test isempty(cyclemat)
 end
 
+### Complex and detailed balance tests
+
+# The following network is conditionally complex balanced - it only 
+
+# Reversible, forest-like deficiency zero network - should be detailed balance for any choice of rate constants. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), A <--> D
+        (k5, k6), A + D <--> E
+        (k7, k8), A + D <--> G
+        (k9, k10), G <--> 2F
+        (k11, k12), A + E <--> H
+    end
+
+    k1 = rand(rng, numparams(rn))
+    rates1 = Dict(zip(parameters(rn), k1))
+    k2 = rand(StableRNG(232), numparams(rn))
+    rates2 = Dict(zip(parameters(rn), k2))
+
+    rcs, D = reactioncomplexes(rn)
+    @test Catalyst.isforestlike(rn) == true
+    @test Catalyst.isdetailedbalanced(rn, rates1) == true
+    @test Catalyst.isdetailedbalanced(rn, rates2) == true
+end
+
+# Simple connected reversible network
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B
+        (k3, k4), B <--> C
+        (k5, k6), C <--> A
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    rates1 = [:k1=>1.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0]
+    @test Catalyst.isdetailedbalanced(rn, rates1) == true
+    rates2 = [:k1=>2.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0]
+    @test Catalyst.isdetailedbalanced(rn, rates2) == false 
+end
+
+# Independent cycle tests: the following reaction entwork has 3 out-of-forest reactions. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), A <--> D
+        (k5, k6), B + C <--> D
+        (k7, k8), A + D <--> E
+        (k9, k10), G <--> 2F
+        (k11, k12), A + D <--> G
+        (k13, k14), G <--> E
+        (k15, k16), 2F <--> E
+        (k17, k18), A + E <--> H
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    k = rand(rng, numparams(rn))
+    p = parameters(rn)
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+
+    # Adjust rate constants to obey the independent cycle conditions. 
+    rates[p[6]] = rates[p[1]]*rates[p[4]]*rates[p[5]] / (rates[p[2]]*rates[p[3]]) 
+    rates[p[14]] = rates[p[13]]*rates[p[11]]*rates[p[8]] / (rates[p[12]]*rates[p[7]])
+    rates[p[16]] = rates[p[8]]*rates[p[15]]*rates[p[9]]*rates[p[11]] / (rates[p[7]]*rates[p[12]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == true
+end
+
+# Deficiency two network: the following reaction network must satisfy both the independent cycle conditions and the spanning forest conditions 
+let
+    rn = @reaction_network begin
+        (k1, k2), 3A <--> A + 2B
+        (k3, k4), A + 2B <--> 3B
+        (k5, k6), 3B <--> 2A + B
+        (k7, k8), 2A + B <--> 3A
+        (k9, k10), 3A <--> 3B
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    @test Catalyst.edgeindex(D, 1, 2) == 1
+    @test Catalyst.edgeindex(D, 4, 3) == 6
+    k = rand(rng, numparams(rn))
+    p = parameters(rn)
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+
+    # Adjust rate constants to fulfill independent cycle conditions. 
+    rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]])
+    rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+    # Should still fail - doesn't satisfy spanning forest conditions. 
+
+    # Adjust rate constants to fulfill spanning forest conditions. 
+    cons = rates[p[6]] / rates[p[5]]
+    rates[p[1]] = rates[p[2]] * cons
+    rates[p[9]] = rates[p[10]] * cons^(3/2)
+    rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]])
+    rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == true
+end
+
 ### Other Network Properties Tests ###
 
 # Tests outgoing complexes matrices (1). 
@@ -637,6 +739,108 @@ let
     @test Catalyst.robustspecies(EnvZ_OmpR) == [6]
 end
 
+
+### Complex and detailed balance tests
+
+# The following network is conditionally complex balanced - it only 
+
+# Reversible, forest-like deficiency zero network - should be detailed balance for any choice of rate constants. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), A <--> D
+        (k5, k6), A + D <--> E
+        (k7, k8), A + D <--> G
+        (k9, k10), G <--> 2F
+        (k11, k12), A + E <--> H
+    end
+
+    k1 = rand(rng, numparams(rn))
+    rates1 = Dict(zip(parameters(rn), k1))
+    k2 = rand(StableRNG(232), numparams(rn))
+    rates2 = Dict(zip(parameters(rn), k2))
+
+    rcs, D = reactioncomplexes(rn)
+    @test Catalyst.isforestlike(rn) == true
+    @test Catalyst.isdetailedbalanced(rn, rates1) == true
+    @test Catalyst.isdetailedbalanced(rn, rates2) == true
+end
+
+# Simple connected reversible network
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B
+        (k3, k4), B <--> C
+        (k5, k6), C <--> A
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    rates1 = [:k1=>1.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0]
+    @test Catalyst.isdetailedbalanced(rn, rates1) == true
+    rates2 = [:k1=>2.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0]
+    @test Catalyst.isdetailedbalanced(rn, rates2) == false 
+end
+
+# Independent cycle tests: the following reaction entwork has 3 out-of-forest reactions. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), A <--> D
+        (k5, k6), B + C <--> D
+        (k7, k8), A + D <--> E
+        (k9, k10), G <--> 2F
+        (k11, k12), A + D <--> G
+        (k13, k14), G <--> E
+        (k15, k16), 2F <--> E
+        (k17, k18), A + E <--> H
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    k = rand(rng, numparams(rn))
+    p = parameters(rn)
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+
+    # Adjust rate constants to obey the independent cycle conditions. 
+    rates[p[6]] = rates[p[1]]*rates[p[4]]*rates[p[5]] / (rates[p[2]]*rates[p[3]]) 
+    rates[p[14]] = rates[p[13]]*rates[p[11]]*rates[p[8]] / (rates[p[12]]*rates[p[7]])
+    rates[p[16]] = rates[p[8]]*rates[p[15]]*rates[p[9]]*rates[p[11]] / (rates[p[7]]*rates[p[12]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == true
+end
+
+# Deficiency two network: the following reaction network must satisfy both the independent cycle conditions and the spanning forest conditions 
+let
+    rn = @reaction_network begin
+        (k1, k2), 3A <--> A + 2B
+        (k3, k4), A + 2B <--> 3B
+        (k5, k6), 3B <--> 2A + B
+        (k7, k8), 2A + B <--> 3A
+        (k9, k10), 3A <--> 3B
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    @test Catalyst.edgeindex(D, 1, 2) == 1
+    @test Catalyst.edgeindex(D, 4, 3) == 6
+    k = rand(rng, numparams(rn))
+    p = parameters(rn)
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+
+    # Adjust rate constants to fulfill independent cycle conditions. 
+    rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]])
+    rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+    # Should still fail - doesn't satisfy spanning forest conditions. 
+
+    # Adjust rate constants to fulfill spanning forest conditions. 
+    cons = rates[p[6]] / rates[p[5]]
+    rates[p[1]] = rates[p[2]] * cons
+    rates[p[9]] = rates[p[10]] * cons^(3/2)
+    rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]])
+    rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == true
+end
+
 ### DEFICIENCY ONE TESTS
 
 # Fails because there are two terminal linkage classes in the linkage class