You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Constructively compute whether a network will have complex-balanced equilibrium
1663
-
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 ofvariable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...].
error("Incorrect number of parameters specified.")
1669
-
end
1670
-
1671
-
pmap =symmap_to_varmap(rs, parametermap)
1672
-
pmap =Dict(ModelingToolkit.value(k) => v for (k,v) in pmap)
1673
-
1674
-
sm =speciesmap(rs)
1675
-
cm =reactioncomplexmap(rs)
1676
-
complexes, D =reactioncomplexes(rs)
1677
-
rxns =reactions(rs)
1678
-
nc =length(complexes); nr =numreactions(rs); nm =numspecies(rs)
1679
-
1680
-
if!all(r->ismassaction(r, rs), rxns)
1681
-
error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.")
1682
-
end
1683
-
1684
-
rates = [substitute(rate, pmap) for rate inreactionrates(rs)]
1685
-
1686
-
# Construct kinetic matrix, K
1687
-
K =zeros(nr, nc)
1688
-
for c in1:nc
1689
-
complex = complexes[c]
1690
-
for (r, dir) in cm[complex]
1691
-
rxn = rxns[r]
1692
-
if dir ==-1
1693
-
K[r, c] = rates[r]
1694
-
end
1695
-
end
1696
-
end
1697
-
1698
-
L =-D*K
1699
-
S =netstoichmat(rs)
1700
-
1701
-
# Compute ρ using the matrix-tree theorem
1702
-
g =incidencematgraph(rs)
1703
-
R =ratematrix(rs, rates)
1704
-
ρ =matrixtree(g, R)
1705
-
1706
-
# Determine if 1) ρ is positive and 2) D^T Ln ρ lies in the image of S^T
0 commit comments