  We consider a problem from biological network analysis of determining regions
in a parameter space over which there are multiple steady states for positive
real values of variables and parameters. We describe multiple approaches to
address the problem using tools from Symbolic Computation. We describe how
progress was made to achieve semi-algebraic descriptions of the
multistationarity regions of parameter space, and compare symbolic results to
numerical methods. The biological networks studied are models of the
mitogen-activated protein kinases (MAPK) network which has already consumed
considerable effort using special insights into its structure of corresponding
models. Our main example is a model with 11 equations in 11 variables and 19
parameters, 3 of which are of interest for symbolic treatment. The model also
imposes positivity conditions on all variables and parameters.
  We apply combinations of symbolic computation methods designed for mixed
equality/inequality systems, specifically virtual substitution, lazy real
triangularization and cylindrical algebraic decomposition, as well as a
simplification technique adapted from Gaussian elimination and graph theory. We
are able to determine multistationarity of our main example over a
2-dimensional parameter space. We also study a second MAPK model and a symbolic
grid sampling technique which can locate such regions in 3-dimensional
parameter space.
