The drive to better understand complex phenomena like cardiac functions has led to mathematical models that are able to describe the multiscale structure and dynamics of the heart. The numerical simulation of these complex cardiac models is very challenging since it requires high-resolution space and time discretizations, as well as numerical tools that can handle natural discontinuities (for instance, disconnected cells). The development of effective and scalable solvers for the solution of these models has grown increasingly in the last decade, and modern computational architectures have increased the possibility to run large-scale simulations within reasonable computational times. In any case, the multiscale systems arising from discretizations of such models has required the development of specific techniques which can balance accuracy in the solution while being computationally competitive (in terms of efficiency and scalability). Indeed, being able to describe and computationally reproduce the many interactions between the macroscopic and microscopic events leads to numerical choices that yield large scale nonlinear algebraic systems of equations, exceeding millions of degrees of freedom. In this talk, we focus on the numerical simulation of the cardiac electrical activity, described either from a macroscopic point of view by means of the homogenized Bidomain model or from a microscopic point of view by considering a cell-by-cell model coupling. We propose ad-hoc preconditioned solvers that are able to describe the physiological features of these models, then we analyze them from a theoretical point of view and validate our analysis through extensive numerical tests, showing scalability and quasi-optimality of these solvers. These results provide a basis for improving parallel solvers for cardiac electrophysiology models that combine parallel efficiency and solution accuracy.