For each optimization task with associated constraints (i.e., RC2, RC3, and RC4), the solver is initialized independently with random states default to the solver. More specifically, RC3 reactions are not found by applying constraints to RC2 results; instead, RC3 reactions are found by applying the RC3 constraints to the same reference species set used for finding RC3 reactions. Similarly, RC4 reactions are not found by further fine-tuning the results of RC3 reactions. Therefore, it is possible for the solver to find more RC3 reactions than RC2 reactions. Nonetheless, RC3 by definition is a subset of RC2, and we refer to RC3 reactions as reactions that are at least as strict as RC3. To solve the MILP, lp_solve is used, although a wrapper to use Pyomo 80,81 is also implemented. In general, we find lp_solve to be significantly faster most of the time, though it is also prone to occasionally not converging to a solution. For each species, our implementation attempts to find at most 5 isodemic reactions and take the median value of ∆ iso f H • t from found reactions. Furthermore, the values of ∆ qm f H • for the reference species used in fitting IDRs include AECs so that the resulting ∆ iso f H • t values are at least somewhat reasonable and close to the reference values.