The change of potency to nucleate cracks in high cycle fatigue (HCF) at a primary nonmetallic inclusion in a martensitic gear steel due to the existence of a neighboring inclusion is computationally investigated using two-and three-dimensional elastoplastic finite element (FE) analyses. Fatigue indicator parameters (FIPs) are computed in the proximity of the inclusion and used to compare crack nucleation potency of various scenarios. The nonlocal average value of the maximum plastic shear strain amplitude is used in computing the FIP. Idealized spherical (cylindrical in 2D) inclusions with homogeneous linear elastic isotropic material properties are considered to be partially debonded, the worst case scenario for HCF crack nucleation as experimentally observed for similar systems (Furuya et al., 2004, “Inclusion-Controlled Fatigue Properties of 1800 Mpa-Class Spring Steels,” Metall. Mater. Trans. A, 35A(12), pp. 3737–3744; Harkegard, 1974, “Experimental Study of the Influence of Inclusions on the Fatigue Properties of Steel,” Eng. Fract. Mech., 6(4), pp. 795–805; Lankford and Kusenberger, 1973, “Initiation of Fatigue Cracks in 4340 Steel,” Metall. Mater. Trans. A, 4(2), pp. 553–559; Laz and Hillberry, 1998, “Fatigue Life Prediction From Inclusion Initiated Cracks,” Int. J. Fatigue, 20(4), pp. 263–270). Inclusion-matrix interfaces are simulated using a frictionless contact penalty algorithm. The fully martensitic steel matrix is modeled as elastic-plastic with pure nonlinear kinematic hardening expressed in a hardening minus dynamic recovery format. FE simulations suggest significant intensification of plastic shear deformation and hence higher FIPs when the inclusion pair is aligned perpendicular to the uniaxial stress direction. Relative to the reference case with no neighboring inclusion, FIPs decrease considerably when the inclusion pair aligns with the applied loading direction. These findings shed light on the anisotropic HCF response of alloys with primary inclusions arranged in clusters by virtue of the fracture of a larger inclusion during deformation processing. Materials design methodologies may also benefit from such cost-efficient parametric studies that explore the relative influence of microstructure attributes on the HCF properties and suggest strategies for improving HCF resistance of alloys.