How does one solve a bivariate normal density function? If the exponent of $e$ of a bivariate normal density is 
$$\frac{-1}{54} *(x^2+4y^2+2xy+2x+8y+4) \\\text{find } \sigma_{1},\sigma_{2} \text{ and } p \text{ given that } \mu_{1} =0 \text{ and } \mu_{2}=-1.  $$

One must use this definition to solve.
A pair of random variables $X$ and $Y$ have a bivariate normal distribution and they are referred to as jointly normally distribed random variables if and if their joint probability density is given by 

If $X$ and $Y$ have a bivariate normal distribution normal distribution the conditional density of $Y$ given $X =x$ is a normal distribution with the mean 
$\mu_{Y|x} = \mu_{2} + \rho \cdot \frac{\sigma_2}{\sigma_1}\cdot (x-\mu_{1})$
and the variance
$\sigma^{2}_{Y|x} = \sigma^{2}_{2}(1-\rho^2)$
and the conditonal density of $X$ given $Y=y$ is a normal distribution with the mean
$\mu_{X|y} = \mu_{1} + \rho\cdot \frac{\sigma_{1}}{\sigma_{2}}\cdot (y-\mu_{2})$
and the variance 
$\sigma^{2}_{X|y} = \sigma^{2}_{1}(1-\rho^2)$
Does one re arrange the formula so that one can just plug in what $\mu_{1} \text{ and } \mu_{2}$ are? Does one substitute like this?
$$f(x,y) = \frac{-1}{54} *(x^2+4y^2+2xy+2x+8y+4) \\ = \frac{1}{-54} \left[ (x-\mu_{1})^2 +2(x-\mu_{1})(y-\mu_{2})+4(y-\mu_{2})^2 \right]$$
EDIT : Thanks to the valiant work of Anatoly I understand one needs to compare coefficients in order to derive $\rho,\sigma_{1},\sigma_{2}$ 
$$\color{lime}{(1)}: \frac{1}{2(1-p^2)\sigma^2_{1}} = \frac{1}{54}$$
However after this one I do not know what to do does anyone know how to do the rest?

If we substitute $\mu_1=0$ and $\mu_2=-1$ in the exponent of $e $ of the joint probability density formula,   we get
$$\displaystyle -\frac 1{2(1-\rho^2)}\left[\left(\frac{x}{\sigma_1}\right)^2
-2\rho   \left(\frac{x}{\sigma_1}\right)   \left(\frac{y+1}{\sigma_2}\right)     +\left(\frac{y+1}{\sigma_2}\right)^2\right]$$
Expanding it, we have
$$-\frac 1{2(1-\rho^2)}   \left[\frac{1}{\sigma_1^2}x^2   +\frac{1}{\sigma_2^2} y^2
-  \frac{ 2\rho }{\sigma_1 \sigma_2}  xy -   \frac{ 2\rho }{\sigma_1 \sigma_2}  x      +  \frac{2y}{\sigma_2^2} +  \frac{1}{\sigma_2^2}  \right]$$
which can be written, if we move a given factor $k^2$ into the brackets, as
$$-\frac {1}{2k^2(1-\rho^2)}   \left[\frac{k^2}{\sigma_1^2}x^2   +\frac{k^2}{\sigma_2^2} y^2
-  \frac{ 2k^2 \rho }{\sigma_1 \sigma_2}  xy -   \frac{ 2 k^2\rho }{\sigma_1 \sigma_2}  x     +  \frac{2k^2y}{\sigma_2^2} +  \frac{k^2}{\sigma_2^2}  \right]$$ 
Comparing this expression with 
$$-\frac {1}{54}   \left[x^2   +4 y^2+2  xy +2 x  +  8y + 4 \right]$$
which is provided in the OP we directly get $\sigma_1=k \,\,$, $\sigma_2=k/2\,\,$ (note that both variances are by definition positive), and $\rho=-1/2\,\,$. So we have 
$$-\frac {1}{2k^2 \left(1-\left (\frac {1}{2} \right)^2 \right)} =-\frac {1}{54} $$
and then $k=\pm 6\,\,$. This leads to  $\sigma_1^2=36\,\,$   and $\sigma_2^2=9\,\,$. 
