Quantitative diagnostics of soft tissue through viscoelastic characterization using time-based instrumented palpation

Although palpation has been successfully employed for centuries to assess soft tissue quality, it is a subjective test


Introduction
Viscoelasticity in biological tissue was first recognized by Bayliss and Robertson (1939) in cat lungs. The time-dependent behavior in tissue has been attributed to the movement of fluid within its microstructure as in cartilage or liver (Huang et al., 2001;Kerdok et al., 2006) and also to the lubricating effect that the proteoglycan matrix exerts between the collagen fibrils that occur for example in blood vessels (Holzapfel, 2000). The mechanical behavior of collagen is also known to vary with strain rate (Bischoff, 2006), which is believed to contribute to the viscoelastic response of tissue such as lung (Ebihara et al., 2000) Chen). (Garteiser et al., 2012), prostate (Phipps et al., 2005) and skin (Xu et al., 2007). Breast tissue has been observed to exhibit viscoelasticity which varies with conditions such as hormonal balance and disease (Balleyguier et al., 2013;Chen et al., 2013;Sack et al., 2013;Sinkus et al., 2005). The viscoelastic behavior of brain tissue (Guo et al., 2013;Rashid et al., 2014), has been shown to vary with pathological conditions such as multiple sclerosis (Streitberger et al., 2012), cancer and ageing (Sack et al., 2009). Liver tissue exhibits viscoelasticity (Kerdok et al., 2006;Taylor et al., 2002) that has been observed to vary, subject to diseases such as fibrosis (Salameh et al., 2007), water content (Kerdok et al., 2006) and preservation time (Yarpuzlu et al., 2014). Arteries have a rate-dependent behavior where the response is often anisotropic and nonlinear . In cartilage, viscoelasticity is a desirable property allowing the stress associated with physiological work to be less damaging (Soltz and Ateshian, 1998). Muscular tissue also exhibits anisotropic viscoelasticity that becomes nonlinear at large strains (Takaza et al., 2013;Van Loocke et al., 2009). Thus, there is huge potential to use viscoelastic measurements in the assessment of the quality of soft tissue, in terms of a number of positive or negative physiological or pathological axes.
The elastic behavior of biological tissue has been thoroughly discussed in reviews (Ali et al., 2010;Beda, 2007;Boyce and Arruda, 2000) and can be described using such models as the 8-chain model (Bergstrom and Boyce, 2001) or other phenomenological continuum models such as Mooney-Rivlin (Palacio et al., 2013) or Ogden (Rashid et al., 2014). Viscoelasticity can be characterized using Prony series (Park and Schapery, 1999), Kelvin-Voigt fractional derivative (Taylor et al., 2002) and recruitment models (Bates, 2007). Of particular note is the work of Holzapfel and co-workers (Gasser and Holzapfel, 2002;Holzapfel, 2000;, who proposed viscoelastic models based on a strain energy density function of fiberreinforced anisotropy that predicts the viscoelasticity of soft tissues such as arterial walls. Peña et al. (2010) also modeled the viscoelastic behavior of vascular tissue including the softening due to damage. Anisotropic viscoelastic models have also been proposed to predict the behavior of heart valves (Anssari-Benam et al., 2011). Bergström and Boyce (2001Boyce ( , 1997 described tissue viscoelasticity using the interaction between two networks, one related to the elasticity and the other to the viscous behavior. Although both networks have the same deformation gradient, the velocity gradient in the viscous network is decomposed into elastic and viscous velocity tensors. The so-called recruitment models consist of multiple Maxwell elements arranged so that, when a given element reaches the allowable strain, the following element starts to deform. Such models have proven successful in lung tissue, for example, where stress relaxation follows a power law (Bates, 2007).
The characterization of viscoelasticity for biological tissue requires the entire history of deformation in order to take into account its rate-dependency. Two different approaches are often adopted: creep/stress relaxation (quasi-static/lowrate) and dynamic tests. The results from quasi-static tests are often correlated by means such as Prony series (Brinson and Brinson, 2008), although other fitting techniques have been used, including a more generalized transfer function (between stress and strain) in Laplace space (Yu and Haddad, 1994), or using splines (Tran et al., 2011). These tests incidentally provide the instantaneous and the long-term moduli. In dynamic tests, the system (usually in a tension-compression configuration) is subjected to a pre-defined displacement function, such as a sinusoid. The amplitude of the force response and phase lag between the displacement and force signals can then be measured over a range of frequencies. In such tests, different techniques can be used to identify the corresponding transfer functions to relate the stress and strain in Laplace space (Renaud et al., 2011;Yu and Haddad, 1994).
It should be noted here that the two types of tests in fact measure the same thing, where the transfer function obtained from a stress relaxation/creep test mathematically dictates the dynamic behavior that would be measured in the frequency test, and vice versa. Therefore, the ranges of time and frequency that are used in any test may influence the results that are measured. For example, using higher frequencies or shorter times allows the detection of the smaller time constants that define the behavior during the early stage of stress relaxation, whereas lower frequencies provide information on the larger time constants that characterize the long-term response.
The time constants obtained from models such as Prony series define the characteristics of amplitude and phase-shift, but overfitting of a single data set using excess time constants can lead to undesired, sometimes entirely artificial, characteristics at certain frequency ranges. Therefore, it is critical to find the balance between the fitting error and the number of fitting parameters. Ideally the results of both dynamic and quasi-static tests should be fitted at the same time to avoid overfitting; however this would increase the complexity which becomes inapplicable for clinical tissue diagnostics.
Various types of medical devices have been reported for measurement of the mechanical properties of soft tissues. Oflaz and Baran (2014) designed a hand-held instrument to measure the in situ stiffness of biological tissues using a constant depth of indentation. A similar approach was proposed by , where the system remained attached to a movable cart. Indentation has also been used to measure the force-displacement behavior of tissues (Carson et al., 2011;Phipps et al., 2005). A less conventional approach was used by Ahn et al. (2012) where a device that mimics the sweeping palpation of the finger was developed. Sangpradit et al. (2011) presented a device for rolling tissue indentation again for application to tissue diagnostics. Some types of tissue have been adopted as candidates for diagnosis using direct or indirect mechanical palpation. Carson et al. (2011) and  used hemispherical indenters to measure the prostate stiffness ex vivo, which was then correlated to tissue quality. Mechanical approaches were also used to characterize swine brain in vivo (Miller et al., 2000), porcine  and human (Carter et al., 2001) liver and skin, as well as soft tissue from lower limb extremities (Pailler-Mattei et al., 2008; Silver-Thorn, 1999).
Techniques for characterizing tissue in vivo have seen significant improvement during the last two decades, for example in prostate cancer diagnostics. Prior to the 1980s, most prostate biopsies were performed using a needle which was directly guided by the clinician (Sartor et al., 2008), largely relying on their experience of the anatomy. In the early 90s the technique was improved by the introduction of transrectal ultrasound (TRUS) guided biopsies and the sextant approach. The procedure is often repeated 10 to 12 times (up to 20) in different areas of the prostate, the samples then undergoing histological examination. However, despite these improvements, such an approach has some disadvantages such as false negatives (as the TRUS does not identify prostate cancer regions and the needle only samples small cores), pain and discomfort for the patient and the possibility of complications including rectal bleeding and hematuria Su (2010). The elastic behavior of the prostate is regarded as an important index in clinical diagnostics (Ahn et al., 2011;Carson et al., 2011;Hoyt et al., 2008;Masuzaki et al., 2007;Murayama et al., 2007), whereas its viscoelasticity has received less attention. Krouskop et al. (1998) reported a phase shift between stress and strain when prostatic tissue was dynamically examined at 0.1, 1 and 4 Hz. Carson et al. (2011) also showed that, when prostatic tissue is subjected to mechanical indentation, the required force depends on the velocity of indentation. Other studies by Phipps et al. (2005) and Hoyt et al. (2008) also indicated that the prostatic tissue (both healthy and diseased) exhibits phase shift and stress relaxation and therefore concluded that viscoelasticity could be of great importance in clinical assessment of tissue pathological conditions in the prostate.
Therefore, for prostate disease diagnosis, there is a clear need for an inexpensive, less invasive, fast and reliable method to improve quantitative assessment of tissue quality. In the more general framework of tissue diagnostics there is also an obvious need for quantitative methods through which the mechanical properties of the tissue can be related to its pathological condition. Some diagnostic techniques, such as elastography, already exploit the fact that relative change in mechanical properties is a good indicator of tissue quality, which could form the basis of an index for quantitative tissue diagnostics and enhance other diagnostic procedures, especially for in vivo recognition of cancerous nodules and pathological analysis.
This paper presents a framework for characterizing tissue viscoelastic properties as a diagnostic index of its pathological condition. The framework is applied to both healthy and cancerous prostates subject to dynamic mechanical palpation, with a view to forming general conclusions on the assessment of tissue condition and the identification of diagnostic parameters of cancerous nodules in soft tissue. As mentioned above, instead of requiring parameters from different viscoelastic material models, a set of generalized and simplified parameters is adopted for quantitative assessment of the tissue without a priori knowledge of its histology and pathology. The methodology and analysis of sensitivity to probe deployment options are explored by simulation in both 1D and 2D prostate models. The proposed method is then applied to a 3D model obtained from an MRI scan of an actual excised prostate.

2.
Materials and methods

Viscoelastic models of biological tissue
The Zener model is one of the simplest rheological models for characterizing both creep and stress relaxation of viscoelastic materials and is adopted here to demonstrate the usefulness of time constants as indicators of rate-dependent behavior for quantitative tissue diagnostics. The model consists of a spring in parallel with a spring-dashpot series couple (Klatt et al., 2007;Renaud et al., 2011), and its stress-strain relationship can be expressed as where E 1 denotes the stiffness of the single spring, and E 2 and η the stiffness and damping coefficient of the spring-dashpot couple, respectively. When the system is subjected to a uniaxial dynamic force F¼sin(w Á t), the stress σ is sin(w Á t)/A, where A is the area, w the angular frequency and t the time. Eq. (1) then becomes To simplify the notation: and K 1 is to be obtained from boundary conditions. It is important to note that, if a dynamic displacement is used as the input, then the resulting differential equation is symmetrical, i.e. the typology of Eq.
(2) is the same except the values of a, b, c and s therefore the solution for stress will be of the same form as Eq. (3). Using Prony series, the transfer function of a viscoelastic material can be described as: which indicate that, to describe the behavior of a viscoelastic material, it is necessary to know the stress history which the material has experienced.

Quasi-static tissue characterization
The mechanical properties of prostatic tissue are known to change in the presence of benign prostatic hyperplasia (BPH) and cancer (Krouskop et al., 1998). For the sake of simplicity, mechanical properties of normal prostatic and cancerous tissue obtained from stress relaxation tests by Hoyt et al. (2008) are adopted in this paper. The relaxation curve is fitted by normalized Prony series in Fig. 1, using least squares approximation with a sufficient number of points based on the Kelvin Voigt Fractional Derivative (KVFD) model  (Hoyt et al., 2008). The least number of terms in Eq. (5) to fit the data adequately was found to be two; using a single term only allows fitting to either short or long time but not both. Table 1 shows the corresponding viscoelastic parameters, where D i and τ i are the constants to be determined in Eq. (5).

Dynamic tissue characterization
Dynamic characterization of viscoelastic tissue is usually performed by constructing phase and amplitude diagrams with respect to a range of frequencies and obtaining the parameters using a certain model (Martinez-Agirre and Elejabarrieta, 2011; Renaud et al., 2011). Alternatively, the storage and loss moduli that relate to the elastic and viscoelastic components can also be obtained in the same way (Madsen et al., 2005;Nayar et al., 2012). In its simplest form, mechanical palpation involves imposing an input displacement from an indenter and measuring the corresponding reaction force. Without loss of generality, a sinusoidal displacement with a compressive mean position is chosen so that separation between tissue and indenter can be minimized. A smoothed mean value of the force signal f(t) is obtained using a weighted local regression algorithm (LOESS) (Cleveland, 1979). By fitting f with a series of exponential functions as described below, a set of parameters: The parameters 'a' and 'b' are related to the material elasticity, b, in particular, being an index for the long-term elastic modulus. It should be mentioned here that the duration of the palpation experiment is critical in determining the reduced set of parameters. From Eq. (6) it can be seen that, for a given time of palpation, one time constant (t i ) will predominantly influence the viscous behavior of the material. Time parameters significantly higher than t i (i.e. n Ziþ1) can be considered constants since their behaviors during short testing times are indistinguishable from an instantaneous modulus. On the other hand, smaller time parameters (i.e. n riÀ 1) will make the exponential term tend to unity which means that those terms will be added to the long-term modulus b i . In both cases, b becomes an important index for apparent elasticity of the tissue.

2.4.
Implementation of viscoelastic models

1D mathematical model
The 1D model was created first to conduct a theoretical analysis that allows a better understanding of tissue viscoelastic behavior in both transient and steady state. A sinusoidal displacement was applied at one end of the model, which is fully fixed at the other. The cross-sectional area is A, Young's modulus E and the length L. For the sake of simplicity, it is assumed in this section that the viscoelasticity of prostatic tissue can be modeled by a Prony series with one time constant, which, after normalization, is given by: When subjected to a certain displacement, the reaction force is described as Once the steady state is achieved the mechanical behavior is then governed by two oscillatory terms: a sine and a cosine, the latter of which accounts for the phase shift between input and output signals. The amplitude of the oscillations in steady state can be calculated as: Fig. 1 -Prony series approximation of the normalized stress relaxation obtained from the model (Hoyt et al., 2008), using one or two terms. It can be observed that results using two terms present better fitting for both short and long term behaviors simultaneously.  where ψ and ξ are the coefficients of the sine and cosine terms in Eq. (8) respectively, normalized with respect to E, A and L. The phase shift is calculated as 2.4.2. Prostate cancer diagnostics: 2D and 3D models Finite element models were created to demonstrate the effectiveness of the proposed method in quantitative tissue diagnostics. A 2D model, representing a cross-section of ex vivo prostate sample, which is constrained at the bottom surface (opposite the indenter), was created first. An indenter with a hemispherical tip was modeled as a rigid solid, to which a sinusoidal displacement was applied. The prostate was then modeled in ABAQUS (Dassault Systemes, Vlizy-Villacoublay, France) and the material behavior was adjusted with a second order Ogden strain energy density function to mimic a linear viscoelastic material with an elastic modulus of 17 kPa for healthy tissue and 34 kPa for cancerous tissue, according to Table 1. The mesh under the indenter was refined to allow a better solution to the contact problem and inertial effects were taken into account using a tissue density of 1 kg/m 3 (Torlakovic et al., 2005). A 3D prostate model was also created to provide a more realistic evaluation of the proposed methodology. The model was based on a prostate specimen, excised using the laparoscopic radical prostatectomy approach, from a patient undergoing surgery for localized prostate cancer. After removal of the prostate, a 7-T magnetic resonance imaging (MRI) was performed on the fresh specimen, using a resolution of 0.31 mm in the sagittal and coronal planes and 1.5 mm in the axial plane. In total, 44 images were obtained to reconstruct the 3D model using Scan-IP (Simpleware Ltd., Exeter, UK). Fig. 2 summarizes the proposed methods for carrying out such a procedure in practice, where Fig. 2(a) shows the workflow to assess tissue quality from the MR imaging via viscoelastic characterization and Fig. 2(b) illustrates the setup of in situ palpation of the reconstructed prostate model in comparison with the ex vivo measurement.

3.
Results and discussion

Tissue characterization: 1D theoretical analysis
Clinical tissue diagnostics is a major challenge due to various complex constraints such as time, available number of samples, patient discomfort and pathological conditions. Therefore, it is critical to select the optimal measures (such as phase-shift and/or amplitude), and to choose the optimal parameters, such as frequency and number of measurement points. The 1D dynamic analysis was carried out using different time parameters across the available frequency scope in order to identify the optimal range of frequency over which a dynamic measurement should be performed. The range of frequency starts at 1 Hz, since lower values would lead to excessively long tests. The upper limit is determined by those which produce undesirable dynamic behavior and, in the extreme, data loss due to indenter lift-off if the sample is unable to recover its original shape when the indenter retracts. Fig. 3 shows the amplitude and phase response from a sinusoidal displacement over a range of frequency using various viscoelastic parameters. It can be seen in Fig. 3(a) that viscosity plays very little role at low frequencies because the low velocity of deformation leads to the strain being effectively time-independent. Under such conditions the behavior can be characterized using the long-term modulus. Increasing the strain rate (frequency) amplifies the force amplitude and also the apparent stiffness until a plateau is reached. Materials with time constants  τ in the range of 1 to 10 s exhibit primarily the instantaneous modulus when subjected to frequencies in excess of 1 Hz. The range of frequency where viscoelastic behavior of such materials exhibits is found to be between 0.001 and 1 Hz, as shown in Fig. 3. Therefore, the differentiation of the tissues which exhibit time constants consistent with the measurements of Hoyt et al. (2008) ðτ Healthy ¼ 5:49 s and τ Cancer ¼ 13:87 sÞ would need to be based on their (quasi-)elastic behavior since frequencies lower than 1 Hz are not clinically practical. Such quasi-elastic behavior can also be observed in the phase lag diagram shown in Fig. 3(b), which tends to zero in the range of higher frequency (41 Hz). To differentiate materials using phase lag their time constants should be at least one order of magnitude different. More importantly, to observe a noticeable phase shift, the material needs to have at least one time constant in the range of 0.001 s to 1 s. Fig. 4 shows the amplitude and phase lag behavior of the prostatic and cancerous tissue when the dynamic behavior is fitted using a Zener model. The parameters of the Zener model are estimated by least squares fitting, using the quasi steadystate solutions to the reaction force from input (displacement) frequency at 0.1, 1, 10 and 100 Hz for both Prony and Zener models. Fitting the results, which are shown in Table 2, over a wide range of frequencies allows characterization of the dynamic behavior with a reduced risk of overfitting.
The viscoelastic behavior predicted by the Zener model is similar to that obtained using a Prony series, where a quasielastic response is evident at the extremes of high and low frequency illustrated by constant amplitude and negligible phase lag. The transition between the two states can be considered as the frequency range where dynamic measurements of viscoelastic properties can be made. The parameter that dictates this behavior is the time constant τ¼ η/E 2 , so it can be seen that the dynamically sensitive range runs between 0.001 Hz and 0.1 Hz for the case of amplitude and 0.001 and 1 Hz for the phase lag.
In this particular case, phase shift only becomes significant at frequencies below about 1 Hz when parameters from Hoyt et al. are adopted. In that case, unfortunately, the reaction force obtained will only depend on the long-term modulus and therefore only elasticity will be measured, losing most of the viscous information that the material presents.

3.2.
Sensitivity analysis: 2D dynamic study In a practical application of tissue diagnostics using mechanical probing, it would be rather difficult to perform palpation across the entire surface. Therefore a grid of test points would normally be used to assess the accessible surface. It is useful to explore the sensitivity of the proposed method with respect to various diagnostic indices, such as position, depth and size of a cancerous nodule. In this section the proposed method is applied to a 2D model with a representative cross-section of prostate sample where size, depth and position of a cancerous nodule vary. The displacement is applied to the prostate through an indenter, shown in Fig. 5 where stress distribution at different stages of indentation are illustrated, including initial contact in Fig. 5(b), mean position of indentation in Fig. 5(c) and lowest indentation point in Fig. 5(d). The stress becomes higher underneath the indenter and around the cancerous nodule when the indentation progresses, which leads to higher stress in a larger area and hence a larger reaction force, although there will be a limit to how deep an indentation can be made before discomfort or damage occurs. Increasing the mean depth of indentation for a given nodule size and position also increases the measured apparent stiffness, b, due to the increased strain and also the larger contact area of the indenter. Finally, excessively deep indentations cause an artificial increase in the measured stiffness due to the mechanical support at the bottom. Fig. 6 illustrates the changes in viscoelastic parameters a, b and τ of a prostate when two different sizes of nodule (large, 6 mm in radius, and medium, 4 mm in radius) both located 12.5 mm under the surface are shifted horizontally in 5 mm increments away from the indenter. As shown, there is little discrimination when the nodule is far away from the indenter, which is to be expected since the disruption to the stress distribution caused by a distant nodule is small. Therefore, in order to assess tissue quality, a sufficient number of indentations needs to be used for lateral resolution/discrimination.
The effects of size and depth of the cancerous nodule are explored in Fig. 7, and found to have a substantial effect on both parameters a and b. The depth and size of a cancerous nodule are, in fact, coupled in the elastic response, which means that a smaller nodule close to the surface will lead to similar force feedback to a larger/deeper nodule. Therefore, in this case, the elasticity alone is unable to discriminate between size and depth, although changes in τ (hereinafter referred to as Tissue Quality Index, TQI) remains less affected in both cases, which implies that using relative changes of both tissue elastic and viscous behaviors may lead to decoupling of depth and size of cancerous nodule, thus making quantitative diagnostics possible.

Quantitative cancer diagnostics in prostate: A 3D study
In this section the proposed method will be applied to the 3D prostate model shown in Fig. 8(a), which is reconstructed from an excised prostate with the cancerous nodule closer to  the posterier side. The palpation area used in the model was on the posterior side of the prostate and was sequentially indented at 8 sites as shown schematically in Fig. 8(b). The indentation sites were selected to maximize the area to be tested across the entire surface whilst maintaining good contact between indenter and prostate. The anterior surface was constrained to mimic the ex vivo boundary conditions during an actual mechanical indentation. Fig. 8(c) and (d) show the recorded force in both healthy and cancerous samples and the fitted smoothed data used to determine the reduced set of parameters for tissue quality assessment.
It is evident that the displacement and force signals are hardly distinguishable and the difference becomes negligible especially during steady-state. This again indicates that only using phase shift and/or amplitude may not be sufficient to assess the tissue quality in a quantitative way. In order to examine the applicability of the proposed method in 3D, two models are created here: one with a cancerous nodule inside which has the same viscoelastic parameters used earlier; another one which is fully healthy (material properties of healthy tissue are assigned to the cancerous nodule). The resulting viscoelastic parameters (i.e. a, b and τ) of the entire prostate are shown for both cases in Table 3. It is important to note here that, even in the healthy case, the reduced set of parameters does not remain constant at all probe points due to the surface curvature which changes the contact area and hence tissue response in the tissue-indenter contact zone, and therefore the force feedback. Nevertheless, at probe points 5, 6 and, 8, close to the nodule, the viscoelastic properties cause a stiffer and slower response for the prostate with the cancerous nodule. It is important to note that, as shown in Fig. 9, the tissue quality index offers a unique capability of quantitatively characterizing the location and extent of cancerous nodules inside the tissue.

Concluding remarks
Based on the evidence from the literature that the viscoelasticity of soft tissue is affected by some pathological changes such as cancer, this paper has used published data for viscoelasticity of cancerous and healthy prostatic tissue to establish a novel framework for quantitative diagnostics for soft tissue using dynamic palpation.
The proposed method features diagnostic procedures that can be used to obtain elastic and viscous behaviors simultaneously, with the viscoelastic parameters being able to characterize the cancerous nodule, thus becoming a more reliable index for quantitative assessment of tissue quality. The method is illustrated on a 3D prostate model reconstructed from an MRI scan of an excised prostate. It is shown that the change of viscoelastic time constant between healthy and cancerous tissue, which defines the tissue quality index, could be a key indicator for quantitative diagnostics of tissue pathological conditions. Furthermore, the method presented here has certain advantages in a clinical context such as reduced duration of examination and less invasiveness. The   proposed method is not limited to any particular material model or scale, and thus becomes useful for tissue where defining a unique time constant is not trivial. It is important to mention that the objective of this methodology is not to provide a very detailed description of the viscoelastic behavior. For that purpose, frequency testing as well as stress relaxation should be undertaken. Instead this study aims at providing a methodology to assess tissue quality in the framework of tissue diagnostics. Assessment of other pathological conditions, such as benign prostate hyperplasia (BPH) and prostatitis in prostate tissue, will depend on their elastic properties and also the difference in viscoelastic time constant in comparison to healthy tissue. Selection of the indenter size and shape as well as indentation depth also requires further investigation to maximize the capability to assess tissue quality with enough resolution whilst still satisfying clinical constraints.