A numerical method to detect soft tissue injuries from tissue displacements

Abstract

A method is introduced to detect soft tissue injuries from internal measurements of displacements taken while a quasi-static load is applied to the tissue. We aimed to quantify and improve reliability of the diagnosis of soft tissue injuries with quasi-static ultrasound elastography. Our method consists of (i) building a model of displacements in healthy tissues and (ii) characterizing injured tissue for how displacements estimated in this tissue differed from the model. The model was built by applying principal component analysis to a data-set of displacements in healthy tissues. A Monte Carlo simulation of Finite Element models was performed to test our method. Our method could detect injuries quantitatively (area under ROC curves equal to $0.86$ with no displacement noise, $0.85$ with $2\%$ noise and $0.81$ with $4\%$ noise) and provided maps with significantly higher contrast-to-noise ratios between damaged and healthy regions than maps provided by current methods, either strain-based and modulus-based ($p<0.001$).

Keywords:

• soft tissue injury
• principal component analysis
• damage displacement
• ultrasound elastography

AMS Subject Classifications:

• 65N21
• 92C10
• 92C55
• 74S30

1 Introduction

Ultrasound elastography (UE) was developed as an alternative to palpation [1, 2] and has emerged as a powerful diagnostic tool for diseases such as cancer [3] and hepatic fibrosis [4]. One approach to UE, known as quasi-static UE, estimates displacements inside the tissue with ultrasound during quasi-static mechanical loading. Estimated displacements have been used to examine injuries in soft tissues structures such as the Achilles tendon,[5–11] plantar fascia,[12, 13] supraspinatus,[10] and tendinous insertions in the elbow.[14] An injury is broadly considered here to be trauma that damages the structure of soft tissue.

With quasi-static UE, a natural question is how can the estimated displacements be used to detect injuries. Displacements are typically used to evaluate mechanical stiffness, which refers to a material’s resistance to loading, and mechanical stiffness is a natural candidate to evaluate, because structural damage alters mechanical stiffness. However, estimated displacements may not uniquely determine mechanical stiffness without some knowledge of tissue stress, boundary conditions, body forces and constitutive equations.[15] Estimated displacements may not even determine strain precisely, because displacement are subject to noise and limited resolution, and are estimated in one or two directions, if not just the axial direction, i.e. the direction of ultrasound wave propagation.[2] In short, mechanical stiffness cannot be precisely determined with limited displacement information and without stress and other mechanical information.

Because estimated displacements alone cannot determine mechanical stiffness, it is common for methods based on quasi-static UE to include mechanical assumptions about the soft tissue. Common methods fall into two categories: strain-based and modulus-based. On the one hand, strain-based methods use strain as a relative indicator of mechanical stiffness. That is, tissue is considered stiffer in a region with lower strain than in a region with higher strain. Consequently, these strain-based methods implicitly assume that tissue stress is approximately the same throughout the tissue.[16] On the other hand, modulus-based methods aim to estimate specific material properties, such as shear modulus by assuming at least one of the following properties: material isotropy, linearity, plane strain and plane stress.[16–21]

There are two main issues with strain-based and modulus-based methods. First, they frequently use inaccurate assumptions about soft tissue. In contrast with the assumptions described above, soft tissues have non-linear material behaviour and soft tissues that are also fibrous have anisotropic material behaviour.[22] Moreover, it is unclear whether certain soft tissues possess geometries that support uniform stress, plane stress or plane strain. For example, plane stress (plane strain) is assumed when the structure has one dimension that is an order of magnitude smaller (larger) than the other two dimensions,[15] but many structures do not satisfy this assumption. Even after assuming any one of the material properties listed above, material properties are not uniquely determined by displacements estimated in one or two dimensions.[23]

The second issue is the lack of quantitative methods based on quasi-static UE to detect injuries from strain or modulus.[24] The lack of quantitative approaches is chiefly due to poor repeatability and reproducibility, whereas measurements (e.g. strain) are difficult to repeat and reproduce, because quasi-static UE involves variable imaging windows, loading conditions and image processing algorithms. Particularly, imaging and loading are often conducted freehand, potentially leading to high variability in examinations. To address the lack of quantitative approaches, researchers aim to standardize examinations with, for example, automatic loading devices [25] or real-time feedback on loading (as offered by many commercial ultrasound systems). They have also attempted to minimize variation among examinations using the ratio of strain in one region and strain in another, rather than strain itself.[7] Even if loading, imaging and image processing algorithms were standardized, there remains variability in mechanical stiffness and anatomy among healthy volunteers. How, then, can we develop a quantitative approach when examinations are also subject to this variability?

In what follows, we propose an alternative strategy to overcome the two issues described above (Section  2). We focus on the case when a specific protocol is applied to a specific soft tissue structure. Our strategy is to supplement information provided by estimated tissue displacements with a data-set acquired prior to the examination, which consists of tissue displacements estimated in healthy tissues. We perform the statistical technique called principal component analysis (PCA) on the data-set to build a model of tissue displacements that arise from healthy tissues. We then introduce a measurement we call damage displacement to measure the difference between the model and tissue displacements measured during the examination. We test our method numerically using Monte Carlo simulation of three-dimensional Finite Element Models. Our goals were to assess the feasibility and accuracy of our method by determining whether it can (i) construct an accurate model of displacements in healthy tissues from a reasonable number of principal components and samples in the data-set; (ii) detect, locate and indicate the extent of an injury and (iii) better characterize injuries than strain-based and modulus-based methods.

2 Method description

Suppose $\mathrm{\Gamma }\subset R^m$ is a region of interest in an ultrasound image, where a point $\mathit{x}\in \mathrm{\Gamma }$ denotes a specific location. Let1 $$\begin{array}{c}\hfill \mathit{u}:\mathrm{\Gamma }\subset {\mathbb{R}}^m\to {\mathbb{R}}^n,\end{array}$$1 be a map of displacements obtained from quasi-static UE, where $m,n$ are positive integers satisfying 1 $\le n\le m\le 3$. That is to say, $\mathit{u}(\mathit{x})$, the image of $\mathit{x}$ under $\mathit{u}$, represents displacements measured at the point $\mathit{x}$. The restriction $1\le n\le m\le 3$ implies that the images are at most three-dimensional and displacements are estimated only in directions that are captured in the image. Finally, we assume that the displacement map $\mathit{u}$ is defined on a finite set of points $\{{\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_p\}$, with $p\in \mathbb{N}$, and was measured using a specific protocol on a specific soft tissue structure. For example, $\mathit{u}$ could be measured by imaging the long-axis of the Achilles tendon while the ankle is extended. It follows that displacements described by $\mathit{u}$ can be represented by an $np$-dimensional vector. Also, we can consider the set of all possible maps that use the specific protocol to examine specific structure when healthy. We designate the latter set as $V$. Although maps in $V$ are generated using a specific protocol and specific structure, these maps may be generated with variable loading conditions and imaging windows and for structures with variable anatomy and material properties.

The relationship between $\mathit{u}$ and $V$ can be exploited to obtain information about the underlying tissue of interest (Figure  1). Specifically, injuries can be identified by simply verifying whether $\mathit{u}$ is not in $V$. Otherwise, it is possible that the displacements could have been estimated in healthy tissues. Since injured tissue can be detected whenever $\mathit{u}$ is not in $V$, we conjecture that the difference between $\mathit{u}$ and $V$ can also help locate the injury and evaluate its extent. To quantify this difference, we introduce the following definition.

Definition 2.1

Given a displacement map $\mathit{u}\in {\mathbb{R}}^{np}$ and the set of displacement maps from healthy tissues $V$, damage displacement $\mathit{w}$ is defined as2 $$\begin{array}{c}\hfill \mathit{w}:=\mathit{u}-arg\underset{\mathit{v}\in \overline{V}}{min}\left|\left|\mathit{u},-,\mathit{v}\right|\right|,\end{array}$$2 where $\overline{V}$ denotes the closure of $V$ and $\left|\left|·\right|\right|$ is the standard Euclidean norm.

We remark on two aspects of this definition. The first is that the name damage displacement is chosen to imply that $w$ does not represent actual displacements, even though, like the original displacements $u$, it has units of displacements and lies in ${\mathbb{R}}^{np}$. The second is that the definition does not presuppose that damage displacement is unique. Additional conditions on $\overline{V}$, such as convexity, are necessary to ensure uniqueness. In what follows, we address this issue by approximating $\overline{V}$ with a set $\widehat{V}$ so that $arg{min}_{\mathit{v}\in \widehat{V}}\left|\left|\mathit{u},-,\mathit{v}\right|\right|$ contains exactly one element.

Figure 1. Damage displacement, $w$, can be represented pictorially as the difference between $u$, a map of estimated displacements and $V$, the set of displacement maps that can be estimated in healthy tissue.

Principal component analysis (PCA) is used to approximate the set $\overline{V}$ defined above. We must approximate $V$ due to the impossibility of measuring all possible displacements maps that may occur in healthy tissues. For example, there are an infinite number of possible imaging planes that could be used during the examination. We use a finite set of vectors $\}{\mathit{u}}_i{\{}_{i=1}^K\subset V$ for ${\mathit{u}}_i\in {\mathbb{R}}^{np}$ and $K\in \mathbb{N}$ to build principal components.[26] Each vector ${\mathit{u}}_i$ represents an event in which the specified protocol is applied to the specified structure that is healthy. Each entry in the vector represents a measured displacement at a point in the image during this event. Let ${\mathit{p}}_i\in {\mathbb{R}}^{np}$ denote the principal component with the $i$th largest principal value ${\mathit{\lambda }}_i$. The first $J$ principal components can explain a fraction of the variation in the set $\}{\mathit{u}}_i{\{}_{i=1}^K$, denoted by ${var}_J$:3 $$\begin{array}{c}\hfill {var}_J:=\left(\sum _{j=1}^J,{\mathit{\lambda }}_j\right){\left(\sum _{k=1}^K,{\mathit{\lambda }}_k\right)}^{-1}.\end{array}$$3 The set $\overline{V}$ is approximated with a subset of the principal components that explains a specific fraction of the variation.

Definition 2.2

Given a threshold for the fraction of explained variation $\mathit{\zeta }\in [0,1]$ and a set of displacements from healthy tissues, $\}{\mathit{u}}_i{\{}_{i=1}^K$. The set $\overline{V}$ can be approximated as4 $$\begin{array}{c}\hfill \widehat{V}:=\left\{\overline{\mathit{u}},+,\sum _{i=1}^L,a_i,{\mathit{p}}_i,:,a_i,\in ,\mathbb{R},,,i,=,1,,,\dots ,,,L\right\},\end{array}$$4 where $\overline{\mathit{u}}:=K^{-1}\left({\sum }_i^K,{\mathit{u}}_i\right)$, ${\mathit{p}}_i$ is a principal component of $\}{\mathit{u}}_i{\{}_{i=1}^K$ with the $i$th largest corresponding eigenvalue, and $L$ is the smallest integer such that ${var}_L>=\mathit{\zeta }$.

We note that each $a_i\in \mathbb{R}$ is often restricted to a bounded interval in other applications of PCA. [27] Additionally, $\widehat{V}$ is non-empty, closed, and convex as long as $K\ge 1$ so that the projection theorem [28, Chapter XII] ensures that damage displacement exists and is unique. More importantly, we obtain a closed-form expression of $\mathit{w}$ given orthonormal principal components ${\mathit{p}}_1,\dots ,{\mathit{p}}_L$:5 $$\begin{array}{c}\hfill \mathit{w}\approx \mathit{u}-arg\underset{\mathit{v}\in \widehat{V}}{min}\left|\left|\mathit{u},-,\mathit{v}\right|\right|=(\mathit{u}-\overline{\mathit{u}})-\sum _{i=1}^L〈{\mathit{p}}_i,\mathit{u}-\overline{\mathit{u}}〉{\mathit{p}}_i,\end{array}$$5 where $〈·,·〉$ is the standard Euclidean dot product.

Lastly, we must specify where to measure displacements in the region of interest $\mathrm{\Gamma }$. When using PCA, it is implicitly assumed that the $k$th displacement in ${\mathit{u}}_i$ corresponds with the $k$th displacement in ${\mathit{u}}_j$ ($i\ne j$, $i,j=1,\dots ,K$ and $k=1,\dots ,np$). We build a correspondence between displacements by building a correspondence between the points at which the two displacements are measured. One possibility is to consider two points at the same location with respect to the imaging coordinates, i.e. axial and lateral coordinates. However, one point could lie inside the structure in one image and outside in another, since imaging windows and anatomy are subject to change among patients and examinations. A better choice is to build a correspondence between points in an image with the technique known as image registration. That way, a point in one image can be located relatively within the same region of the tissue as a point in another image. Lets suppose that the specified protocol captures the same region of the specified structure. We then have the following definition.

Definition 2.3

Suppose $U$ is some region in ${\mathbb{R}}^m$ and $\}{\mathit{q}}_i{\{}_{i=1}^p$ is a set of points in $U$. Let $\mathit{m}$ define a map that registers $U$ to a region of interest in the image $\mathrm{\Gamma }$ satisfying $\mathit{m}(U)\subset \mathrm{\Gamma }$. Then, the $k$th displacement is measured at the point $\mathit{m}({\mathit{q}}_k)$ for $k=1,\dots ,p$.

Of course, how we define $U$ and $M$ will depend on the image registration technique, for which there are numerous options (see [29] for a list of options). In the next section, we use affine mappings to register sets.

3 Numerical study

We used a numerical study to assess the feasibility and accuracy of our method. Our assessment was comprised of three specific aims: determining whether our method can (i) construct an accurate model of displacements in healthy tissues from a reasonable number of principal components and samples in the data-set (Section 3.2), (ii) detect, locate and indicate the extent of an injury (Section 3.3) and (iii) better locate injuries than strain-based and modulus-based methods (Section 3.4). We accomplished the three aims by performing Monte Carlo simulation, where each simulation event represented an examination of musculoskeletal soft tissue with quasi-static UE. Monte Carlo involves randomly generated parameters, which we used to capture the variation between examinations and patients. We begin this section by constructing a model of each simulation event and describing the random parameters associated with the event (Section  3.1).

3.1 Event model

We restricted the events (examinations) to the case when a specific quasi-static UE protocol is applied to a specific soft tissue structure. The specific structure was modelled after a bone–ligament–bone complex, such as the Patellar ligament, and will be referred to simply as a ligament. The specific protocol consisted of stretching the ligament along its long-axis while imaging the ligament with a linear transducer and then estimating axial displacements in the imaging window.

We calculated axial displacements using Finite Element Analysis to estimate three-dimensional displacements associated with stretching of the ligament. Our Finite Element Analysis is detailed in Appendix 1 and implemented with Matlab code, which can be found at https://sites.google.com/site/amylouisecochran/code. Next, we calculated axial displacements relative to the imaging window associated with each event and then used an affine mapping to register the axial displacements in the imaging window to the square $U:=[0,1]\times [0,1]$. This procedure generated a map of axial displacements $\mathbf{u}$ for each event.

To represent such complex events, we used 17 random parameters to specify properties of the tissue (anatomy and material properties) and the examination (loading conditions and imaging windows). For injured tissues, we used an additional four random parameters to specify properties of the injury (location and extent of the damage). All parameters were uniformly and independently distributed (Table 1). Although random variables are often normally distributed in Monte Carlo simulations, we used uniformly distributed random variables to avoid sampling undesirable values, such as negative material properties. When it is necessary to refer to a location, the variables are defined with respect to the Cartesian coordinates $(x,y,z)\in {\mathbb{R}}^3$ given in Figure 2.

Table 1. Ranges of uniform random variables that were used for simulation.

Symbol	Description	Range
$l$	Length	$3.13$–$5.26$ cm^a
$s_h$	Width : length	$0.64$–$0.71$^a
$s_w$	Height : length	$0.19$–$0.21$^a
$r_1,r_2$	Bone radii	$3.6$–$4.4$ cm
$h_1,h_2$	Asymmetry constant	$-0.2$–$0.2$ cm
${\mathit{\theta }}_1,{\mathit{\theta }}_2,{\mathit{\theta }}_3$	Rotation	$-2^{\circ }$–$2^{\circ }$
$t_x,t_y,t_z$	Translation	$-0.2$–$0.2$ cm
$c_1$	Matrix stiffness	$7.5$–$8.2$ MPa^b
$c_3$	Fibre stiffness	$0.24$–$0.26$ MPa^b
$c_4$	Exponent constant	$57.4$–$63.4$^b
$\mathit{\gamma }$	Stretch	$0.02$–$0.04$
$(x_0,y_0,z_0)$	Damage centre	Within the ligament
$\mathit{\alpha }$	Damage level	$0.25$–$0.75$

a Based on dimensions of patellar ligaments.[30]

b Based on material constants of medial collateral ligaments.[31]

Anatomy resembled a bone–ligament–bone complex and was represented by seven random variables: $l,s_w,s_h,r_1,r_2,h_1,h_2$ (Figure 2). Variables $l$, $s_w$ and $s_h$ designated the ligament’s length, width-to-length ratio, and height-to-length ratio, respectively. Bony insertions were modelled as cylinders with radii $r_1$ and $r_2$. We used different radii for each cylinder so that the anatomy did not have an unrealistic symmetry about the $yz$-plane. The cylinders were translated away from the $x$-axis so that the anatomy was also asymmetric about the $xz$ plane. One insertion was translated in the $y$-direction by $h_1$; the other was translated in the $y$-direction by $h_2$. The ranges for length are exactly those found for patellar ligaments; average value for width : length are based on the ratio of the average length and average width for patellar ligaments; and average value for width : length are based on twice ratio of the average length and average width for patellar ligaments.[30] (The latter ratio was multiplied by two to avoiding thin structures and any resulting numerical stability.)

Figure 2. Cross-sections of a bone–ligament–bone complex in the three Cartesian planes. Seven random variables represented ligaments of various sizes, proportions and shapes.

Imaging windows were represented by six random parameters: ${\mathit{\theta }}_1,{\mathit{\theta }}_2,{\mathit{\theta }}_3,t_x,t_y,t_z$. We use this representation in order to capture the six degrees of freedom that are inherent to any plane. The length of the simulated transducer was aligned almost parallel to the long-axis of the ligament. To be specific, the imaging window was initialized as the $xy$-plane. The normal to the imaging window was then rotated, first around the $z$-axis by ${\mathit{\theta }}_1$ degrees, second around the $y$-axis by ${\mathit{\theta }}_2$ degrees, and finally around the $z$-axis by ${\mathit{\theta }}_3$ degrees. The imaging window was then translated by $t_x$, $t_y$ and $t_z$, along the $x$, $y$ and $z$ directions, respectively.

Material behaviour was modelled as hyperelastic, transversely isotropic (fibrous), and nearly incompressible.[32] Suppose$$J:=\sqrt{det(\mathbf{C})},{\overline{I}}_{\mathsf{1}}:=\mathrm{tr}(\overline{\mathbf{C}}),\overline{\mathit{\lambda }}:={\left({\mathbf{m}}_{\mathsf{0}},·,\overline{\mathbf{C}},·,{\mathbf{m}}_{\mathsf{0}}\right)}^{1/2},$$which are invariants of the right Cauchy-Green strain tensor $C$ and its isochoric portion given by $\overline{\mathbf{C}}=J^{-2/3}\mathbf{C}$. The vector ${\mathbf{m}}_0$ defines the fibre direction and $\mathrm{d}et$ and $\mathrm{t}r$ denote determinant and trace, respectively. Material behaviour was captured by a strain energy function split into isochoric (first two terms) and volumetric portions (last term):6 $$\begin{array}{c}\hfill c_1({\overline{I}}_1-3)+W_f(\overline{\mathit{\lambda }})+\mathit{\kappa }(J-1-lnJ).\end{array}$$6 The first term represents an isotropic contribution, known as the contribution from the matrix, whereas the second isochoric term, given by the function $W_f$, captures the contribution from the fibres. We use the piecewise strain energy function to define $W_f$ [32]:7 $$\begin{array}{c}\hfill \overline{\mathit{\lambda }}{W}_f(\overline{\mathit{\lambda }})=\left\{\begin{array}{cc}0\hfill & \overline{\mathit{\lambda }}<1\hfill \\ c_3\left(exp,\left[c_4,\left(\overline{\mathit{\lambda }},-,1\right)\right],-,1\right)\hfill & 1\le \overline{\mathit{\lambda }}<{\mathit{\lambda }}^{\ast }\hfill \\ c_5\left(\overline{\mathit{\lambda }},+,c_6\right)\hfill & {\mathit{\lambda }}^{\ast }\le \overline{\mathit{\lambda }}<\mathit{\lambda }\hfill \end{array}\right.\end{array}$$7 Material behaviour was represented simply with only three random parameters: $c_1,c_3,c_4$. The parameters $\mathit{\kappa }$ and ${\mathit{\lambda }}^{\ast }$ were fixed at $3,750$ MPa and $1.05$, respectively, and the parameter $c_5$ was chosen such that $c_3:c_5$ was fixed at $0.25:307.5$. Parameters $\mathit{\kappa }$ and ${\mathit{\lambda }}^{\ast }$ and the ratio $c_3:c_5$ were chosen to agree with those found in [31]. The remaining material parameter $c_6$ was chosen so that $W_f(\mathit{\lambda })$ is continuous. Additionally, fibre direction $m_0$ was aligned with the long-axis of the ligament, given by the $x$-axis in the figure. Loading conditions were modelled using one random parameter, $\mathit{\gamma }$, which defined a uniaxial stretch along the structure’s long-axis.

Examinations of injured tissue were modelled using the same random variables that represented loading, imaging, material properties and anatomy for healthy tissue. Four additional random variables: $\mathit{\alpha },x_0,y_0,z_0$, represented a region of weakened mechanical integrity given by a sphere of radius $0.5$ cm and centre $(x_0,y_0,z_0)$. The constants: $c_1$, $c_3$ and $c_5$ were decreased in the damaged region by a fraction $\mathit{\alpha }$. Soft tissue damage has previously been modelled by decreasing material constants by a fraction.[33] We note that injured regions and imaging windows were randomly located in the tissue. As a result, the injured regions did not always overlap with the imaging window.

3.2 Accuracy and feasibility of PCA Model

Using the event model described above, we simulated 80 events (examinations) involving healthy tissue and calculated a map of axial displacements $\mathit{u}$ for each event. These displacements were then used to construct 11 different PCA models. The first four models were constructed using $80$ samples, but with a varying number of principal components: two, four, six and eight. The next four models were constructed using a fixed number of samples, but with a varying number of samples: $20$, $40$, $60$ and $80$. The fixed number of principal components was determined based on accuracy of the first four PCA models, and the fixed number of samples was determined based on accuracy of the next four PCA models.

The remaining three models were constructed for a fixed number of events and a fixed number of principal components, but with a varying level of displacement error, which we modelled as Gaussian white noise (i.e. normally distributed with mean zero and standard $\mathit{\sigma }$, to each displacement estimate). The standard deviation $\mathit{\sigma }$ was defined for each map $u$ as $\mathit{\zeta }\overline{\mathit{u}}$, where $\overline{\mathit{u}}$ is the average displacement magnitude of displacements over all $80$ events and $\mathit{\zeta }$ was varied between $0$, $0.02$, and $0.04$. These PCA models are used for the remaining tests.

We simulated an additional four groups of $50$ events involving healthy tissue to test accuracy of the PCA models. Each group was assigned one of the first four PCA models and one of the next four PCA models. With only three noise levels, only three of the four groups were assigned one of the last three PCA models. For each event and each of its assigned PCA model, we measured axial displacement $\mathit{u}$, damage displacement $\mathit{w}$ and the relative magnitude of damage displacement given by $||\mathit{w}||/||\mathit{u}||$. The relative magnitude was a measure of relative error in the PCA model, since theoretically $\mathit{w}$ should be zero for healthy tissue. Whenever the assigned PCA model incorporated noise, we added noise at the same level to axial displacements $\mathit{u}$.

3.3 Method testing

We tested the ability of our method to detect, indicate the extent of and locate an injury. First, we designated a ligament as injured whenever the norm of damage displacement was larger than a threshold value. Sensitivity and specificity were measured for a range of threshold values, and a receiver operating characteristic (ROC) curve was generated.

Second, we plotted the norm of damage displacement against the damage level ($\mathit{\alpha }$), and tested for positive correlation between the two variables. We used $150$ of the same events involving healthy tissue and simulated an additional $150$ events involving injured tissue. For the new events, we measured axial displacement and added noise at the three levels; we then measured damage displacement and the relative magnitude of damage displacement.

Last, we measured contrast between injured and healthy regions from images of damage strain. Damage strain was defined as the partial derivative of damage displacement with respect to the axial coordinates. Damage displacement was linearly interpolated onto a grid with $80\mathrm{\mu }\mathrm{m}$ spacing in the axial and lateral directions. Interpolated displacements were used to calculate damage strain via a linear least squares approach [34] with a 15-by-15 sample kernel to calculate damage strain, which was interpolated linearly back onto the initial points. Contrast between injured and healthy regions was measured using the contrast-to-noise ratio (CNR) [35]:8 $$\begin{array}{c}\hfill \mathrm{CNR}:=\frac{2{({\mathit{\tau }}_{in}-{\mathit{\tau }}_{out})}^2}{{\mathit{\sigma }}_{in}^2+{\mathit{\sigma }}_{out}^2},\end{array}$$8 where $\mathit{\tau }$ and $\mathit{\sigma }$ denote mean and standard deviation of damage strain, respectively, measured both inside and outside the injured region denoted by subscripts $in$ and $out$.

Contrast-to-noise ratio (CNR) is commonly used to compare how well techniques can identify damaged regions and is believed to be well correlated with visual impression.[36] For phantom and in vivo studies, regions inside and outside the pathological region are manually defined by the viewer and used to calculate CNR (see e.g. [37]). The entire pathological region may be difficult to delineate, but a region clearly inside the pathological region and clearly outside the pathological region can frequently be determined and used to calculate CNR. For simulation studies, the damage region is actually known precisely a priori. Hence, this information can be explicitly used to calculate CNR and then used to compare methods. Indeed, this is exactly how we calculated CNR. Because CNR is well defined only when the damaged region is captured in the image, we simulated additional events involving injured tissue until we had three groups of $50$ events for which CNR was well defined and could then be compared for different noise levels.

3.4 Method comparison

Using the same test described above, we compared the ability of our method to locate an injury against that of the strain-based and modulus-based methods. We did not compare the ability of our method to detect and indicate the extent of an injury, because there is not a universally accepted scalar measurement of strain or modulus that is taken entirely within the tissue of interest and that is used for diagnosing injuries. To compare methods, we used the same events from the previous sections, and whenever independence was needed for a statistical comparison, we divided events into groups of $50$. The strain-based method measured axial strain $a$, the partial derivative of axial displacement with respect to the axial coordinates. Axial strain was calculated in the same manner as damage strain except with damage displacement replaced by axial displacement. The modulus-based method measured shear modulus $\mathit{\mu }$, which was calculated using the adjoint method [17, Appendix 2].

3.5 Statistics

We used the statistical software JMP (SAS Institute, Cary, NC) to perform statistical tests, with the exception of Games–Howell post-hoc tests. We used SPSS software (SPSS, Chicago, IL) to perform the latter test, because JMP does not have Games–Howell post-hoc tests. We determined the appropriate statistical test for differences between means by checking for normality using a Shapiro–Wilks test and for equal variances using four available tests: O’Brian, Brown–Forsythe, Levene and Barlett. In the case of non-normality, we identified significant differences between means with Kruskal–Wallis’s test, and in the case of unequal variances, we identified significant differences between means with Welch’s test. Additionally, multiple comparisons were performed with a Dunn’s test in the case of non-normality and with Games–Howell’s test in the case of unequal variances. Positive correlation was tested using the linear fit tool available in JMP. Statistical significance was considered $p<0.05$.

4 Numerical results

4.1 Accuracy and feasibility of PCA Model

Through statistical comparisons, we found that the PCA model could achieve sufficient accuracy with six principal components (Figure 3). To arrive at this conclusion, we first found that the relative magnitude of damage displacement had significantly unequal variances between the PCA models with two, four, six and eight principal components ($p<0.0001)$. Welch’s test found significant differences between means among the four PCA models ($p<0.0001$), and Games–Howell tests found statistical differences between each pair of model ($p<0.001$), except between the PCA model with six principal components the PCA model with eight principal components ($p=0.787$). Because there were no statistical differences in accuracy when using six principal components rather than eight, even though there were statistical differences in accuracy when using six principal components rather than two or four, we considered six principal components to be sufficient, and hence, we fixed the number of principal components at six in subsequent models.

With the number of principal components fixed at six, we found that the PCA model could achieve sufficient accuracy with $40$ samples (Figure  3). We found that the relative magnitude of damage displacement satisfied the hypothesis of equal variances among the PCA models with $20$, $40$, $60$ and $80$ samples. ($p>0.29$), but rejected the hypothesis that the relative magnitude of damage displacement was normal ($p=0.0047$). A Kruskal–Wallis test revealed significant differences between means for these four PCA models ($p=0.0047$) and multiple comparisons found that the significant differences arose between the PCA model with $20$ samples and the PCA model with $60$ samples ($p=0.0135$), as well as between the PCA model with $20$ samples and the PCA model with $80$ samples. The PCA models with $40$ samples did not give rise to any significant difference when compared to the other PCA models. Therefore, we concluded that $40$ samples was sufficient and, subsequently, fixed for the remaining models.

When testing the effect of noise level on the PCA model, we found statistical differences between the PCA models with different noise levels: $0$, $2$ and $4\%$ ($p<0.0001$; Figure 3). Statistical significances were tested using a Welch’s test due to rejecting the hypothesis of equal variances ($p<0.03$). Games–Howell test revealed statistical differences between each pair of the three PCA models. The relative magnitude of damage displacement dropped from $2.7$ to $3.4\%$ when noise level of $2\%$ was added to the displacements and then dropped to $4.6\%$ when noise level reached $4\%$.

Figure 3. Relative magnitude of damage displacement of PCA models with different (a) number of principal components, (b) number of samples and (c) noise levels. Relative magnitude was a measure of the relative error in the PCA models. There was no noise in models used in (a) and (b). The number of principal components were fixed at six for models used in (b) and (c). The number of samples were fixed at 80 for models used in (a) and at 40 for models used in (c).

4.2 Method testing

Magnitude of damage displacement was able to detect injured tissues, but its accuracy decreased with each increase in noise level (area under ROC curve of $0.86$ for $0\%$ noise, $0.85$ for $2\%$ noise and $0.81$ for $4\%$ noise; Figure 4). Injured tissue was most accurately detected when the damaged region was imaged ($N=81$; area under ROC curve of $0.93$ for $0\%$ noise, $0.92$ for $2\%$ noise and $0.90$ for $4\%$ noise), but was still detected when the damaged region was not imaged ($N=69$; area under ROC curve of $0.77$ for $0\%$ noise, $0.76$ for $2\%$ noise and $0.71$ for $4\%$ noise).

Figure 4. Receiver operator characteristic (ROC) curves generated from the magnitude of damage displacement for (a) different noise levels and (b) different types of injury in the absence of noise: all injuries, damaged region is captured ($N=81$) or damaged region is not imaged ($N=69$). Detection accuracy improved with decreasing noise level and when the damaged region is imaged.

Relative magnitude of damage displacement was positively correlated with the damage level at each noise level ($N=150$ and $p<0.0001$ for each noise level) (Figure  5). Despite positive correlation, an increase in damage produced only a relatively small increase in magnitude, as indicated by the relatively flat regression lines (slopes $=0.102$ for $0\%$ noise, $0.097$ for $2\%$ noise and $0.088$ for $4\%$ noise). Additionally, the regression line explained only a small fraction of the data variation (adjusted $R^2$ equals $0.21$ for $0\%$ noise, $0.20$ for $2\%$ noise and $0.19$ for $4\%$ noise).

Figure 5. Damage level versus relative magnitude of damage displacement at different noise levels: (a) no noise, (b) $2\%$ noise and (c) $4\%$ noise. Damage level was positively correlated with relative magnitude of damage displacement at each noise level ($N=150$ and $p<0.0001$ for each noise level).

Damaged regions were more easily located in maps of damage strain than maps of damage displacement (Figure  6). We found that CNR for maps of damage strain did not satisfy the hypothesis of normality when there was no noise ($p=0.0023)$. Tests for unequal variances were inconclusive ($p=0.016-0.1040$ among the four available tests) so that we tested for differences in means using both Welch’s test and Kruskal–Wallis test. Both tests concluded that CNR values were not significantly different for the different noise levels ($p=0.27$ for Welch and $p=0.35$ for Kruskal–Wallis; Figure 7). Each group had $50$ samples.

Figure 6. Maps of damage displacement and maps of damage strain. damage displacement was calculated from axial displacements with (a) no noise, (b) $2\%$ noise and (c) $4\%$ noise; damage strain was then calculated from damage displacement in each case: (d) no noise, (e) $2\%$ noise and (f) $4\%$ noise. Damaged regions are marked by the dotted line, and data-to-colour transformations were defined based on the mean and standard deviation of values in each map to remove bias towards maps with larger values.

Figure 7. Contrast-to-noise ratio (CNR) with no noise, $2\%$ noise, and $4\%$ noise. Mean values of CNR did not differ significantly between each noise level ($p\ge 0.27$). Each group had $50$ samples.

4.3 Method comparison

Our method was better able to locate an injury than the strain-based and modulus-based methods. Specifically, maps of damage strain had significantly better contrast than maps of axial strain and maps of damage strain (Figure  8). When there was no noise, maps of damage strain had an average CNR of $7.1$ compared to $3.2$ for maps of axial strain and $0.25$ for maps of shear modulus. CNR did not satisfy the hypothesis of normality for maps of damage strain ($p=0.0023$ for no noise, $p=0.0036$ for $2\%$ noise and $p=0.0084$ for $4\%$ noise) and had significantly unequal variances between all three types of maps ($p<0.0001$ for each noise level). Welch’s test found significant differences between means of CNR among the three methods ($p<0.0001$), and Games–Howell test found that significant differences between each pair of methods ($p<0.001$ at each noise level). Damaged regions were difficult to delineate in maps of shear modulus and experienced more negative values in maps of axial strain (Figure  9).

Figure 8. Contrast-to-noise ratios for our method, strain-based method, and modulus-based method at each noise level: (a) $0\%$, (b) $2\%$ and (c) $4\%$. Our method produced significantly better contrast-to-noise ratios than the other two methods ($p<0.001$ at each noise level).

Figure 9. Maps of axial strain and shear modulus. Axial strain was calculated from axial displacements with (a) no noise, (b) $2\%$ noise and (c) $4\%$ noise; shear modulus was calculated from axial displacements with (d) no noise, (e) $2\%$ noise and (f) $4\%$ noise. Damaged regions are marked by the dotted line. Damaged regions are marked by the dotted line, and data-to-colour transformations were defined based on the mean and standard deviation of values in each map to remove bias towards maps with larger values.

5 Discussion

We presented a new method for detecting injuries from internal measurements of displacements taken while a quasi-static load is applied to the tissue. We wanted to overcome issues with two standard approaches used with quasi-static UE: the strain-based method and the modulus-based method. Chief among these issues was the lack of quantitative methods for detecting injuries, likely caused by variation in displacements among patients and examinations. Rather than controlling this variation, we aimed to capture this variation in a PCA model of displacements acquired from healthy tissues. We then characterized injured tissue for how its displacements differed from the PCA model; this difference was captured in a measurement we called damage displacement. We tested feasibility and accuracy of our method with a simulation of quasi-static UE examinations. Our goals were to determine numerically whether our method can (i) construct an accurate model of displacements in healthy tissues from a reasonable number of principal components and samples in the data-set; (ii) detect, locate and indicate the extent of an injury and (iii) better locate injuries than strain-based and modulus-based methods.

We found that we could construct an accurate model of displacements in healthy tissues with $40$ samples and $6$ principal components. Adding either samples or principal components saw no significant improvement to accuracy ($p\ge 0.47$), and so, we fixed the number of samples at $40$ and the number of principal components at $6$ for later tests of our method. Between the number of samples and the number of principal components, we anticipate that the former is a more pressing concern. The number of principal components is associated with the number of vector inner products and additions, a rather small amount of computation when compared to calculating strain and modulus. We consequently believe that six principal components would not greatly limit the feasibility of our method. In contrast, how many samples are needed from healthy volunteers may highly determine whether our method is feasible. Although more than $40$ samples are likely needed in practice, needing only $40$ samples here leaves open the possibility that the samples could be reasonably obtained in practice. In contrast to the number of samples and principal components, noise level degraded accuracy significantly with each increase in this parameter. However, even at the highest noise level ($4\%$), our PCA model had an average error of $4.6\%$.

We also found that our method could accurately detect, indicate the extent of, and locate injuries. Best detection accuracy was achieved when there was no noise and when the damaged region was imaged; the area under the ROC curve was high as $0.93$ in this case. Even when there was noise and when the damaged region was not imaged, our method could detect injuries with the area under of the ROC curve dropping to $0.71$. Additionally, we found that the damage level was positively-correlated with the relative magnitude of damage displacement at each noise level ($p<0.0001$), but the associated regression line was relatively flat ($\approx 0.1$) and explained only a small fraction of the variation (adjusted $R^2=0.19-0.21$). These last two results may prevent our method from accurately predicting the extent of the injury. Lastly, an injury was easier to locate in maps of damage strain than damage displacement, and the contrast in the former maps did not significantly degrade with increasing noise levels.

Our method was also better able to locate injuries than the strain-based and modulus-based methods, as indicated by significantly higher contrast-to-noise ratios between the healthy and damaged regions ($p<0.001$). With quasi-static UE, strain-based and modulus-based methods are currently used to detect injuries from estimated displacements, with injuries detected primarily through a visual inspection of maps of strain and modulus.[5–7, 9, 10, 10, 12–14, 38, 39] Visual inspection has led to an diagnostic accuracy of $94-97\%$.[6, 14] One approach to visual inspection is to examine strain “patterns”, a broad term presented in [7] which refers to how strain appears to be spatially distributed. Their concept was that strain patterns of healthy volunteers differed from strain patterns of injured tissues; that way, injured tissues could be identified for these differences. In this regard, their approach is similar to ours, since we assume that displacement maps of healthy volunteers differ from displacement maps of injured tissue and that these differences can help identify injured tissues. However, our method explicitly defines how a variable is spatially distributed and also detects differences quantitatively rather than qualitatively. Even ignoring its ability to detect injuries quantitatively, our method provided maps in our numerical study with significantly better contrast than maps of the strain-based and modulus-based methods. This result alone suggests that our method may potentially facilitate the visual inspection of injuries.

Whereas visual inspection is commonly used to detect injuries, quantitative approaches are less common for applications of quasi-static UE. Therefore, a potential benefit of our method is that it provides a scalar measurement taken within the tissue of interest (relative magnitude of damage displacement), which can be used to detect and indicate the extent of the injury quantitatively. In contrast, strain-based and modulus-based methods do not provide a reliable scalar measurement taken within the tissue of interest. For reference, we were able to confirm this limitation by testing the magnitudes of axial strain and shear modulus in a similar way we tested relative magnitude of damage displacement. When there was no noise, the result was an ROC curve of $0.56$ for the strain-based method and $0.51$ for the modulus-based and a lack of positive correlation between either magnitude and damage level ($p>0.16$).

Two studies have suggested a scalar measurement based on strain outside the tissue of interest, in addition to strain within the tissue. The first study evaluated normal Achilles tendons using the strain index, i.e. the ratio of strain within the tendon to strain in surrounding tissue,[7] whereas the second evaluated muscles using the ratio of strain within the muscle to strain within an reference gel, located between the skin and transducer.[40] We did not simulate surrounding tissue, and so, we were unable to directly compare our scalar measurement against either of these two measurements. However, our method may help in situations where these two measurements are limited. Both of these two scalar measurements are ideal when the transducer applies uniform pressure; certain structures, such as the supraspinatus, cannot be examined by applying transducer pressure.[41] Additionally, the strain index often requires adipose tissue located outside the tissue of interest and may not reflect relative stiffness in two regions if stress is unequal in the two regions. For example, Drakonaki et al. [7] reported average strain ratios greater than one when comparing the Achilles tendon to adipose tissue beneath the tendon. They also reported relatively high variability in the strain index between examiners and between examinations (coefficient of variation $=30\%-39\%$).

Our numerical study evaluated strain-based and modulus-based methods when the unhealthy tissue lay outside the imaging window. This setting highlights how detecting an injury may be incongruent with estimating modulus. Namely, when the injury lay outside the imaging window, only healthy tissue is imaged, and so, the modulus would reflect healthy tissue. If the modulus is estimated correctly, then the injury would go undetected. Only when the modulus is estimated imperfectly may the injury be diagnosed. Clearly, it is better to detect an injury rather than estimate modulus in the imaging window. In general, tissue displacements are affected by both local and non-local changes to the modulus. Our method focuses on differences in displacements between healthy and injured tissue rather than the modulus. The outcome was an ability to detect injuries when the damaged region lay outside the imaging window.

This paper considers the problem of detecting injuries from internal measurements of displacements taken while a quasi-static load is applied to the tissue. However, injuries can be evaluated from internal measurements of displacements using other elasticity imaging techniques–Magnetic Resonance Elastography,[42–46] shear wave ultrasound imaging,[47–50] transient elastography,[51] among others.[52] These alternative imaging techniques use dynamic loading of the tissue rather than quasi-static loading and promise quantitative diagnoses of injuries, as evidenced in [42]. However, like quasi-static UE, there are also concerns about reproducibility and variability in elasticity among healthy volunteers.[43, 49] Recent mathematical and numerical frameworks: small volume expansions,[53] multistatic response-based imaging,[54–56] topological derivatives,[57] and time-reversal algorithms [58] may help to locate pathological regions in soft tissue. In particular, multistatic response-based imaging and topological derivatives have performed well in the presence of limited apertures and noisy measurements.[54–57] It would be interesting to extend these frameworks to the case of quasi-static loading.

Along with the methods described above, there are other methods that are similar to our proposed method. That is, they either (i) measure how displacements are affected when a pathological region is added to an otherwise normal region [53–57] or (ii) use a data-set acquired beforehand into elasticity imaging or tracking of tissue motion.[59–64] The mathematical techniques: small volume expansions, multi-static imaging, and topological derivatives use the governing equations of tissue motion to determine how adding a pathological region affects displacement measurements. The method based on small volume expansions, in particular, calculates the difference between measured displacements and displacements in the absence of any anomaly.[53] Regarding the use of a data-set, a data-set of normal and malignant tissue was used to develop a classifier for data obtained from their in-house designed system, the Breast Mechanical Imager.[59] Also, a data-set, as well as principal component analysis, was used to model pathological motion of ventricle walls in the heart.[61, 62] In contrast to these applications that use a data-set, we used a data-set acquired only from healthy tissues. This lends to certain advantages: training samples are easier to obtain, data variation is smaller without injured tissue in the data-set, and training does not require initial classification of the data-set.

We conclude this paper by discussing three primary limitations of our study. First, noise levels of displacements were lower than levels found experimentally. Our study used low noise levels ($2\%$ and $4\%$) comparable to the $1\%$ noise level used to test the adjoint method.[17] However, internal measurements of displacements are expected to have noise levels much higher. For example, Walker and Trahey [65] remark that a lower bound of $2.09$ ns in the standard deviation of measured displacements would accompany an expected measured displacements of $13$ ns. This leads to a noise level of 16%. Our study saw that increasing noise levels had a measurable effect on the ability to diagnose injuries quantitatively and an insignificant effect on CNR. Hence, we expect that increasing the noise level to $16\%$ will also degrade diagnostic ability, while affecting CNR to a lesser degree.

A second limitation is that we did not present any experimental results. We expect there will be greater difficulty to acquire the data-set and register the images experimentally. Additionally, there will be additional sources of variation in the experiment, such as irregular material properties, boundaries and correlated noise. Consequently, a large training set may be needed to capture the displacement variation among healthy tissues. To address this variation, we may need to model displacements in healthy tissues using a statistical technique other than PCA, such as robust PCA.[66] Robust PCA is more robust in the presence of data error. Additionally, we could use the model of displacements in healthy tissue to define initial search spaces during displacement estimation.[27, 36, 67, 68] The final limitation is that we did not compare our method against other modulus-based methods, which may perform better than the modulus-based method we tested by including either non-linearity [21, 69] or anisotropy.[20]

In summary, we provided the first step in validating a new method for detecting injuries of soft tissue. We incorporated principal component analysis into quasi-static UE resulting in better ability to contrast than strain-based and modulus-based methods during a numerical study. Future study is needed to test our method experimentally and in the presence of realistic noise levels.

Appendix 1

Axial displacements

Tissue displacements were calculated for each event using three-dimensional finite element analysis (FEA). This analysis was performed in Matlab (Mathworks, Natick, MA) using two publicly available toolboxes: distMesh [70] and iFEM.[71, 72] We used distMesh to construct tetrahedral meshes for the anatomy associated with each instance and used iFEM to perform FEA using sparse matrices and adaptive mesh refinement. This last toolbox was modified to handle non-linear material behaviour and vector-valued FEA.

Maps $u\in {\mathbb{R}}^p$ were constructed from axial displacements estimated from the FEA. We let $U:=[0,1]\times [0,1]$, $p:=101\times 501=50,601$, and $\}{rmq}_i{\{}_{i=1}^p$ be given by:$${\mathrm{q}}_i:=\left[\genfrac{}{}{0pt}{}{q_x^i}{q_y^i}\right]:=\left[\genfrac{}{}{0pt}{}{{\displaystyle b+\overline{ı}\frac{1-2b}{500}}}{{\displaystyle b+\left(i,-,151,\overline{ı},-,1\right)\frac{1-2b}{100}}}\right],$$for $i=1,\dots ,50601$ and $\overline{ı}:=\lfloor (i-1)/101\rfloor$ (Definition  2.3).

The set $\}{\mathit{q}}_i{\{}_{i=1}^p$ defines a grid in $U$. The constant $b$ was set to $0.01$ to ensure measurements within the ligament. The ratio of the grid’s spacing in the $x$ direction to the grid spacing in the $y$ direction was set to the mean value of $s_w(=0.2)$ so that each point after registration represented nearly a square region in the imaging plane. Furthermore, point spacing after registration was approximately $5\mathrm{cm}/500\mathrm{points}=100\mathrm{\mu }\mathrm{m}/\mathrm{point}$ in the imaging coordinates, which is within the range of axial resolution in ultrasound applications.

For each event, an affine mapping $\mathit{m}$ was defined in three steps (please refer to Figure A1). In the imaging window, we first considered the line that halved the structure across its height and measured the distance, call it $D$, from one bone to another along this line. Next, we marked the two points on the line located a distance $0.025D$ from each bone and drew perpendiculars to this line at these two points. The perpendiculars intersected the boundary of the structure at four new points, denoted by their lateral and axial coordinates by $l_i$ and $a_i$ , respectively, with $i=1,\dots ,4$. The four points were enumerated clockwise with respect to Figure A1, starting from the bottom left corner. The corners of $U$ were registered to the four points in the image via the mapping:

0$$M(u_x,u_y):=\left[\genfrac{}{}{0pt}{}{l_4-l_1{a}_4-a_1}{l_2-l_1{a}_2-a_1}\right]\left[\genfrac{}{}{0pt}{}{u_x}{u_y}\right]+\left[\genfrac{}{}{0pt}{}{l_1}{a_1}\right]\mathrm{for}\left[\genfrac{}{}{0pt}{}{u_x}{u_y}\right]\in U$$

Figure 10. Affine mapping $\mathit{m}$ from the unit square to the image of the ligament. On the right side, the solid line delineates the ligament and the dashed line delineates the range of $\mathit{m}$.

Appendix 2

Modulus-based method

Following the adjoint method in [17], we defined shear modulus, $\mathrm{\mu }$, in terms of a dummy variable $\mathit{\xi }$: $\mathrm{\mu }:=1+{\mathit{\xi }}^2$. Suppose $\mathrm{\Omega }$ is the convex hull of points that were defined for $u$, $u_m$ is the function given when $u$ is linearly interpolated over $\mathrm{\Omega }$, $u(\mathit{\xi })$ is axial displacements when calculated from $\mathit{\xi }$, and ${\left|\left|·\right|\right|}_{L^2(\mathrm{\Omega })}$ denote the $L^2$-Norm over $\mathrm{\Omega }$. Then, we optimized the following function with respect to $\mathit{\xi }$:B.1 $$\begin{array}{c}\hfill {\left|\left|u,(\mathit{\xi }),-,u_m\right|\right|}_{L^2(\mathrm{\Omega })}+\mathit{\alpha }{\left|\left|1,+,{\mathit{\xi }}^2\right|\right|}_{L^2(\mathrm{\Omega })},\end{array}$$B.1 where $\mathit{\alpha }$ is numerical parameter.

Displacements $u(\mathit{\xi })$ were estimated using two-dimensional Finite Element Analysis. The region $\mathrm{\Omega }$ was meshed with linear, triangular elements, and the tissue was assumed to undergo plane strain and have material behaviour that is linear, isotropic and incompressible. Tractions were set to zero at the boundary parallel to the structure’s long-axis, and displacements were set to actual axial and lateral displacements on the rest of the boundary. Equation (B1) was optimized with the adjoint method and interior-point method with limited memory BFGS.[17] Optimization was performed in Matlab (Mathworks, Natick, MA). The adjoint method involves two numerical parameters, $\mathit{\alpha }$ and $\mathit{\beta }$, which we set to ${10}^{-8}$ and ${10}^3$, respectively.