In multiaxial low cycle fatigue, the directional distribution of initiated cracks affects the subsequent growth of a dominant crack. So, it is important to clarify the characteristics of crack initiation to improve the life evaluation in multiaxial fatigue. In this work, a simulation pr ocedure of fatigue crack initiation under multiaxial stresses was first established based on the modeling of crack initiation in slip-band systems assumed in polycrystalline metals. Especially, the directional distribution of initiated cracks was focused on in the analysis. By introduci ng a new parameter “intensity of slip-band crack”, simulated directional distributions of initiated cracks were compared with experimental ones observed in biaxial fatigue tests using tubular specimens of pure copper. The qualitative tendency in the material with homogeneous microstructu re is well simulated by using the proposed procedure. Finally, the procedure was modified by modeling inclusions, and simulations were carried out to estimate directional distribution of cracks in a material with inclusions. The applicability of the modified modeling was also discussed.