@@ -21,29 +21,31 @@ function isconcentrationrobust(rn::ReactionSystem; p::VarMapType = Dict())
2121 S, D = CatalystNetworkAnalysis. removespec (rn, i)
2222 (Catalyst. deficiency (rn) == CatalystNetworkAnalysis. deficiency (S, D)) &&
2323 begin
24- possibly_ACR = true
25- break
26- end
24+ possibly_ACR = true
25+ break
26+ end
2727 end
2828 (possibly_ACR == false ) && return :NO_ACR
2929 end
3030
31- # Convert parameter values to rational values.
31+ # Convert parameter values to rational values.
3232 try
3333 ftype, stype = eltype (p). parameters
3434 stype <: Rational || (p = Dict {ftype, Rational} (p))
3535 catch InexactError
3636 return :INEXACTPARAMS
3737 end
3838
39- # Compute steady state ideal.
39+ # Compute steady state ideal.
4040 sfr_f = eval (C. SFR (rn; p = p))
4141 R,
42- polyvars = polynomial_ring (QQ,
42+ polyvars = polynomial_ring (
43+ QQ,
4344 vcat (
4445 map (s -> Symbolics. tosymbol (s, escape = false ), species (rn)),
4546 map (p -> Symbolics. tosymbol (p), parameters (rn))
46- ))
47+ )
48+ )
4749 sfr = sfr_f (polyvars... )
4850 specs = polyvars[1 : numspecies (rn)]
4951
@@ -64,7 +66,7 @@ function isconcentrationrobust(rn::ReactionSystem; p::VarMapType = Dict())
6466 return isempty (p) ? :GLOBAL_ACR : :MASS_ACTION_ACR
6567 end
6668
67- # Iterate over elimination ideals I ∩ Q[x_i] and check for unique positive root.
69+ # Iterate over elimination ideals I ∩ Q[x_i] and check for unique positive root.
6870 # Create a dummy univariate ring
6971 r_ξ, ξ = polynomial_ring (QQ, :ξ )
7072
@@ -74,7 +76,7 @@ function isconcentrationrobust(rn::ReactionSystem; p::VarMapType = Dict())
7476 iszero (g) && continue
7577
7678 # Generate the univariate polynomial corresponding to the generator
77- coeff_pos = [exponent_vector (g, j)[i]+ 1 for j in 1 : length (Oscar. terms (g))]
79+ coeff_pos = [exponent_vector (g, j)[i] + 1 for j in 1 : length (Oscar. terms (g))]
7880 coeffs = zeros (QQFieldElem, first (coeff_pos))
7981 coeffs[coeff_pos] .= Oscar. coefficients (g)
8082 poly = r_ξ (coeffs)
@@ -109,10 +111,10 @@ function linearelements(I::Ideal, numspecies::Int)
109111 return linearidxs
110112end
111113
112- # TODO : Other methods.
114+ # TODO : Other methods.
113115# Compute decompositions of the positive-restriction ideal
114116# Numerically checking for ACR
115- # Compute the steady-state parameterization, check for constant components.
117+ # Compute the steady-state parameterization, check for constant components.
116118# Algorithm 6.1 for finding all pairs of ACR candidates
117119
118120"""
@@ -128,40 +130,40 @@ function robustspecies_δ1(rn::ReactionSystem)
128130 error (" This algorithm currently only checks for robust species in networks with deficiency one." )
129131 end
130132
131- # A species is concentration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes
132- # belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some
133- # multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B).
133+ # A species is concentration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes
134+ # belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some
135+ # multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B).
134136
135137 if ! nps. checkedrobust
136138 tslcs = terminallinkageclasses (rn)
137139 Z = complexstoichmat (rn)
138140
139- # Find the complexes that do not belong to a terminal linkage class
141+ # Find the complexes that do not belong to a terminal linkage class
140142 nonterminal_complexes = deleteat! ([1 : length (complexes);], vcat (tslcs... ))
141143 robust_species = Int64[]
142144
143145 for c_s in nonterminal_complexes, c_p in nonterminal_complexes
144146
145147 (c_s >= c_p) && continue
146- # Check the difference of all the combinations of complexes. The support is the set of indices that are non-zero
148+ # Check the difference of all the combinations of complexes. The support is the set of indices that are non-zero
147149 suppcnt = 0
148150 supp = 0
149151 for i in 1 : size (Z, 1 )
150152 (Z[i, c_s] != Z[i, c_p]) && (suppcnt += 1 ; supp = i)
151153 (suppcnt > 1 ) && break
152154 end
153155
154- # If the support has length one, then they differ in only one species, and that species is concentration robust.
156+ # If the support has length one, then they differ in only one species, and that species is concentration robust.
155157 (suppcnt == 1 ) && (supp ∉ robust_species) && push! (robust_species, supp)
156158 end
157159 nps. checkedrobust = true
158160 nps. robustspecies = robust_species
159161 end
160162
161- nps. robustspecies
163+ return nps. robustspecies
162164end
163165
164166# #### DEFICIENCY METHODS
165- # These are based on the computation of "robust ratios" between complexes - complexes y, y' such that x^y / x^y' is constant for all positive steady states.
167+ # These are based on the computation of "robust ratios" between complexes - complexes y, y' such that x^y / x^y' is constant for all positive steady states.
166168# For def 0 networks, this is true for all networks in the same linkage class
167- # For def 1 networks, this is true for all
169+ # For def 1 networks, this is true for all
0 commit comments