The $SO(4)$ Plebanski action yields a first order formulation of Euclidean General Relativity as a constrained (topological) BF-theory. It depends on a $so(4)$ connection 1-form $\omega^{IJ} = \omega_{\mu}^{IJ}dx^{\mu}$ through the curvature $F^{IJ}(\omega) = d_{\omega}\omega^{IJ}$, with $d_{\omega}(\cdot)$ the exterior covariant differential, on a $so(4)$ 2-form $B^{IJ} = B^{IJ}dx^{\mu}\wedge dx^{\nu}$ and a scalar symmetric traceless matrix $\phi_{[IJ][KL]}$ ($\epsilon^{IJKL}\phi_{IJKL} = 0$) with components acting as Lagrange multipliers. Assuming zero cosmological constant ($\Lambda$ = 0), the action reads: \begin{equation} S[\omega; B, \phi] = \int_{M} B^{IJ}\wedge F_{IJ}(\omega) - \frac{1}{2}\phi_{IJKL}B^{IJ}\wedge B^{KL} \end{equation}
Varying this action with respect to the multipliers $\phi_{IJKL}$ is supposed to provide the well-known **simplicity constraints** \begin{equation} B^{IJ}\wedge B^{KL} = e\epsilon^{IJKL}, ~~~ \text{with}~~ e = \frac{1}{4!}\epsilon_{IJKL}B^{IJ}\wedge B^{KL} \end{equation} imposing the on-shell antisymmetric form of the field $B^{IJ}$. 
However, I don't seem to understand how the variation is to be performed in order to obtain those constraints. I understand it is necessary to take into account the symmetry properties of the fields $\phi_{IJKL}$, but none of the papers I have read regarding this topic detail that calculation. 

Does anyone know how this variation should be performed or any papers with further details?