Differential operator in Mathematica

I am fairly new to Mathematica and I wanted to know if I can define the following operator

$$E:=(x^2\sqrt{1-K(x^2+y^2+z^2}+y^2+z^2)\gamma^{1}\partial_{1}+(y^2\sqrt{1-K(x^2+y^2+z^2}+x^2+z^2)\gamma^{2}\partial_{2}$$

where the gamma's are actually matrices, but that isn't pertinent to the calculation; what is, however, is the position in which they are since they are three possible values that they can assume when multiplied together:

$$\mu=\nu=0 \rightarrow -2$$

$$\mu=\nu\in \{1,2,3\} \rightarrow 2$$

$$\mu\neq \nu \in\{0,1,2,3\} \rightarrow 0$$

since they satisfy the anticommutation relation, $$\{\gamma^{\mu},\gamma^{\nu}\}=\gamma^{\mu}\gamma^{\nu}+\gamma^{\nu}\gamma^{\mu}=2\eta^{\mu\nu}I_{4}$$, such that I can then square the operator and apply it to some function $$\psi(t,x,y,z)$$, $$E^2\psi$$. I've been trying to use Nest and have even been looking at other posts for some ideas, but I can't seem to wrap my finger on how I can get this to work. Any help????

https://mathematica.stackexchange.com/questions/72433/polynomial-expansion-of-operatorfunct