A simulation method to track the time evolution of coherent vorticity and current density, called coherent vorticity and current density simulation (CVCS), is developed for three-dimensional (3D) incompressible magnetohydrodynamic (MHD) turbulence. The vorticity and current density fields are, respectively, decomposed at each time step into two orthogonal components, the coherent and incoherent fields, using an orthogonal wavelet representation. Each of the coherent fields is reconstructed from the wavelet coefficients whose modulus is larger than a threshold, while their incoherent counterparts are obtained from the remaining coefficients. The two threshold values depend on the instantaneous kinetic and magnetic enstrophies. The induced coherent velocity and magnetic fields are computed from the coherent vorticity and current density using the Biot-Savart kernel. In order to compute the flow evolution, one should retain not only the coherent wavelet coefficients but also their neighbors in wavelet space, and the set of those additional coefficients is called the safety zone. CVCS is performed for 3D forced incompressible homogeneous MHD turbulence without mean magnetic field for a magnetic Prandtl number equal to unity and with 256^3 grid points. The quality of CVCS is assessed by comparing the results with a direct numerical simulation. It is found that CVCS with the safety zone well preserves the statistical predictability of the turbulent flow with a reduced number of degrees of freedom. CVCS is also compared with a Fourier truncated simulation using a spectral cutoff filter where the number of retained Fourier modes is similar to the number of the wavelet coefficients retained by CVCS. It is shown that the wavelet representation is more suitable than the Fourier representation, especially concerning the probability density functions of vorticity and current density.