Alternating subduction polarity along suture zones has been documented in several orogenic systems. Yet the mechanisms leading to this geometric inversion and the subsequent interplay between the contra-dipping slabs have been little studied. To explore such mechanisms, 3D numerical modelling of the Wilson cycle was conducted from continental rifting, breakup and oceanic spreading to convergence and self-consistent subduction initiation. In the resulting models, near-ridge subduction initiating with the formation of contra-dipping slab segments is an intrinsically 3D process controlled by earlier convergence-induced ridge swelling. The width of the slab segments is delimited by transform faults inherited from the rifting and ocean floor spreading stages. The models show that the number of contra-dipping slab segments depends mainly on the size of the oceanic basin, the asymmetry of the ridge and variations in kinematic inversion from divergence to convergence. Convergence velocity has been identified as a second-order parameter. The geometry of the linking zone between contra-dipping slab segments varies between two end-members governed by the lateral coupling between the adjacent slab segments: (1) coupled slabs generate wide, arcuate linking zones holding two-sided subduction; and (2) decoupled slabs generate narrow transform fault zones against which one-sided, contra-dipping slabs abut.