A new model of sulfur solubility in mafic and/or ultramafic silicate magmas, which accounts for the effects of pressure, temperature, oxygen fugacity, major element, and Ni contents in the silicate melt and the coexisting sulfide liquid, is presented in this paper. The model postulates the existence of positively charged Fe-Ni sulfide complexes in the melt of a general formula (FeyNi1−y)zS2(z−1)+, which are formed as a result of complexation reactions between the sulfide-forming ions (Fe2+, Ni2+, S2−) and (Fe,Ni)S species in the silicate liquid. The new model can explain both the anomalously high S solubility in iron-enriched silicate systems and the “parabola-like” dependence of S contents in silicate melts on their Fe content. The proposed mechanism of sulfide solubility was calibrated on a dataset of 213 anhydrous experimental glasses (both Ni free and Ni bearing) and 53 S-saturated MORB glasses, and incorporated into a new version of the COMAGMAT (v. 5) magma crystallization model. The COMAGMAT-5 model can estimate sulfur concentration at sulfide saturation (SCSS) in a wide range of experimental and natural compositions, including Fe/Ni variations in silicate melts and coexisting sulfides. Despite relatively low concentrations, nickel is shown to have a pronounced effect on S solubility, causing significant variations in the onset of sulfide immiscibility in melts with otherwise similar major element compositions. An application example of the new SCSS model to “B-1 magma” proposed as parent for the Lower and Lower Critical zones of the Rustenburg Layered Suite, Bushveld Complex, is discussed.