Identification of tumor nodule in soft tissue: An inverse finite‐element framework based on mechanical characterization

Identification and characterization of nodules in soft tissue, including their size, shape, and location, provide a basis for tumor identification. This study proposes an inverse finite‐element (FE) based computational framework, for characterizing the size of examined tissue sample and detecting the presence of embedded tumor nodules using instrumented palpation, without a priori anatomical knowledge. The inverse analysis was applied to a model system, the human prostate, and was based on the reaction forces which can be obtained by trans‐rectal mechanical probing and those from an equivalent FE model, which was optimized iteratively, by minimizing an error function between the two cases, toward the target solution. The tumor nodule can be identified through its influence on the stress state of the prostate. The effectiveness of the proposed method was further verified using a realistic prostate model reconstructed from magnetic resonance (MR) images. The results show the proposed framework to be capable of characterizing the key geometrical indices of the prostate and identifying the presence of cancerous nodules. Therefore, it has potential, when combined with instrumented palpation, for primary diagnosis of prostate cancer, and, potentially, solid tumors in other types of soft tissue.

success of such devices, there is a need to develop the detailed solid mechanics of the tumor in situ, in order that improvements in sensitivity, specificity, and resolution can be achieved.
Many studies have concluded that cancerous prostate tissue has a higher elastic modulus than its healthy counterpart. 5,6 However, even equipped with properties for cancerous and healthy tissue, it is still challenging to use these for quantitative tumor identification and assessment. This is because the inclusion of a cancerous nodule with different mechanical properties from its surroundings will result in a disturbance to the stress distribution in the tissue, when it is subject to external loading, such as instrumented palpation with a prescribed probing depth. In principle, it is possible to identify the presence of tumor nodules by understanding how the mechanical response reflects the disturbance to the stress field introduced by the tumor nodule. A number of studies have addressed this "inverse problem" in the scenarios of liver, 7 myocardium, 8 prostate 9 and breast 10 tissues, and porcine liver and kidney. 11,12 In addition, studies by Kim et al 13 and Ahn et al 14 used a palpation device on phantoms with similar mechanical properties to a cancer nodule in prostate, comparing the results to FE calculations. However, the method required a knowledge of both geometry of the tissue samples and the elastic parameters for nodule and matrix. Similarly, the inverse model developed by Liu et al 15 for the estimation of nodule depth in prostate required a priori knowledge of the force feedback of and identical sample without a nodule, which may not be realized in practice. Notably, surface displacement as results of indentation was also used to inversely derive the distribution of shear modulus in soft tissue with stiff inclusions 16,17 and the results have shown promising capability of such methods even when high noise levels, although artificially created, were present. Other inverse methods, such as those based on artificial neural networks, [18][19][20] also require a priori knowledge of the nodule geometry or the stress distribution in the tissue.
This study aims to establish an FE-based inverse method to characterize the size of a soft tissue sample and detect the presence of a cancerous nodule in it, without a priori knowledge of the nodule geometry and location. Prostate tissue is used, without loss of generality, as an exemplar system. The method is based on knowing the reaction forces obtained from mechanically probing a prostate sample and those derived from an equivalent FE model, which is optimized iteratively toward the target solution by minimizing an error function between the two sets of reaction forces. By doing so, it is possible to estimate not only the volume of the prostate, but, more importantly, the embedded tumor nodule can also be identified from its influence on the stress field resulting from mechanical probing. The effectiveness of the proposed method is verified using a realistic prostate model reconstructed from a magnetic resonance (MR) image.

| Problem formulation
The prostate is surrounded by a number of other structures in vivo, including the bladder, and the rectum as well as being embedded in connective tissue, as illustrated by an MR image in Figure 1A. Digital Rectal Examination (DRE) involves accessing the palpable surface trans-rectally, which will deform these other structures to accommodate the probe before palpation actually commences. To avoid excessive complication, and to focus on the problem in hand, the FE model was simplified, as shown in Figure 1B, by introducing constraints at the boundaries representing the interactions with the F I G U R E 1 The schematic of the simplified "prostate" FE model. A, An MR image of the male pelvis, illustrating the anatomical features. The posterior surface can be probed through the rectum; B, The simplified FE model which represents the prostate tissue and its thickness values (h) at probing points. The constraints were based on the anatomical features shown in A. At all probing points (p), a displacement of depth (d) was applied along the anterior-posterior axis and the reaction force calculated connective tissue and bladder and other structures. The rectum wall is not considered in the model, as its effect has been proven to be negligible when high probing depths are used. Instrumented palpation was modeled using a rigid "probe," with prescribed displacement (d) along the anterior-posterior axis repeated at a number of points (p).

| Material parameters and FE models
The point-wise probing was considered to be quasi-static, using a strain rate lower than 0.01 second −1 , allowing the prostate tissue to be considered as a (hyper)-elastic material. 21 To allow for the high local strain that may occur in FE simulations, a neo-Hookean hyperelastic model was used for both cancerous and noncancerous tissue, with strain energy function expressed as: where I 1 is the first invariant of the right Cauchy-Green deformation tensor and J e the elastic volume ratio. The remaining material parameters, C 1 and D 1 , are related to bulk modulus (K 0 ), initial shear modulus (μ 0 ), and Poisson's ratio (ν), as shown in Equation (2). The prostate tissue was modeled as an nearly incompressible material, 22 and the Poisson's ratio of both tissue types is 0.49. Due to lack of published data of direct hyperelastic measurements for prostate tissues, mechanical properties of prostate tissue were adopted from the work by Hoyt et al, 23 who derived the Young's moduli of the noncancerous and cancerous tissues from ex vivo radical prostatectomy samples using Kelvin-Voigt Fractional Derivative (KVFD) model. After fitting the neo-Hookean model in Equation (1) with the elastic properties reported in Hoyt et al, 23 the noncancerous tissue (reported to have an elastic modulus of 17 kPa) had C 1 = 2.85 kPa and D 1 = 7067 kPa and the cancerous tissue (an elastic modulus of 42.5 kPa) had C 1 = 7.12 kPa and D 1 = 2827 kPa. For both materials, the parameter of D 1 is significantly greater than C 1 , indicating near incompressibility. The probe was modeled as a rigid material of diameter 10 mm, and its contact with the prostate tissue was assumed to be frictionless. 24 A quasi-static displacement along the anterior-posterior axis was applied to the probe, with depths ranging between 2 and 8 mm. The FE model was meshed with 3-node linear triangular elements (mesh refinement was conducted to ensure cost-effective convergence of FE simulations) and solved in ABAQUS (Dassault Systemes, Vlizy Villacoublay, France) using implicit quasi-static analysis. CPS6M elements (with hourglass control and full integration with selectively reduced integration) for 2D cases and C3D8R elements (reduced integration and hourglass control with additional mesh refinement) for 3D cases were used to ensure no hourglassing or volumetric locking in all models. Standard hourglass control method in ABAQUS was adopted, based on the classical method proposed by Flanagan and Belytschko, 25 where additional artificial stiffness values that are relatively small when compared to the actual stiffness of the material were introduced. This approach, usually accompanied by mesh refinement, can effectively achieve hourglass control in elements with reduced integration. 26

| Inverse FE method: Prostate volume estimation and tumor nodule identification
For the proposed inverse method, an FE model with "initial guess" of the prostate "thickness," h i , was established, as illustrated above in Figure 1B, and an iterative process was executed to drive the initial guess to the target (ie, the true solution), based on the difference between the reaction forces obtained from probing the prostate and the FE model being iteratively optimized. In this section, details of the proposed inverse method will be given, which aims to estimate the volume of the prostate and identify the presence of the cancerous nodule.

| Levenberg-Marquardt algorithm
To construct the optimization algorithm for estimating the prostate size, the Levenberg-Marquardt (LM) method was adopted. 27 It iteratively minimized an error function of reaction forces between the current FE model and the "target" model: which denotes the difference between the sums of target reaction forces (of "true" solution), RF Target, and those from the FE model being optimized, RF FE , at all probing points (p n ). The parameters, that is, the thicknesses (h), were those to be optimized against the first derivatives of the error functions. In the 2D formulation illustrated in Figure 1B, the Jacobian matrix (J) can be numerically calculated by perturbing each variable, h n , individually: where i represents the i th iteration, and F n , h n , and Δh n denote the error function, the thickness variable and the perturbations in the thickness variable at n th probing point, respectively. The maximum change of Δh was set to be 0.1 mm to ensure stability of the algorithm and accuracy of solving the first derivatives for the Jacobian matrix. The updated parameters, h i + 1 used in the following iteration step (i + 1) were obtained from: where H is the Hessian matrix (ie, the square matrix of second-order partial derivatives), I is the identity matrix, δ the increment of the variables h, and λ is the damping factor, which was used to ensure the robustness of the optimization algorithm and will be discussed in the following section. It is worth noting that the proposed optimization algorithm only estimates the size of the prostates at each probing points. As a result, the anterior surface of the prostate model is estimated with a cubic spline function using the thickness values, h i .

| Efficiency and accuracy of the LM algorithm
The damping factor, λ, is the key parameter that controls the efficiency and accuracy of the convergence of the LM algorithm. 28 If the error function in Equation (3) increases between two consecutive iterations, the updated parameter h i + 1 will not be accepted. As a result, the damping factor will increase in order to drive the algorithm toward the direction of gradient descent, which may lead to slower but guaranteed convergence to mimima. Otherwise, the damping factor is reduced leading to faster convergence. Thus, the algorithm is capable of alternating between a slow descent and a fast convergence, depending on the evolution of the error function. The initial value of the damping factor, λ, was defined as the maximum value among the diagonal elements in the initial Hessian matrix H, as it can be related to the eigenvalues. 29 The parameter that dictates if the damping factor is to be increased or reduced is the gain ratio, ρ, which is defined as the ratio between the actual reduction in the error and the predicted reduction in error, based on the chosen increment δ: whose absolute value is between 0 and 1, although it can be either positive (indicating a need to decrease the damping factor) or negative (to increase the damping factor): where the variable v was used to increase the value of the damping factor, of a multiplier of 2, when the gain ratio was negative.

| Inverse analysis procedure
The proposed procedure of inverse analysis is illustrated in Figure 2, aimed to minimize the difference between the target reaction forces (eg, from an instrumented palpation measurement) and their counterparts from the FE model, based on quasi-static point-wise palpation over the posterior surface of the prostate through the rectum. The iterative process is terminated when the convergence criterion is met, where the difference between two sets of reaction forces: is smaller than the convergence threshold, 0.1%. The essential premise of the technique is that the presence of tumor nodule will cause distortion of the stress field, hence different reaction forces at one or more probe points. Since the presence of a tumor nodule is not known at the time of probing, the prostate size estimated from the inverse analysis is either (a) the "true" size, if there is no tumor nodule present, or (b) a "distorted" size, if a tumor nodule is present. Therefore, tumor nodule identification requires a robust method to distinguish between these two cases.
Here, it is hypothesized, and later validated, that it is possible to identify the presence of the tumor nodule by analyzing the reaction force data using multiple probing depths. This hypothesis is based on the likely interplay between the probing depth and the existence of a tumor nodule, whereby any stress distortion is likely to be "amplified" by deeper probing. Therefore, by increasing the probing depth, one could examine the estimated area of the prostate with respect to increasing probing depth and conclude that, if the area estimation varies, a tumor nodule is present.

| RESULTS AND DISCUSSION
This section is divided in three areas. First, a sensitivity study is carried out for prostate volume determination against a target of a tumor-free gland of idealized (rectangular) shape. Next, the hypothesis of tumor detection using variable probe depth is tested in the rectangular shape. Finally, tumor detection is demonstrated in a real prostate shape with an idealized (although realistic) tumor.

| The framework of inverse analysis: Sensitivity study on tumor-free model
Two dimensions of sensitivity were tested in the idealized tumor-free model shown in Figure 1B; the effects of the initial guess of prostate size and number of probe points, and the effect of depth of probing.

| Initial guess and spacing of probing
As shown in Figure 2, an initial guess in the size of the prostate h 1 …h n is required to initiate the inverse algorithm, therefore the sensitivity of the algorithm to the initial guess needs to be explored. For simplicity, all target values of h are identical at 35 mm (ie, a rectangular target) and the size range of the initial guess was chosen to be between 10 and 60 mm, a little wider than that of adult prostate (20-50 mm). Probing to a 5 mm depth, at 3, 5, or 7 points (ie, 15, 6, 4.3 mm horizontal spacing) was carried out on the rectal surface of the model. F I G U R E 3 Estimation errors in prostate size for a 35 mm target, when the initial guess and the number of probing points are varied Figure 3A shows the relative errors in estimated areas of prostate using the proposed inverse FE method, as the initial guess is varied. Convergence is reached (with relative errors less than 0.1%) regardless of the initial guess chosen, although the rate of convergence is influenced by how far the initial guess is from the true solution, as might be expected. Figure 3B shows the effect of number of probing points, that is, spatial "resolution." Here, convergence is achieved in all three cases, although it is worth noting that the size of the Jacobians in Equation (4) increases with the power of the number of probe points with consequences for computational cost. Figure 4 shows the effect of varying the "target" thickness from 20 to 95 mm and the probing depth from 1 to 9 mm. Five probing points were used in all cases, with an initial thickness of 10 mm. All 56 cases are presented in Figure 4A, where the vertical bars illustrate the relative error in estimated area of prostate. In general, the larger the "target" prostate model is ( Figure 4B) or the greater the probing depth is ( Figure 4C), the higher the estimation error, although all cases reach convergence, most with under 1% error. Thus, the proposed method can estimate the total area of a tumorfree prostate, irrespective of the probing depth used.

| Identification of tumor nodule
In this section, the proposed methodology will be applied to cases where a tumor nodule is present in the model, in order to test the feasibility of nodule identification. Figure 5A shows a schematic of the model, where a tumor nodule F I G U R E 4 Relative errors in area estimation when the target size and probing depth are varied of diameter 15 mm and with Young's modulus a factor of 2.5 higher than the "matrix" is placed at a depth of 10.65 mm from the top (posterior) surface. This represents a readily palpable tumour 30 so is a reasonable test case for the ability to distinguish the effects of a local increase in stiffness from an increase in thickness in the idealized prostate shape.
Probing was carried out at three points (p 1 , p 2 , and p 3 ), using a probing depth of 8 mm. A comparison between the reaction force profiles of models with or without the tumor nodule, is shown in Figure 5B. The values of reaction forces are significantly higher when the stiff tumor nodule is present. This is also reflected in the stress distributions under probing at location p 2 ( Figure 5C,D).
Using the same models as Figure 5A as the "target," the method was applied with a second probing depth of 2 mm. The results for the four cases, that is, with and without nodule and 2 and 8 mm probing depths, are illustrated in Figure 6. When no tumor nodule is present in the target, Figure 6A, the estimated prostate has a much flatter profile at the anterior surface than when a nodule is present ( Figure 6B). As hypothesized, increasing the probing depth amplifies the effect of the tumor nodule on the stress field, leading to an even greater underestimation in prostate area. This is illustrated by the error in area estimation illustrated in Figure 6C. In tumor-free cases, when the probing depth changes from 2 to 8 mm, the estimated area of prostate does not change significantly (errors of 0.08% and 0.06%, respectively). However, when the tumor nodule is present, the error in estimated area of prostate increases drastically, and even more so when a deeper probing is used. This demonstrates the potential of using the proposed inverse method for tumor nodule identification by using a range of probing depths.
It should be noted here that the prostate model in this study is simplified to have two distinct tissue types "mechanically," that is, "noncancerous" and "cancerous." Noncancerous tissues often include healthy tissue, benign prostate hyperplasia (BPH), etc., whilst cancerous tissue may be at different stages of malignancy. There has been a very limited number of published papers, including one from the author group, 3 that showed BPH often does not  (3), was normalized using its value at the first iteration. Thickness of initial guess: 15 mm significantly increase the tissue elasticity in quasi-static measurements. For cancerous tissue at various stages, it is reasonable to assume the mechanical properties could vary. Therefore, the material property of tumor nodule was reduced (soft-equivalent to 25.5 kPa) and increased (hard-equivalent to 59.5 kPa), by 40%, from the benchmark value (equivalent to 42.5 kPa), to explore the sensitivity of the proposed method to the nodule stiffness. The results, as shown in Figure 6D, indicated that the nodule stiffness could indeed influence nodule detectability. When the nodule became softer, that is, with a Young's modulus of 25.5 kPa, in contrast to17 kPa for noncancerous tissue, probing at a depth of 2 mm led to insufficient sensitivity (4.2% difference in prediction error). However, using the probing depth of 8 mm could significantly improve the sensitivity, in both cases of "softer" and "harder" nodules.
Furthermore, the nodule depth from the posterior surface was also investigated. The results shown in Figure 6E demonstrated, rather intuitively, that the greater the nodule depth, the poorer the detectability becomes. These results suggest a limit on tumor detectability with respect to the nodule depth and adopting a great probing depth would significantly improve such a limit and the sensitivity in nodule detectability.
In order to further demonstrate the capability of the proposed method, a simplified 3D model was created, as illustrated in Figure 7A. It should be highlighted that, due to the nature of the optimization algorithm proposed, the method seeks to minimize the reaction forces between the "target" and the FE results of the "model" being iteratively optimized, regardless of whether the problem itself is in 2D or 3D. Material properties and the size/depth of the tumor nodule remained the same as the 2D case shown in Figure 5A for comparison. Four probing points, using a probing depth of 8 mm, were deliberately placed in an array of 2 × 2, away from the nodule at the center. Nevertheless, similar results to those from the 2D models were obtained in Figure 7B-when the target values of reaction forces were reached, the inclusion of the stiff tumor nodule could lead to significant "error" in predicting the volume of the prostate, whilst the error of volume prediction in the "healthy" case remains low. Such a contrast in volume predictions could indicate the presence of the tumor nodule.

| Prostate characterization and tumor identification: A feasibility study
To demonstrate the feasibility of the proposed framework in a realistic situation, a prostate model was reconstructed (using Scan-IP, Simpleware, Mountain View, California) from the MR image shown in Figure 8A and the same palpable nodule as used in Section 3.1.2 superimposed onto the model ( Figure 8B). The boundary conditions were applied as before, considering the anatomical features surrounding the prostate. Probing was performed vertically at three locations on the posterior surface, as indicated in Figure 8B, using probing depths ranging from 4 to 8 mm with and without the palpable tumor. This arrangement is a reasonable facsimile of actual patients on whom instrumented palpation measurements and FE simulations have been made in vivo. 31 Figure 9 summarizes the estimated areas of the MRI-based prostate model for the five probing depths, with and without the tumor nodule. In the tumor-free models (in black), the estimated areas of prostate are close to the target with little variation between probing depths. By contrast, even when the probing depth is low the presence of tumor nodule leads to underestimated prostate sizes, worsening as the probing depth increases. It is also interesting to note that most of the converged models (even those where a tumor is present) reflect the asymmetry of the target (h 1 > h 3 ). Most importantly, Figure 9 demonstrates that, when assessing the response of prostate tissue to probing from the posterior side, one can vary the probing depth and use the obtained reaction forces to estimate the "apparent" prostate size using the proposed method. If the estimated prostate size remains relatively constant with little variation when the probing depth is varied, the examined prostate sample can be regarded as homogeneous, meaning that there is no palpable lesion. Otherwise, the consistent variation in the estimation of prostate size, due to the disturbance in the stress field caused by the presence of a stiff nodule, could serve as a primary indicator of a tumor.

| CONCLUSIONS
There is a need of robust and effective methods for the quantitative interpretation of data from mechanical probing in order to identify inhomegenities in soft tissue. This work has focused on the identification of tumor nodules in the prostate by instrumented probing.
A finite-element based inverse method has been proposed to estimate the size of the prostate using reaction force data at several probe points on the posterior surface of the prostate. A sensitivity analysis showed that the initial guess of the shape and the interval between probing points have little influence on the estimated prostate size, showing the robustness of the proposed method.
Importantly, it was hypothesized, and later validated, that the presence of tumor nodule that is of higher modulus can be identified by conducting probing at various depths, without a priori anatomical knowledge of the prostate such as volume and geometry. This was possible because the estimate of prostate size changes significantly as the probing depth is varied if a tumor is present, but remains constant when no tumor nodule is present. The sensitivity to the material property and depth of tumor nodule was also explored and their influences characterized.
The accuracy of the results is not affected by the initial conditions which means that the inverse procedure can be applied irrespective of the size and geometry of the organ under investigation. Increasing the number of indentation points would improve the quality of the diagnosis, although it is necessary to explore the accuracy and sensitivity of the method to features of the cancer nodule, such as its diameter, depth, and the uncertainty of the probing location when acquiring the force data from patients in vivo. More importantly, although the examples given in this study demonstrated the feasibility of the proposed method where a single, large, nodule is present in prostate, the presence of multiple nodules is likely to have certain influence. This will be investigated in the future work, where soft computing methods such as fuzzy logic and artificial neural networks could be employed to help quantify the uncertainty given rise by the presence of multiple nodules toward a more robust inverse algorithm for detection of tumor nodules. In addition, boundary conditions at the anterior surface are simplified in this study to represent interaction with the surrounding tissues. Despite the previous study by the author group 31 where the patient-specific modeling of pelvis was introduced to in vivo instrumented palpation, more accurate representation of the interaction between prostate and surrounding tissues need to be further investigated. Furthermore, the current framework assumes two "materials" (cancer and matrix) to be homogeneous and isotropic. The mechanical characterization of soft tissue may be more complex in reality 32 and such effects on the proposed method will be investigated in future studies.