The complete expression for the integral is
$$
D(E) = \frac{1}{4\pi^2}\int_{E(k_x,k_y)=E}
\frac{d\mathbf{l}}{\left|\nabla_k E(k_x,k_y)\right|}
$$
This is the integration over a curve which in your case is described by the following equation:
$$
E = E_0 + a_x k_x^2 - a_y k_y^2
$$
Depending on what energy value you substitute, you should select either $k_x$ or $k_y$ as the integration variable. The variable selected should change continuously from $-\infty$ to $+\infty$ on the curve.
This is just for convenience. One variable is convenient on the one side of the singularity, another is on the other side.
If you select say $k_x$ the integral will look as follows:
$$
D(E) = 4\cdot\frac{1}{4\pi^2}\int_0^\infty
\frac{1}{\left|\nabla_k E\bigl(k_x,k_y(k_x)\bigr)\right|}
\left|\frac{d\mathbf{l}(k_x)}{dk_x}\right| dk_x
$$
This is a usual one-dimensional integral.
The factor $4$ is here because there are four identical parts of the curve while we are integrating over only one.
In calculating the electron dispersion you probably obtained the diagonalized Hamiltonian in the momentum space
$$
H=\sum_\mathbf{k}\left[c^{\dagger}_{\mathbf{k}A},c^{\dagger}_{\mathbf{k}B}\right]\left[\begin{array}{cc}0 & \Delta(\mathbf{k})\\ \Delta^{\dagger}(\mathbf{k}) &0\end{array}\right]\left[\begin{array}{c}c_{\mathbf{k}A} \\ c_{\mathbf{k}B}\end{array}\right].
$$
If you you chose your $x$ axis along the zigzag direction (arXiv:1004.3396), the two nonequivalent Dirac valleys are $\mathbf{K}_\kappa=\left(\kappa\frac{4\pi}{3\sqrt{3}a},0\right)$, $\kappa=\pm1$ and $\mathbf{K}_{-1}=\mathbf{K}^{\prime}$, where $a$ is the C-C distance. Then $\Delta(\mathbf{k})=-t\left(1+e^{-i\mathbf{k}\cdot\mathbf{a}_1}+e^{-i\mathbf{k}\cdot\mathbf{a}_2}\right)$, where $t$ is the hopping term, and $\mathbf{a}_1=\left(\sqrt{3}a/2,3a/2\right)$ and $\mathbf{a}_2=\left(-\sqrt{3}a/2,3a/2\right)$ are the lattice vectors.
Taylor expanding $\Delta(\mathbf{k})$ up to linear terms around those two points you obtain
$$
\Delta(\mathbf{k})=\kappa\frac{3ta}{2}q_x-i\frac{3ta}{2}q_y
$$
where $\mathbf{q}$ is the displacement momenta from the $\mathbf{K}_\kappa$ point. Promoting these displacement momenta to operators you obtain the Hamiltonian
$$
H=\hbar v_F\left[\begin{array}{cccc}0 & q_x-iq_y & 0 & 0\\q_x+iq_y & 0 & 0 & 0\\0 & 0 & 0 & -q_x-iq_y\\0 & 0 & -q_x+iq_y & 0\end{array}\right]
$$
where $v_F=\frac{3ta}{2\hbar}$ is the Fermi velocity. This is in $\left[\Psi_{A\mathbf{K}},\Psi_{B\mathbf{K}},\Psi_{A\mathbf{K}^{\prime}},\Psi_{B\mathbf{K}^{\prime}}\right]^T$ basis, if you rearrange your basis as $\left[\Psi_{A\mathbf{K}},\Psi_{B\mathbf{K}},\Psi_{B\mathbf{K}^{\prime}},\Psi_{A\mathbf{K}^{\prime}}\right]^T$ you get the compact form
$$
H=\hbar v_F\tau_z\otimes\boldsymbol{\sigma}\cdot\mathbf{k}
$$
where $\tau_z$ acts in the valley space. This is similar to the Dirac-Weyl equation for relativistic massless particles, where instead of $v_F$ you get the speed of light
$$
H=\pm\hbar c\boldsymbol{\sigma}\cdot\mathbf{k}
$$
where $+$ denotes right-handed antineutrions, and $-$ denotes left-handed neutrions. The differences are that $\boldsymbol{\sigma}=\left(\sigma_x,\sigma_y\right)$ for graphene acts in pseudospin space and $\boldsymbol{\sigma}=\left(\sigma_x,\sigma_y,\sigma_z\right)$ for neutrinos acts in real spin space.
Best Answer
So I do not know if this is still relevant, but maybe someone else can profit from my answer?
To analyse the behaviour around the critical points $\mathbf k^*$, you can just Taylor expand $\varepsilon_\mathbf{k}$ around these points and calculate the DOS with your or this $$ \rho(\varepsilon) = \frac{1}{(2\pi)^d} \int_{BZ}d^dk\,\delta(\varepsilon - \varepsilon_\mathbf k) $$ formula (valid in $d$ dimensions). For a square 2D lattice, this will result in a logarithmic van-Hove singularity at $\varepsilon = 0$ (ie. $\mathbf k^* = (0,\pi),(\pi,0) $) and no divergence at $\varepsilon = \pm 4t$ (ie. $\mathbf k^* = (0,0),(\pi,\pi)$).
Doing so, you will see that the integrable zeroes of the gradient of the dispersion relation do not result in singularities in the DOS. For example, the integral of $ 1/x^2 $ over $[-\infty,\infty]$ is finite, even though the integrand diverges at 0.
However, the value of $\varepsilon$ does not matter, so the fact that $\varepsilon = 0$ does not tell us anything about a potential singularity. Also the last mentioned points do not result in a zero gradient! (See comments!)