Model validation for a noninvasive arterial stenosis detection problem

A current thrust in medical research is the development of a non-invasive method for detection, localization, and characterization of an arterial stenosis (a blockage or partial blockage in an artery). A method has been proposed to detect shear waves in the chest cavity which have been generated by disturbances in the blood flow resulting from a stenosis. In order to develop this methodology further, we use both one-dimensional pressure and shear wave experimental data from novel acoustic phantoms to validate corresponding viscoelastic mathematical models, which were developed in a concept paper [8] and refined herein. We estimate model parameters which give a good fit (in a sense to be precisely defined) to the experimental data, and use asymptotic error theory to provide confidence intervals for parameter estimates. Finally, since a robust error model is necessary for accurate parameter estimates and confidence analysis, we include a comparison of absolute and relative models for measurement error. Mathematics Subject Classification: 62F12; 62F40; 65M32; 74D05.


Introduction
Coronary artery disease (CAD) is an increasingly prevalent medical condition, and is often a precursor to and cause of a patient experiencing cardiac arrest. Current methods for detection of arterial stenoses (blocked arteries) include the angiogram and CT scans. Angiograms are viable but quite invasive, while CT scans are expensive, introduce radiation into the patient, and can only detect hard plaques (blockages). A desirable new detection method would be noninvasive and less expensive but still effective. To this end, using acoustic waves generated by stenoses has been proposed as a possible detection method [1,3,21,24,25]. This would place sensors on the surface of the chest to listen for sounds from coronary arteries, with the hope of detecting and then localizing any blockages.
In keeping with [3], we note that the entire system couples two processes: (1) the generation of pressure and shear waves transmitted into the body by the arterial wall as a result of the turbulent blood flow generated by a stenosis, and (2) the propagation of pressure and shear waves through the chest to sensors attached to the chest wall. The first process is not completely understood, though some basic ideas are present in the literature. The current understanding [3,21] is that turbulent flow produces normal forces on the vessel walls at and downstream of a stenosis, which then exert pressure on the vessels wall causing a small displacement in the surrounding soft tissue. Previous work (e.g. [1,21,22,25]) has demonstrated the existence of such sounds, and [24] discusses past work done on building devices to detect sounds from coronary arteries. As has been noted in previous work (see, e.g., [3,21]), the focus in practice is on detecting shear waves, due in part to the fact that faster pressure waves have a wavelength that is too long for the distance scales in the body and the corresponding observations that the speed and frequency of shear waves generated by a stenosis are all sufficiently low enough to make shear sound detection plausible.
In terms of the second process, propagation through the chest cavity, it is also known that sounds from coronary stenoses have been sufficiently strong to be detected; although rare, [24] notes that sounds from coronary stenoses have been heard with a stethoscope. The difficulty in detecting these sounds is that they are weaker than others such as those caused by the heart valves and larger blood vessels. The sounds are also attenuated and scattered during their passage through the body to the chest wall, which makes them difficult to detect and localize without a deeper understanding of the physics underlying the wave propagation through the chest (hence the need for mathematical models). The approach taken by the work summarized in [24] is treating the problem as a signal-to-noise problem where the signal is weak relative to the noise. In a line of work by Banks, et al., [3,13,14,15,16,20,23], the sound propagation problem was approached from a mathematical modeling standpoint. Though signal processing would likely be an important component of a diagnostic device, the approach taken in the Banks, et al., line of work was to build models to describe the underlying physics of sound propagation. These models allow for a characterization of coronary stenoses, which will assist in uncovering these particular coronary artery sounds from the noisy background in the body. Initial experiments were conducted where a gel mold was built with a tube running through the middle; cases where the tube was unblocked were compared to those with partial blockages, and the results suggested that there were significant differences in sound generation between the blocked and unblocked cases. Unfortunately, this line of work ended before experimental data could be incorporated and fitted to models. The current work picks up the Banks, et al., ideas, starting again with a one dimensional model and experimental setup.
In this document, we continue the work of our concept paper [8] by focusing on wave propagation through a viscoelastic medium. Here we develop a slightly more general constitutive relationship than in our concept paper, use this constitutive relationship in one-dimensional pressure and shear wave dynamical models, and demonstrate successful inverse problem results for the one-dimensional case using experimental data from a tissue-mimicking gel mold. This data comes from novel acoustic phantoms built and tested at Queen Mary, University of London (QMUL) and Barts Health Trust (BHT) in England. We will examine both the pressure and shear cases to provide evidence of the fidelity of our model fit to data.

Experimental Setup
We begin with a discussion of the experimental setup. Two separate novel experiments have been devised at QMUL to gather pressure and shear data; though we will discuss them together, the experiments are run at completely different times and with slightly different phantoms. Devices have been designed (see left panes of Figures 1-2), in which an agar gel mold phantom (homogeneous, 97% water, density ρ = 1010 kg/m 3 ) is loaded into the rig, a weight is attached applying stress to the phantom, and then the weight is released, causing the material to oscillate. The displacement motion of the material throughout the experiment is measured with a laser device. The choice of loading and a quick release is designed to produce dynamic data; the idea was inspired in part by the impacts the stenosed vessel wall experiences with each heartbeat and also by past success in gathering shear data for filled rubber elastomers using an initially loaded rubber sample which then underwent an impulsive hammer hit (see e.g. [13,14]). This yields one dimensional pressure data along the vertical axis in the pressure case (right pane of Figure 1) and, in the second experiment, shear data in the radial direction perpendicular to the vertical axis (right pane of Figure 2).
In this work, we focus on results generated from a load weight of 264 g, as this weight level produces a well-defined response. In the future, we may also incorporate data from smaller weight levels in order to examine the different phantom responses to other weight levels. The gel phantoms were stored in water when not in use, which keeps the gel at the desired 97% water composition.
When the experiment is conducted, data like those depicted in Figure 3 are produced. The material is at rest, a weight is added and allowed to settle, then the string holding the weight is rapidly cut with a flame to allow the material to freely oscillate. Once oscillations have died out the material relaxes back toward a stable state.
The key pieces that will be modeled are the loading profile (loading begins at t = Γ 1 and lasts until the weight begins to be released at t = Γ 3 ), which we will model as instantaneous loading to position A, and the oscillations after weight release (free oscillations begin at t = Γ 4 = 0) which are the main object of investigations here.
With the setup of the experiment in mind, we can turn to our mathematical model of wave propagation. The model will be developed to take into account all features of the data, including the loading profile and the relaxation present in data. For more information on the experimental setup, interested readers may refer to [18].  Figure 3: Sample one-dimensional data. Loading of the material (initially at rest) begins at t = Γ 1 , and the material is loaded and continues to relax for t ∈ (Γ 2 , Γ 3 ). At time Γ 3 the load is cut which takes roughly 10-15 ms. The gel is then freely oscillating at Γ 4 = 0, and oscillations continue for a period of time dependent on the loading weight and wave type (pressure/shear). The value A is the displacement of the material at the beginning of free oscillations. The overall displacement scale of the data is on the order of 10 −4 m, while the oscillations immediately after the weight release are on the order of 10 −5 m.

Model Development and Constitutive Equation
Since our phantom is cylindrical, the model development begins with three-dimensional equations of motion in cylindrical coordinates. These are given in [23, p.20], and also in [20], and are derived from momentum and mass balance principles. Using the fact that the gel is homogeneous in both the pressure and shear cases and that there are symmetries in the experimental design, these three-dimensional equations can be reduced to simplified onedimensional models for both cases. In the pressure displacement case, the governing partial differential equation where ρ is the density of the material, the stress tensor σ is given by the constitutive relationship for the material (the form of which will be discussed later), g(t) is a function that describes the loading process (again, to be discussed later), and the material is initially at rest. The value u(x, t) represents the displacement of the material at position x and time t, with x ∈ (0, L) and t > Γ 1 . The time Γ 1 is chosen as the beginning of any stress-strain history in the material; we are assuming the material has been at rest long enough that it is only affected by displacements for t > Γ 1 , where Γ 1 is the time when we start modeling the material history. For our device, L = 0.0518 m is the height of the phantom.
In the shear displacement case, the governing equation becomes where ρ, σ, and g(t) are analogous to the pressure case and where r ∈ (r min , r max ) for t > Γ 1 . For our device, r min = 0.0105 m and r max = 0.054 m. Throughout this work, we will use r as the spatial variable when the model is for shear displacement and x as the spatial variable for pressure displacement.
In order to complete these models, we must provide a form for σ. This is the constitutive relationship, also called the stress-strain law since it relates strain ( ∂u ∂x or ∂u ∂r ) and/or the strain rate to stress σ. The next sections discuss this aspect of the model.

Constitutive equation
We incorporate the previous modeling ideas together into a new constitutive equation for the pressure (1) and shear (2) wave PDEs. Throughout this section, the constitutive relationship form is the same for the pressure and shear cases, so x and r are interchangeable unless otherwise noted; for notational convenience, we use x as the spatial variable in the discussion which follows. We develop the constitutive equation which will be used in this work assuming that we will be solving the model starting at t = Γ 1 and thus incorporating both the loading process and oscillations into our dynamic equations, for the time being. In Section 3.1.5, we will make an approximation to the loading process which will allow us to focus on the dynamic oscillations of the material after the weight is released, which is our true interest.

Fung model
Some of the initial investigation into the viscoelastic nature of tissue was completed by Fung (see [7] and the references therein). His work is of particular interest because it was validated in actual tissue. Fung developed a "quasi-linear" model with a kernel of the form Within (3), λ represents the stretch of a material (λ = 1 + u x ) and σ e describes the elastic response to the elongation λ, given by (see [3]) σ e (λ) = −β + βe αux where α and β are constants to be estimated (and where we use u r in the shear case). The parameters τ 1 and τ 2 are lower and upper bounds, respectively, on relaxation times, which describe the ways in which the material responds to imposed stresses and strains. This model incorporates a continuum τ ∈ [τ 1 , τ 2 ] of relaxation times, which Fung found to be necessary in order for his model to match the response of tissue, as well as a constant term in the kernel. This Fung kernel will serve as a baseline which we will refer back to when developing the model for this paper.

Linearized constitutive equation
One could keep nonlinearities in the constitutive equation (3). However, we found (as we shall see later) that a linear constitutive relationship gives a reasonable approximation to the data provided by QMUL and BHT. To that end, we will use the first two terms of the Taylor expansion of e αux to approximate where we have combined γ = βα into a single parameter to be estimated; γ will be incorporated into other parameters later in model development. We can then linearize (3) by using (5), add a Kelvin-Voigt damping term (a common linear viscoelastic damping model [7]), and obtain where G(t) is a kernel to be specified. To an extent, the Kelvin-Voigt term describes the overall nature of the damping present in the material, while the kernel G(t) will incorporate different material responses at both the macroscopic and microscopic levels.

Existence and uniqueness for pressure and shear models
Before moving on to the specific form of the constitutive equation kernel, we first establish existence and uniqueness for the pressure (1) and shear (2) equations with the constitutive equation (6). To that end, we set up a similar framework as in the concept paper [8] and connect those results to the current work. We will require that the following assumptions hold: (A1) The boundary condition function satisfies g ∈ L 2 (Γ 1 , T ); (A2) The kernel G is differentiable with respect to t ∈ R + and with constants c 1 and c 2 such that |G(t)| ≤ c 1 and |Ġ(t)| ≤ c 2 for all t ∈ R + .

Pressure case
The pressure PDE (1) with constitutive equation (6) are of the same form as those in Section 2 of [8], except that here we have the initial time denoted as t = Γ 1 instead of t = 0 and slightly different variable names (inconsequential changes). Let H = L 2 (0, L), V = {φ|φ ∈ H 1 (0, L), φ(0) = 0}, and V * denote the topological dual space of V. We identify H with its topological dual H * and thus obtain V ֒→ H = H * ֒→ V * as a Gelfand triple [2,26]. The notation ·, · denotes the inner product in H, and ·, · V * ,V represents the duality pairing between V * and V. Let Then a weak solution u ∈ L T for the pressure equation must satisfy for any t ∈ [Γ 1 , T ] and η ∈ L T . Here and elsewhere u(t) and η(t) denote the functions u(·, t) and η(·, t), respectively. With these definitions, we still have that the following theorem (a restatement of Theorem 2.2 in [8]) holds:

Shear case
This requires a bit more consideration. The shear domain is Ω = [r min , r max ], and is solved on the time frame t ∈ [Γ 1 , T ]. We must slightly redefine the spaces from above to fit the shear model. Let H = L 2 (r min , r max ), V = {φ|φ ∈ H 1 (r min , r max ), φ(r max ) = 0}, and V * denote the topological dual space of V. We identify H with its topological dual H * and thus again obtain V ֒→ H = H * ֒→ V * as a Gelfand triple. Let C w (Γ 1 , T ; V) denote the set of weakly continuous functions in V on [Γ 1 , T ], and Then a weak solution u ∈ L T for the shear equation must satisfy for any t ∈ [Γ 1 , T ] and η ∈ L T and where ·, · is the usual inner product. Since r min > 0, there are no singularities in the final term in (8), and the kernel integral in the numerator of that term will converge in the same manner as the preceding kernel integral. Thus, the arguments from [8] for the pressure case apply in the shear case, and we have the following theorem:

Form for constitutive equation kernel G(t)
We will now state the particular kernel used for this current work, and then manipulate it into a form that gives more physical insight and which will later allow for a conceptual framework using internal variables. We develop this kernel from a different perspective than that given in [8], but the resulting form will be quite similar. Using the notation and parameter conventions of [7], we define the kernel in this work to be where κ r is a constant representing an instantaneous relaxation modulus (justified by the fact that our gel phantom acts partly as a solid) and K(t; P ) = T exp(−t/τ )dP (τ ) represents a continuum of distributed relaxation times with T = [τ 1 , τ 2 ] ⊂ (0, ∞) and where P (τ ) is a probability measure on T . Note that this form for G satisfies |G(t)| ≤ c 1 with G clearly differentiable and |Ġ(t)| ≤ c 2 for some constants c 1 , c 2 so that assumption (A2) is satisfied. It is also worth noting here that our proposed kernel form (9) is similar to that in Fung's model (4), as we see that κ r serves as an analog to the constant portion of Fung's kernel (i.e., 1 1 + c ln(τ 1 /τ 1 ) ) and the K(t; P ) portion is similar to the the continuous relaxation spectrum in Fung's model (i.e., c τ2 τ1 ).
We substitute (9) into (6) and manipulate the form of the stress with integration by parts, noting that u x (Γ 1 ) = 0 since the material is initially at rest and using the fact that that K(0; P ) = 1: where E = κ r γ and with slightly more detail in [9]. This equation (10) is the general form of the constitutive equation used here. The value E 0 = E + γ can be considered to be a dynamic analog to the static Young's modulus in the pressure case or the static shear modulus in the shear case; this also makes clear the fact that Hooke's Law is incorporated into our model. We have already discussed that E 1 is the bulk damping parameter for the Kelvin-Voigt damping term. The final integral represents a history term which describes the relaxation of the material in response to an applied stress/strain. We will ultimately turn to a discretized distribution model (using a discrete measure P (τ )), and connect it to the continuum model through a probability measure approximation as in [4]. This will allow us to develop a computationally feasible inverse problem, and also give insight into the underlying material mechanics. But first we briefly discuss a method for approximating the loading process.

Approximating the loading process
Recalling Figure 3, the loading profile is relatively long compared with the oscillatory period; since our concern is with modeling the oscillations, solving the model from Γ 1 is much longer than necessary. Also, early experimentation with the model indicated that the parameters governing the loading and resting process may differ from those governing the very dynamic post-release oscillatory process.
We address these concerns by modeling the loading as instantaneous from at rest to a displacement of A at position x = L or r = r min . Since the material is linear, this would then mean the phantom has the profile u(x, t) = A L x in the pressure case and u(r, t) = A(r max − r) r max − r min in the shear case, up until the time of the weight release. Since this is an approximation, we will neglect the times t ∈ (Γ 3 , Γ 4 ), the weight release time period, since that time frame is small relative to the loading and settling time from Γ 1 to Γ 3 . We also incorporate a time parameter Υ which will represent our approximation of the time when loading begins. In the formulation here we will use the same relaxation times during the loading process as during the oscillation period, which means that Υ has no meaning other than as a tuning parameter that we must estimate. Thus, we assume the given loading profiles for t ∈ (Υ, 0) since Γ 4 = 0 in our convention. This also means that Υ < 0. We incorporate this loading approximation into our model by manually integrating the constitutive relationship (10). For the purposes here, we will callσ the full constitutive relationship for t > Υ that is described by (10) (where we now use Υ in the place of Γ 1 ), and σ the constitutive relationship for t > 0. We do this for notational simplicity in the final model, at the expense of some minor notational confusion at the current stage.
For the pressure case, we compute (noting that u( incorporates the remaining terms and represents the constitutive relationship for t > 0. For the pressure setup, we then have the following: • The original stress boundary condition isσ(L, t; P ) = 0. Using the preceding development, this corresponds with 0 = σ(L, t; P ) − F (t; Υ, A, P ) which allows us to write the boundary condition for a model solved for t > 0 as σ(L, t; P ) = F (t; Υ, A, P ).
Since A, Υ < 0 and γ > 0, we know that K(t; P ) − K(t − Υ; P ) > 0. Hence, since L > 0, we have a compressive boundary stress, which is what we would expect in the pressure case.
The shear case is similar. For the loading profile u(r, t) = A(rmax−r) rmax−rmin for t ∈ (Υ, 0), we have u r (r, t) = − A rmax−rmin which is incorporated when integrating the history in the same way as the pressure case. We then find the corresponding loading stress to be F (t; Υ, A, P ) = −ζ A rmax−rmin (K(t; P )−K(t−Υ; P ), where ζ is the shear analog to γ. Also, we haveσ r = σ r as in the pressure case. However, we have the termσ r = σ r − F (t; Υ, A, P ) r , which will result in a time-dependent forcing term in the shear PDE.
We make two comments before discussing the internal variable forms. First, if we assume, for example, a single relaxation time and that its value is small, say on the order of 10 −1 , then the term K(t − Υ; P ) = exp(−(t − Υ)/τ 1 ) ≈ exp(−10(t − Υ)) attains its maximum value exp(10Υ) when t = 0. Note that for, say, Υ < −1, this term is negligible. Relaxation times on this order are what we can later obtain in the inverse problem, which would imply that in our case the material is at rest after being loaded sufficiently long that it "forgets" its loading history by the time the weight is released. This is good from an experimental standpoint, since the loading process will never be quite uniform. It is also good to know from a computational perspective; we can limit Υ to being greater than some value, such as −20 < Υ < 0, which will keep the optimization algorithm from marching off unnecessarily (which occurred in some of our early inverse problem tests). Second, since we have integrated out the loading history, we now start the model at the time of weight release, t = 0. This means that the material is considered at rest just prior to the release; thus, in the history integrals we will discuss in the next section, all the history now starts at t = 0 since the history before that point will be incorporated into the initial loading profile and an initial stress condition.

Internal variable formulation
In the previously noted work on this stenosis problem, the double integrals that resulted from using the continuum of relaxation times in the stress equation were computationally intractable so another approach was required. The idea was to use a discrete number of internal variables. As will be noted, these gave rise to a differential form which was an improvement computationally since it led to purely differential equations in the model rather than inclusion of an integro-differential equation. With the advances in desktop computation abilities since that time, the integral form is now reasonable to use in a dynamic model. However, internal variables are still attractive in that they provide a formulation that indicates some of the internal material dynamics. Physically, if we assume that the molecules within the biological tissue are on a microscopic scale then the portion of the material which is represented by each internal variable or internal strain ǫ j is being driven by the overall strain and has a response that varies depending on the value of the corresponding relaxation time τ j .
Previous work [3,20,23] assumed a discrete internal variable form as an approximation to the Fung kernel, using the nonlinear constitutive equation (3). This discrete form assumed the kernel G(t) was in an exponential form, with the effects brought together as a discrete sum in the constitutive equation. The results in [3,20,23] demonstrate that the internal variable approach is valid and does appear to work as well as the continuum of times in the Fung kernel. A connection between the Fung continuum model and the discrete kernel is provided by the work in [15]. The authors there form the kernel G(t) = T q(t; τ )dP (τ ) where T = [τ 1 , τ 2 ] ⊂ (0, ∞) is the set of admissible relaxation times, P (τ ) is a probability measure on T , and q(t; τ ) is a continuous function of relaxation times. If we take q(t; τ ) = exp(−t/τ ), this corresponds with the kernels previously discussed. The authors showed existence and uniqueness results for this kernel in the nonlinear constitutive equation (3). Then, a result from [4] allows one to approximate any measure P (τ ) with a discrete measure. This discrete measure approximation leads us back to the case with a sum of exponentials, but from the probabilistic framework we know conclusively that we are approximating the continuous spectrum of Fung's kernel and that this approximation has been viable when implemented.
With this background on previous work using internal variables in hand, we move forward by modifying our current model. We manipulate the form of Equation (11) as follows: where we let Rather than the integral form for ǫ 1 , we can use the differential form which is then solved simultaneously with the rest of the model dynamics. This is then an internal variable or internal strain, driven by the overall strain u x (t), which is the continuous form of the internal variable formulation. We now may finally make the discrete assumption where ∆ τj is the Heaviside function with step at τ j and p j are the proportions of the material subject to relaxation time τ j . By substituting this discrete P into the form for σ as developed in (12), we obtain the discrete, internal variable form of the constitutive relationship with internal variables obeying (for j = 1, 2, . . . , N p ) and where we have defined γ j = γp j so that γ = Np j=1 γ j . Note that we assume E > 0, since the agar gel acts at least partly as a solid, and that ǫ j = ǫ 1 (·; t j ).

Final PDE pressure and shear models
We now put together the pressure (1) and shear (2) wave equations, using the constitutive equation (14) but with the loading history approximation incorporated as discussed in Section 3.1.5. Recall also that the discrete assumption for P and the form of K gives us γK(t; P ) = γ Np j=1 p j exp(−t/τ j ) = Np j=1 γ j exp(−t/τ j ) where γ j = γp j . These equations are just manipulated versions of the general equations of Theorems 3.1-3.2, so we still know a unique weak solution exists on any finite time interval.

Pressure Model
The pressure equations, solved for t > 0 which is the release time, are then where with the internal variables subject to (for j = 1, 2, . . . , N p ) The parameter ρ still represents the material density and E 1 the damping parameter. E represents an instantaneous relaxation modulus. The γ j values are weightings on relaxation times τ j ; also, we can write E 0 = E + Np j=1 γ j as the viscoelastic analog to Young's modulus.

Shear Model
We next present the shear equations. In order to more easily distinguish between pressure and shear model parameters, we will use G and G 1 in place of E and E 1 and ζ j instead of γ j . We find where with the internal variables subject to (for j = 1, 2, . . . , N p ) We note that G 0 = G + Np j=1 ζ j is the dynamic analog of the shear modulus.

Numerical method
We use the same numerical implementation for both the pressure (15) and shear (16) models, which is also the same as in the concept paper [8]. In time, we use a discontinuous Galerkin method composed of normalized Legendre polynomials (of order 4). In space, we use a continuous spectral finite element method composed of Lagrange basis functions on Gauss-Lobatto nodes (also of order 4). This allows the higher order (4-5th order) elements in space while controlling dispersion error. Under this scheme, the system matrices are diagonalizablethis could be lost if we did not use normalized Legendre polynomials in time. Also, the solver time for the shear model is slower than that for the pressure model due to the time-dependent right hand side. Further details are in a forthcoming BICOM report [19].

Inverse Problem
With models in hand, we now turn to matching the model output to data. We will use two common methods in order to estimate model parameters. One is ordinary least squares (OLS) and the other is generalized least squares (GLS). These will be defined later in Section 4.1.
As discussed in Section 2, separate experiments have been designed to gather one-dimensional pressure and shear data. Measurements in our experiment are taken at x = L for the pressure case and r = r min for the shear case, and will be denoted u j . Corresponding pressure or shear model solutions at the same spatial location will be denoted u(t j ; 10 θ ), where the measurement location has been suppressed so we can retain a general pressure/shear model solution notation and where θ represents a vector of the base-10 logarithm of each parameter (the same idea used previously [8] to reduce parameter scaling issues). Each data set has been trimmed to the dynamic oscillations after the release, and thus the time frame for pressure data is roughly 150 ms while that for shear is 200 ms. The data were sampled at a rate of 2048 Hz; however, this high rate proved to make the inverse problem difficult and computationally intractable because that many data points resulted in the inverse problem being over determined. Instead, we will use every other data point from the larger data set for a sampling rate of 1024 Hz which will we later refer to as the "every data point" set. We take n to be the total number of data points for a particular data set, and thus can describe the measurement time points for the full "every data point" set as t j = j/1024 where j = 0, 1, . . . , n − 1. There will also be a reduced data set where we take every other data point starting with t 0 = 0, which corresponds with a data sampling rate of 512Hz. For pressure data we will report results for both data sets for contrast, while for shear the results were sufficiently similar between data sets that we will report only the results corresponding with the 512 Hz sampling rate.
Since some of the data points were near zero in absolute value, we found that those points resulted in scaling problems when using the GLS model to estimate model parameters (since the corresponding cost functional divides by the model value as we will see later when this method is defined). To account for this, we removed from consideration any data points u j (and their corresponding model solutions at that time point) where |u j | < 5 × 10 −6 . This value was chosen by examining the data, noting that the data is on the order of 10 −5 and that the "jitter" one can see in Figure 3 has a magnitude of roughly 5 × 10 −6 during the times before loading up to Γ 1 , then during the settling period from Γ 2 to Γ 3 , and again in the settling period after the oscillations have died out. Thus, our threshold level is below the level of noise in the data. This level also eliminated only a few data points, while providing significantly improved GLS robustness. The number of data points n is then reduced according to how many thresholded data points were removed.
Before going into the setup and results for the inverse problem, we note that the forward (i.e., direct) problems where we solve for displacement (using the method discussed in Section 3.3) are as follows: • Pressure forward problem: Given E, E 1 , τ j and γ j for j = 1, 2, . . . , N p , Υ, A, L, and ρ, solve model (15) for displacement u(x, t) at each position x ∈ [0, L] for t ∈ [0, T ].
The inverse problems we will develop here are as follows: • Pressure inverse problem: Given pressure displacement data at x = L and a corresponding forward problem solver for displacement, along with specified values for ρ and L, find values for the constants E, E 1 , τ j and γ j (for j = 1, 2, . . . , N p ), A, and Υ which provide the best fit to the data (in a manner which will be defined shortly).
• Shear inverse problem: Given shear displacement data at r = r min and a corresponding forward problem solver for displacement, along with specified values for ρ, r min , and r max , find values for the constants G, G 1 , τ j and ζ j (for j = 1, 2, . . . , N p ), A, and Υ which provide the best fit to the data (again, in a manner which will be defined shortly).
We assume for both the pressure and shear cases that the parameters lie in some admissible set Q ⊂ R κ , where Q is assumed to be compact and κ is the number of parameters requiring estimation. Throughout the remainder of this work, we will denote the log-scaled parameter vector for pressure (for N p = 1 and κ = 6) as θ = (log 10 (E), log 10 (E 1 ), log 10 (γ 1 ), log 10 (τ 1 ), log 10 (−A), log 10 (−Υ)) (17) with a similar vector for shear where we use the shear parameters parameters G, G 1 , and ζ 1 in place of E, E 1 , and γ 1 , respectively. Thus, as long as we define our cost function to be a continuous function of the parameters, we know the inverse problem has a solution (minimizing a continuous function on a compact parameter space). One could broaden this parameter estimation framework to the distributional case if desired, taking an admissible parameter space as a compact subset of Euclidean space (including all parameters excuding relaxation times) along with with the space of probability measures, and use the Prohorov metric framework (see, e.g., [10,Sec. 4]) and the approximation results of [4]. This again leads to minimizing a continuous function of the parameters over a compact space. Either way, the inverse problems we will shortly define will have solutions.

Statistical Models and Parameter Estimators
In order to carefully define the way in which we will measure the closeness of the data to model values, we must first discuss underlying statistical models for the error present in the data. A proper error model is also key to correctly determining parameter confidence intervals. Much of the discussion here is similar to that in [8], with background on ordinary least squares (OLS) and generalized/weighted least squares (GLS or WLS) given in [17], for example. We will assume the errors E j are independent, identically distributed with mean zero (E[E j ] = 0) and constant variance var(E j ) = σ 2 0 ; this process has realizations ε j . Note that we do not assume we know the underlying distributions from which the errors come; we only know the first two central moments as specified. We use this error process in proposing two error models and corresponding parameter estimators.
• Absolute error: Here we have the error process U j = u(t j ; 10 θ0 ) + E j , with realizations where θ 0 is some hypothesized "true" parameter value (see [17]). We use the ordinary least squares cost function The corresponding inverse problem for the logged parameters is then This function minimizes the distance between the data and model where all observations are considered to have equal importance (weight). Since u(t j ; 10 θ ) is a continuous function of θ, J ols is also a continuous function of θ, which means we are minimizing a continuous function of θ over a compact set Q, and thus this inverse problem has a solution.
For this case, we construct the generalized (weighted) least squares cost function (as per, e.g., [17]) where we define the weights w j = u(t j ; 10 θ ) −1 . In this case, since we are examining a relative error model (20), these weights take into account the unequal quality of observations; dividing by the function value has a "normalizing" effect on the errors, accounting for the scale differences which may be present in the errors at larger versus smaller model values.
We now wish to find θ such that J gls (θ) is minimized. We can either solve this directly, or by using an iterative procedure in order to estimateθ gls (since the weights must also be estimated). We will use an iterative method, described as follows (see [17] and references therein for convergence details): 1. Defineθ 0 =θ ols , and set k = 0.
4. Set k = k + 1 and return to Step 2. Terminate when successive estimates forθ gls are sufficiently close, or when one has iterated 20 times. For our problem, our "sufficiently close" criterion was found by determining if ||θ k+1 −θ k || ∞ ≤ 10 −3 , where ||θ|| ∞ is the maximum component of the given vector θ. The parameter values being estimated are all log-scaled, and are thus on the order of [10 −1 , 10 1 ]. This puts the stopping criterion at two orders of magnitude less than the smallest log-scaled parameter value, which is sufficient in our problem.
Even though we have removed all data points with absolute value under 5 × 10 −6 , we still account for the (now unlikely) possibility that some model values may still end up small in absolute value. Thus, we incorporate thresholding on the weights to keep from dividing by zero. We take a weight threshold value of 1 × 10 −10 , as this is almost certainly below the threshold of significance in terms of the model displacements.
With weight thresholding, we are assured that the iterative process is possible numerically. Thus, similar to the ordinary least squares case, at each step k in the iterative GLS estimation process we are minimizing a continuous function of θ over a compact parameter space Q, and thus the inverse problem in each iteration will have a solution. Also, as long as the iterative process is carried out sufficiently many times, under certain conditions the weights will convergeŵ j → u(t j ; 10θ gls ) −1 (see, e.g., [17]).

Optimization considerations
As in [8], we used the Matlab routine lsqnonlin for our optimization routine to solve forθ ols andθ gls . We used the trust-region-reflective (TRR) algorithm that is built in; as our previous effort in [8] demonstrated, the Levenburg-Marquardt option was slower than TRR and did not give us better results. Since we are using at least one relaxation time, we do not consider fmincon which we have shown to be ineffective in estimating relaxation times. In order to start the optimization routines for computingθ ols , we must provide initial parameter values (for θ gls we use the estimated value forθ ols as our initial guess). From a perusal of the viscoelastic materials literature, our experience with the previous conceptual work, and from some manual examination on the current data sets, we developed pressure and shear initial values.

Residuals
We will also include residual plots to assist in analysis of the model fit to data, and to indicate which error model best describes the error in the data. Residuals give a sense for the model fit to data, but more importantly the residuals can give an indication [17] regarding the appropriateness of our error model. If the absolute residuals seem to be randomly dispersed around the horizontal axis and form a horizontal band around that axis, then the absolute error model may be correct. On the other hand, if the (modified) relative residuals seem to be randomly dispersed, then the relative error model may be correct. We define the following: • Absolute residuals are computed as r j = u j − u(t j ; 10θ), whereθ is the particular parameter estimate being considered.
• Relative residuals are computed as r j =ŵ j (u j − u(t j ; 10θ)) whereŵ j = u(t j ; 10θ) −1 and theŵ j are thresholded in the same manner as discussed earlier.

Inverse problem results, N p = 1
We now demonstrate the ability of our model to match data. For this purpose, we will take a single relaxation time (N p = 1) as that is enough to show model fidelity to data. We run both the absolute (OLS) and relative (GLS) error models on a sample data set using a 264 g loading weight, separately for both pressure and shear data. We will report parameter estimates, standard errors and confidence intervals, plots of model fits to data, plots of residuals versus time, and plots of residuals versus model values. We use these elements in order to recommend error models for the pressure and shear cases. Standard errors (and corresponding confidence intervals) are computed using asymptotic error theory. For the absolute error model, the process is the same as that which we used in [8], and is also described in [17,Ch. 3]; for the relative error model, the corresponding asymptotic error methodology is discussed in [17,Ch. 3]. Since the theory is common enough, we do not reiterate it here and refer interested readers to the aforementioned references.
We have also examined parameter estimation using data sampled at different rates. This allowed for a study of whether the parameter estimates and associated confidence intervals remain consistent as the number of data points is reduced. Using fewer data points is also a way of decreasing computational times for the inverse problems. We ran the inverse problem on each data set and using each error model using all the data points (1024 Hz) and using every other data point (512 Hz), as discussed at the beginning of Section 4. We report results from both sampling rates for the pressure case, while only presenting the 512 Hz results for the shear case since the results at each sampling rate were similar.
Before showing results, we must note that the parameter Υ is a special case. As discussed in Section 3.1.5, if the relaxation times are small (which they will be in our results) then Υ will not have much of an effect on the model output once it becomes sufficiently negative. Thus, even though the model output is at least somewhat sensitive to Υ [9], particularly for the pressure model, we may still obtain large confidence intervals for this parameter once sufficiently negative since further increases in the parameter will not effect model output. We will see these larger confidence intervals for Υ in the results, and it should be noted that this is not a major concern since it is an artificial parameter designed to approximate the loading process.

Pressure data results
The pressure parameter estimates and confidence intervals (see [5,6,8,11,17] for information on computing confidence intervals), obtained from the inverse problems using OLS and GLS on using both data sampling frequencies, are shown in Tables 1-4. Model fits as well as residual plots are shown in  In all cases, model fits to data are good. Comparing the OLS results in Tables 1 and 3 with their GLS counterparts in Tables figures, the residuals versus model plots are also not much different when comparing the OLS cases to the GLS cases. From a residual analysis standpoint, then, either model appears reasonable. Thus, since OLS is a simpler error model and since we have a higher degree of confidence in the parameter estimates due to their consistency for different data sampling rates and smaller standard errors, we recommend the OLS model when using pressure data. Since there was not a significant difference between OLS and GLS in terms of residual analysis, we examined other error models [9]. The results were better than GLS, but were not improved enough to supplant OLS as the recommended error model for pressure data.
We do note that the standard error for log 10 (τ 1 ) and log 10 (E 1 ) for both OLS and GLS cases is in general larger relative to the parameter estimate itself than for the other parameters; this is expected since the pressure model output is less sensitive to τ 1 and E 1 . Thus, we do not have as much confidence in the estimate for τ 1 and E 1 as in the other parameter estimates. See [9,Sec. 4.1] for a broader discussion of the sensitivity of the pressure model output to parameters.

Shear data results
Shear results are shown in Tables 5-6 and Figure 6. We report only the results for data sampled at 512 Hz, since the results were consistent with those found for data sampled at 1024 Hz (all results are available in [9]). In all cases, the model fits to data are good. The standard errors for the GLS case is larger than in the OLS case. In the bottom row of Figure 6, the residuals versus model plots are not noticeably different between the OLS and GLS cases. The initial indication is that we have more confidence in the OLS results. However, the residual versus time plots (middle row of Figure 6) raise cause for concern. In the OLS residual versus time plots, there is a noticeable "fan" structure for early times. However, for the GLS error model, the residual versus time plots do not show a fan structure and are fairly randomly distributed. Since this indicates that the OLS error model may not be correct, we are inclined to recommend the GLS error model in the shear case so that we do not mistakenly overstate our confidence in the parameter estimates, which we could do if we used the parameter estimates from the possibly-wrong OLS case.
For the shear case we see that the standard error for G 1 is on the same order of magnitude as the parameter estimate itself for both the OLS and GLS results, and the standard error for τ 1 in the GLS case is also on the same order of magnitude as the estimate for τ 1 . This is consistent with the sensitivity results, where in the shear case the model is less sensitive to G 1 and τ 1 than to G, ζ 1 , and A. Again, more details on the shear model output sensitivity results are available in [9].

Discussion and Future Work
We have developed an updated one-dimensional viscoelastic model for tissue and have used experimental data from a simple homogeneous gel phantom to test the ability of our model to describe wave propagation in the medium. The data were generated from a drop experiment designed to produce oscillations in the gel of a magnitude comparable to that produced by blood flow in a stenosed coronary artery impacting the vessel wall, a disturbance which results in pressure and shear waves propagating away from the vessel walls downstream of the blockage. In our inverse problem results as discussed in Section 4.2, we have shown an ability to consistently model the wave propagation using different error models and at different data sampling frequencies, obtaining good fits to data in all of our inverse problems. In addition to a good fit, though, we also examined statistical properties of the parameter estimators as well as residual plots to gain more insight into the proper error model for the pressure and shear data sets. This is necessary, since a correct error model is essential in order to apply the asymptotic error theory properly and thus obtain correct confidence intervals. For the pressure case, we prefer the absolute error model (OLS) over the relative error model (GLS) since the residual plots indicated no noticeable difference between the two models while the OLS parameter estimates were more consistent and had generally smaller corresponding standard errors. For the shear case, we recommend taking the more conservative route and using the GLS parameter estimates; even though the GLS estimates had larger standard errors, there were indications from the residual versus time plots for OLS that the OLS model is not correct.
We are currently examining a two dimensional model and corresponding experimental configurations. Experiments are currently in progress to produce a two dimensional wave from different points in the medium and with different detection points along the outer wall of the phantom. It is conceivable that the one dimensional parameters could be used as a rough first approximation in a corresponding two dimensional code, which would allow us to focus on trying to determine the location of the wave generation in the medium. Also, these parameter values could be used in a model of wave propagation in another conceptual device designed to mimic a constricted artery and the waves that result from passing fluid through a constricted pipe in the center of the medium. Therefore in the slightly longer term, we will also likely need to conduct an inverse problem using a two dimensional model and corresponding data. These one dimensional results will provide a starting point for parameters in that inverse problem, hopefully decreasing runtime and the time it takes to find viable parameters. The same issues discussed here (sensitivity to parameters, data frequency, number of relaxation times) will again be of concern for the two dimensional problem. Future efforts will also involve scaling up all these experiments to larger phantoms and then to some sort of actual tissue sample experiments.
Overall, we have successfully demonstrated the ability of the pressure and shear mathematical models to accurately describe the data from laboratory experiments. A linear viscoelastic constitutive relationship, i.e., (15b) and (16b), was adequate. This is a significant achievement, as all the work previously discussed was limited to inverse problems on simulated data or data that was not from the impulse-type experiments.