Cisneros [@Dombrowski; @Cisneros] inspired the numerical work aimed at understanding the influence of hydrodynamics in the formation of coherent structures of micro-organisms, simulated as spherical squirmers [@Ishikawa2006] in 3D [@Ishikawa2007a; @Ishikawa2008JFM] and quasi-2D [@Ishikawa2008PRL]. The features of a 3D squirmer suspension have been numerically characterized using the smoothed profile method [@molinaSoftMatt13; @molina_PRE16; @JPS_yamamoto] and Lattice Boltzmann (LB), either in bulk [@alarcon2013] or in response to a steady Couette flow [@LlopisSoftMatter13]. While Evans et al. [@evansPoF] used Stokesian Dynamics to detect the instability of the isotropic state in a 3D suspension of weak pullers, Delmotte and coworkers [@Delmotte] applied a force-coupling method to investigate, by large-scale simulations, the emergence of a polar order in a 3D squirmer suspension.

More recently, experiments by Ref. [@SoftMatter_35; @SoftMatter_33; @PRL_bacteriaFastest; @Cecille; @Thutupalli] have inspired numerical simulations of 2D and quasi-2D suspensions of squirmers interacting via a repulsive potential. When a suspension (simulated by means of multi-particle collision dynamics) was confined between parallel walls (quasi-2D geometry) [@Zottl_prl], active particles were forming clusters whose size depended on the system’s concentration with the largest being made of pullers. The suspension prepared in such confined geometry would undergo phase separation into a dilute and a dense phase [@Zottl_softmatter16], differently 