We develop a stochastic model for the velocity gradients dynamics along a Lagrangian trajectory in isotropic and homogeneous turbulent flows. Comparing with different attempts proposed in the literature, the present model, at the cost of introducing a free parameter known in turbulence phenomenology as the intermittency coefficient, gives a realistic picture of velocity gradient statistics at any Reynolds number. To achieve this level of accuracy, we use as a first modelling step a regularized self-stretching term in the framework of the Recent Fluid Deformation (RFD) approximation that was shown to give a realistic picture of small scales statistics of turbulence only up to moderate Reynolds numbers. As a second step, we constrain the dynamics, in the spirit of Girimaji & Pope (1990), in order to impose a peculiar statistical structure to the dissipation seen by the Lagrangian particle. This probabilistic closure uses as a building block a random field that fulfils the statistical description of the intermittency, i.e. multifractal, phenomenon. To do so, we define and generalize to a statistically stationary framework a proposition made by Schmitt (2003). These considerations lead us to propose a non-linear and non-Markovian closed dynamics for the elements of the velocity gradient tensor. We numerically integrate this dynamics and observe that a stationary regime is indeed reached, in which (i) the gradients variance is proportional to the Reynolds number, (ii) gradients are typically correlated over the (small) Kolmogorov time scale and gradients norms over the (large) integral time scale (iii) the joint probability distribution function of the two non vanishing invariants Q and R reproduces the characteristic teardrop shape, (iv) vorticity gets preferentially aligned with the intermediate eigendirection of the deformation tensor and (v) gradients are strongly non-Gaussian and intermittent, a behaviour that we quantify by appropriate high order moments. Additionally, we examine the problem of rotation rate statistics of (axisymmetric) anisotropic particles as observed in Direct Numerical Simulations. Although our realistic picture of velocity gradient fluctuations leads to better results when compared to the former RFD approximation, it is still unable to provide an accurate description for the rotation rate variance of oblate spheroids.