In this work heme models with four [Fe(II)(P)], five [Fe(II)(P)Im], [Fe(II)(P)(Im)O2] and six ligands [Fe(II)(P)(Im)O2], where P = porphyrin, with different spin states (ms =5, 3 and 1) of the iron atom were investigated using relativistic‐corrected quantum chemistry methods (PW6B95‐D3‐DKH/jorge‐TZP‐DKH). Dependence of the iron‐ligand bond properties on (i) spin state and (ii) number of ligands were analyzed using natural bond orbital analysis, electron density topology, electrostatic potential and electron localization function. It is shown that reversible binding of O2 is possible in case of formation of semicoordination bond between Fe(II) and imidazole. Binding of the fifth and sixth ligand from the energetic and orbital points of view is more favorable for the triplet Fe(II) state. At the same time for the six‐coordinated complex [Fe(II)(P)(Im)O2] interconversion of Fe(II) electrons of valent 3d orbital from quintet to triplet and vice versa is possible under thermal fluctuations (energy barriers less than 2 kcal/mol).