Parameterized ordinary differential equation systems are crucial for modeling in biochemical reaction networks under the assumption of mass-action kinetics. Various questions concerning the signs of multivariate polynomials in positive orthant arise from studying the solutions' qualitative behavior with respect to parameter values. In this work, we utilize circuit polynomials to find symbolic certificates of nonnegativity in order to provide further insight into the number of positive steady states of the n-site phosphorylation cycle model. This is a joint work with Elisenda Feliu, Nidhi Kaihnsa and Timo de Wolff.