Among the available tools for the early diagnosis of breast cancer, the elastographic technique based on ultrasounds has many advantages such as the noninvasive measure, the absence of ionizing effects, the high tolerability by patients, and the wide diffusion of the ecographic machines. However this diagnostic procedure is strongly affected by many subjective factors and is considered not reliable enough even to reduce the number of biopsies used to identify the nature of lesions. Therefore in the literature experimental and numerical simulations on physical and virtual phantoms are presented to test and validate procedures and algorithms and to interpret elastosonographic results. In this work, first a description of the elastographic technique and a review of the principal finite element (FE) models are provided and second diagnostic indexes employed to assess the nature of a lump mass are presented. As advances in FE simulations of elastosonography, axisymmetric phantom, and anthropomorphic models are described, which, with respect to the literature, include some features of breast mechanics. In particular deterministic analyses were used to compare the various details of virtual elastograms and also to investigate diagnostic indexes with respect to the regions where strains were considered. In order to improve the reliability of the elastosonographic procedure, univariate and multivariate sensitivity analyses, based on a probabilistic FE approach, were also performed to identify the parameters that mostly influence the deformation contrast between healthy and cancerous tissues. Moreover, synthetic indicators of the strain field, such as the strain contrast coefficient, were evaluated in different regions of interest in order to identify the most suitable for lesion type assessment. The deterministic analyses show that the malignant lesion is characterized by a uniform strain inside the inclusion due to the firmly bonding condition, while in the benign inclusion (loosely bonded) a strain gradient is observed independently from the elastic modulus contrast. The multivariate analyses reveal that the strain contrast depends linearly on the relative stiffness between the lesion and the healthy tissue and not linearly on the interface friction coefficient. The anthropomorphic model shows other interesting features, such as the layer or curvature effects, which introduce difficulties in selecting a reference region for strain assessment. The results show that a simple axisymmetric model with linear elastic material properties can be suitable to simulate the elastosonographic procedure although the breast curvature and layer distinction play a significant role in the strain assessment.