Metasurfaces represent one of the most vibrant fields of modern science and technology. A metasurface is a complex electromagnetic structure, that is typically deeply subwavelength in thickness, electrically large in transverse size and composed of subwavelength scattering particles with extremely small features; it may generally be bianisotropic, spacevarying and time-varying, nonlinear, curved and multiphysics. With such complexity, the design of a metasurface requires a holistic approach, involving synergistic synthesis and analysis operations, based on a solid model. The Generalized Sheet Transition Conditions (GSTCs), combined with bianisotropic surface susceptibility functions, provide such a model, and allow now for the design of sophisticated metasurfaces, which still represented a major challenge a couple of years ago. This paper presents this problematic, focusing on the computational analysis of metasurfaces via the GSTC-susceptibility approach. It shows that this analysis plays a crucial role in the holistic design of metasurfaces, and overviews recently reported related frequency-domain (FDFD, SD-IE, FEM) and time-domain (FDTD) computational techniques.