A time-changed discretization for the Dirac equation is proposed. More precisely, we consider a Dirac equation with discrete space and continuous time perturbed by a time-dependent diffusion term $σ^2Ht^{2H-1}$ that seamlessly describes a latticizing version of the time-changed Fokker-Planck equation carrying the Hurst parameter $0<H<1$. Our model problem formulated on the space-time lattice $\mathbb{R}_{h,α}^n\times [0,\infty)$ ($h>0$ and $0<α<\frac{1}{2}$) preserves the main features of the Dirac-Kähler type discretization over the space-time lattice $h\mathbb{Z}^n\times [0,\infty)$ in case of $α,H \rightarrow 0$, and encompasses a regularization of Wilson's approach [Physical review D, 10(8), 2445, 1974] for values of $H$ in the range $0<H\leq \frac{1}{2}$ (limit condition $α\rightarrow \frac{1}{2}$). The main focus here is the representation of the solutions by means of discrete convolution formulae involving a kernel function encoded by (unnormalized) Hartman-Watson distributions -- ubiquitous on stochastic processes of Bessel type -- and the solutions of a semi-discrete equation of Klein-Gordon type. Namely, on our main construction the ansatz function $\widehat{\varPsi}_H(y)$ appearing on the discrete convolution representation may be rewritten as a Mellin convolution type integral involving the solutions $\varPsi(x,t|p)$ of a semi-discrete equation of Klein-Gordon type and a Lévy one-sided distribution $L_H(u)$ in disguise. Interesting enough, by employing Mellin-Barnes integral representations it turns out that the underlying solutions of Klein-Gordon type may be represented through generalized Wright functions of type ${~}_1Ψ_1$, that converge uniformly in case that the quantity $α+\frac{1}{2}$ may be regarded as an lower estimate for the Hurst parameter in the superdiffusive case (that is, if $α+\frac{1}{2}\leq H<1$).