Background: Estimation of tissue stiffness for detecting tissue abnormalities has been a major research interest in the last decades. Some existing methods start with defining a region-of-interest (ROI) and formulating the elasticity reconstruction of ROI as an inverse problem (IP). These commonly utilize the Finite Element Method (FEM) and compute Youngâs modulus distribution in ROI from displacement field measurements in response to tissue excitations. For solving the FEM IP, boundary conditions at the outer boundaries of ROI are necessitated [1]. Conventionally, these are set arbitrarily, as zero-displacement or -force constraints at chosen ROI locations, which may lead to errors in IP solutions [2]. Aims: Our goal is to reconstruct accurate elasticity parameters with harmonic excitations by alternatively approximating boundary constraints from measured displacements. Methods: Our proposed method starts with a displacement field tracking step obtained from harmonic excitation at angular frequency Ï, using e.g. cross-correlation [3] or optical flow [4]. We then start our iterative process with an initial assumption for Youngâs modulus and estimate compliance boundary conditions (CBC) as proposed in [5]. Afterwards, we use estimated CBC to reconstruct Youngâs modulus using FEM IP [2]. We iterate this process until updated Youngâs modulus converges or maximum number of iterations is reached as shown in (a). Note that our method has two nested iterations: elasticity IP iteration for reconstructing Youngâs modulus with FEM IP, which checks convergence for displacements, and CBC update iteration to update CBC using estimated Youngâs modulus, which converges for change of Young's modulus. Results: Simulations were performed using a 2D model with a circular inclusion embedded in a numerical phantom as in (b). Our simulation results over exterior iterations are demonstrated for Ï = 2Ïf at f = 70 Hz in (c). Sensitivity of IP parameter reconstruction for two empirically-set baseline comparisons, bottom and all fix are compared to our method CBC for different f in (d). RMSE error of these are shown in (e). Conclusions: The circular inclusion is recognizable for all angular frequencies and ROI choices using CBC, unlike empirically-set boundary conditions bottom and all fix. Our results show that using CBC as boundary conditions yields Young's modulus reconstruction outputs with less RMSE errors.