Skip to main content

Advertisement

Log in

Sonic Estimation of Elasticity via Resonance: A New Method of Assessing Hemostasis

  • Published:
Annals of Biomedical Engineering Aims and scope Submit manuscript

Abstract

Uncontrolled bleeding threatens patients undergoing major surgery and in care for traumatic injury. This paper describes a novel method of diagnosing coagulation dysfunction by repeatedly measuring the shear modulus of a blood sample as it clots in vitro. Each measurement applies a high-energy ultrasound pulse to induce a shear wave within a rigid walled chamber, and then uses low energy ultrasound pulses to measure displacements associated with the resonance of that shear wave. Measured displacements are correlated with predictions from finite difference time domain models, with the best fit corresponding to the modulus estimate. In our current implementation each measurement requires 62.4 ms. Experimental data was analyzed using a fixed-viscosity algorithm and a free-viscosity algorithm. In experiments utilizing human blood induced to clot by exposure to kaolin, the free-viscosity algorithm quantified the shear modulus of formed clots with a worst-case precision of 2.5%. Precision was improved to 1.8% by utilizing the fixed-viscosity algorithm. Repeated measurements showed a smooth evolution from liquid blood to a firm clot with a shear modulus between 1.4 and 3.3 kPa. These results show the promise of this technique for rapid, point of care assessment of coagulation.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11

Similar content being viewed by others

References

  1. Bernal, M., J. L. Gennisson, P. Flaud, and M. Tanter. Shear wave elastography quantification of blood elasticity during clotting. Ultrasound Med. Biol. 38:2218–2228, 2012.

    Article  PubMed  Google Scholar 

  2. Bernal, M., J.-L. Gennisson, P. Flaud, and M. Tanter. Two Dimension (2D) elasticity maps of coagulation of blood using SuperSonic Shearwave Imaging. Acoustics 2012. 2012.

  3. Bux, J. Transfusion-related acute lung injury (TRALI): a serious adverse event of blood transfusion. Vox Sang. 89:1–10, 2005.

    Article  CAS  PubMed  Google Scholar 

  4. Engoren, M. C., R. H. Habib, A. Zacharias, T. A. Schwann, C. J. Riordan, and S. J. Durham. Effect of blood transfusion on long-term survival after cardiac operation. Ann. Thorac. Surg. 74:1180–1186, 2002.

    Article  PubMed  Google Scholar 

  5. Ferraris, V. A., S. P. Saha, J. H. Oestreich, H. K. Song, T. Rosengart, T. B. Reece, et al. 2012 update to the society of thoracic surgeons guideline on use of antiplatelet drugs in patients having cardiac and noncardiac operations. Ann. Thorac. Surg. 94:1761–1781, 2012.

    Article  PubMed  Google Scholar 

  6. FibScreen, Platelet Function Test for ReoRox, MRX 1917. Nyköping, Sweden: MEDIROX AB, 2012.

  7. Franchini, M., M. Montagnana, E. J. Favaloro, and G. Lippi. The bidirectional relationship of cancer and hemostasis and the potential role of anticoagulant therapy in moderating thrombosis and cancer spread. Semin. Thromb. Hemost. 35:644–653, 2009.

    Article  CAS  PubMed  Google Scholar 

  8. Ganter, M. T., and C. K. Hofer. Coagulation monitoring: current techniques and clinical use of viscoelastic point-of-care coagulation devices. Anesth. Analg. 106:1366–1375, 2008.

    Article  PubMed  Google Scholar 

  9. Hess, J. R., and J. H. Lawson. The coagulopathy of trauma versus disseminated intravascular coagulation. J. Trauma 60:S12–S19, 2006.

    Article  PubMed  Google Scholar 

  10. Hoffman, M., and D. M. Monroe. A cell-based model of hemostasis. Thromb. Haemost. 85(6):958–965, 2001.

    CAS  PubMed  Google Scholar 

  11. Huang, C. C., P. Y. Chen, and C. C. Shih. Estimating the viscoelastic modulus of a thrombus using an ultrasonic shear-wave approach. Med. Phys. 40:042901, 2013.

    Article  PubMed  Google Scholar 

  12. Jensen, J. A., and N. B. Svendsen. Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 39:262–267, 1992.

    Article  CAS  PubMed  Google Scholar 

  13. Lang, T., A. Bauters, S. L. Braun, B. Pötzsch, K.-W. von Pape, H.-J. Kolde, and M. Lakner. Multi-centre investigation on reference ranges for ROTEM thromboelastometry. Blood Coagul. Fibrinolysis 16(4):301–310, 2005.

    Article  PubMed  Google Scholar 

  14. Lorente, J. A., L. J. García-Frade, L. Landín, R. de Pablo, C. Torrado, E. Renes, and A. García-Avello. Time course of hemostatic abnormalities in sepsis and its relation to outcome. Chest 103:1536–1542, 1993.

    Article  CAS  PubMed  Google Scholar 

  15. Mauldin, F. W., F. Viola, T. C. Hamer, E. M. Ahmed, S. B. Crawford, D. M. Haverstick, et al. Adaptive force sonorheometry for assessment of whole blood coagulation. Clin. Chim. Acta 411(9–10):638–644, 2010.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Mauldin, F. W., F. Viola, and W. F. Walker. Reduction of echo decorrelation via complex principal component filtering. Ultrasound Med. Biol. 35:1325–1343, 2009.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Murphy, S., and F. H. Gardner. Platelet preservation: effect of storage temperature on maintenance of platelet viability—deleterious effect of refrigerated storage. N. Engl. J. Med. 280:1094–1098, 1969.

    Article  CAS  PubMed  Google Scholar 

  18. Nahirnyak, V. M., S. W. Yoon, and C. K. Holland. Acousto-mechanical and thermal properties of clotted blood. J. Acoust. Soc. Am. 119:3766–3772, 2006.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Orescanin, M., Y. Wang, and M. F. Insana. 3-D FDTD simulation of shear waves for evaluation of complex modulus imaging. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 58(2):389–398, 2011.

    Article  PubMed  Google Scholar 

  20. Sarvazyan, A. P., O. V. Rudenko, S. D. Swanson, J. B. Fowlkes, and S. Y. Emelianov. Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics. Ultrasound Med. Biol. 24(9):1419–1435, 1998.

    Article  CAS  PubMed  Google Scholar 

  21. Scarpelini, S., S. G. Rhind, B. Nascimento, H. Tien, P. N. Shek, H. T. Peng, et al. Normal range values for thromboelastography in healthy adult volunteers. Braz. J. Med. Biol. Res. 42(12):1210–1217, 2009.

    Article  CAS  PubMed  Google Scholar 

  22. Scharf, Louis L. Statistical Signal Processing. Addison-Wesley: Reading, 1991.

    Google Scholar 

  23. Schmitt, C., A. Hadj Henni, and G. Cloutier. Characterization of blood clot viscoelasticity by dynamic ultrasound elastography and modeling of the rheological behavior. J. Biomech. 44(4):622–629, 2011.

    Article  PubMed  Google Scholar 

  24. Shung, K. K., and G. A. Thieme. Ultrasonic Scattering in Biological Tissues. Boca Raton: CRC Press, 1993.

    Google Scholar 

  25. Smith, F. B., A. J. Lee, F. G. Fowkes, J. F. Price, A. Rumley, and G. D. Lowe. Hemostatic factors as predictors of ischemic heart disease and stroke in the edinburgh artery study. Arterioscler. Thromb. Vasc. Biol. 17(11):3321–3325, 1997.

    Article  CAS  PubMed  Google Scholar 

  26. Sneddon, I. N. Fourier Transforms. New York: Dover Publications, 1995.

    Google Scholar 

  27. Solomon, C., H. Schöchl, M. Ranucci, U. Schött, and C. J. Schlimp. Comparison of fibrin-based clot elasticity parameters measured by free oscillation rheometry (ReoRox ®) versus thromboelastometry (ROTEM ®). Scand. J. Clin. Lab. Investig. 75:239, 2015.

    Article  Google Scholar 

  28. Torr, G. R. The acoustic radiation force. Am. J. Phys. 52(5):402–408, 1984.

    Article  Google Scholar 

  29. Viola, F., M. D. Kramer, M. B. Lawrence, J. P. Oberhauser, and W. F. Walker. Sonorheometry: a noncontact method for the dynamic assessment of thrombosis. Ann. Biomed. Eng. 32:696–705, 2004.

    Article  PubMed  Google Scholar 

  30. Viola, F., F. W. Mauldin, X. Lin-Schmidt, D. M. Haverstick, M. B. Lawrence, and W. F. Walker. A novel ultrasound-based method to evaluate hemostatic function of whole blood. Clin. Chim. Acta 411:106–113, 2010.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Viola, F., and W. F. Walker. A comparison of the performance of time-delay estimators in medical ultrasound. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 50(4):392–401, 2003.

    Article  PubMed  Google Scholar 

  32. Virieux, J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics 51(4):889–901, 1986.

    Article  Google Scholar 

  33. Walker, W. F., and G. E. Trahey. A fundamental limit on delay estimation using partially correlated speckle signals. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 42:301–308, 1995.

    Article  Google Scholar 

  34. Yarnell, J. W., I. A. Baker, P. M. Sweetnam, D. Bainton, J. R. O’Brien, P. J. Whitehead, and P. C. Elwood. Fibrinogen, viscosity, and white blood cell count are major risk factors for ischemic heart disease. The caerphilly and speedwell collaborative heart disease studies. Circulation 83(3):836–844, 1991.

    Article  CAS  PubMed  Google Scholar 

  35. Yee, K. S. Numerical solution of initial boundary value problems involving maxwells equations in isotropic media. IEEE Trans. Antennas Propag. 14(3):302–307, 1966.

    Article  Google Scholar 

Download references

Acknowledgments

The authors would like to thank Andrew Homyk his assistance in the preparation of Figs. 2, 3, and 4. We thank Kiev Blasier, Elisa Ferrante, Caroline Wang, and Bob Fehnel for performing the human blood experiments. We thank Cindy Lloyd for formulating the kaolin. We gratefully acknowledge design work by Andy Homyk, Tim Higgins, and Aaron Buchannan at HemoSonics, and Frank Reagan, Lei Zong, and Jeff Gunnarsson at Key Technologies Inc. We acknowledge Francesco Viola for his contributions to the architecture of both the instrument and the cartridge. We thank Timothy J. Fischer and Thomas B. Givens for their helpful suggestions on the manuscript. The authors also thank Dr. Michael F. Insana and Ms. Yue Wang of the University of Illinois for sharing their insights into modeling shear wave propagation via the Finite Difference Time Domain method. This work was supported by NIH Grants 1R43CA180318-01A1 and 2R44HL103030-02A1, and the investors of HemoSonics, LLC.

Conflict of interest

Author William F. Walker is a major shareholder in and employee of HemoSonics, LLC, the sole licensee and assignee of patents protecting the technology described in this manuscript. Author F. Scott Corey is an employee and shareholder in KeyTech Incorporated, a design partner of and shareholder in HemoSonics, LLC.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to William F. Walker.

Additional information

Associate Editor Jane Grande-Allen oversaw the review of this article.

Appendices

Appendix 1: Finite Difference Formulation

Our formulation of the FDTD model is built upon a velocity – stress formulation of the shear wave equation.19,32 We expand upon traditional formulations by including a term to account for viscous loss.

$$\rho \frac{\partial }{\partial t}\vec{v}\left( {\vec{x},t} \right) = \nabla \cdot \sigma \left( {\vec{x},t} \right) + \vec{F}\left( {\vec{x},t} \right)$$
(5)
$$\frac{\partial }{\partial t}\sigma \left( {\vec{x},t} \right) = \left( {\mu + \eta \frac{\partial }{\partial t}} \right)\left[ {\left( {\nabla \vec{v}\left( {\vec{x},t} \right)} \right) + \left( {\nabla \vec{v}\left( {\vec{x},t} \right)} \right)^{T} } \right].$$
(6)

Our analysis assumes an axisymmetric test chamber, which suggests the use of cylindrical coordinates. Representing Eq. (5) in cylindrical coordinates and expanding the vector velocity into its constituent components yields:

$$\frac{{\partial \sigma_{rr} }}{\partial r} + \frac{1}{r}\frac{{\partial \sigma_{r\phi } }}{\partial \phi } + \frac{{\partial \sigma_{rz} }}{\partial z} + \frac{1}{r}\left( {\sigma_{rr} - \sigma_{\phi \phi } } \right) + F_{r} = \rho \frac{{\partial v_{r} }}{\partial t}$$
(7)
$$\frac{{\partial \sigma_{r\phi } }}{\partial r} + \frac{1}{r}\frac{{\partial \sigma_{\phi \phi } }}{\partial \phi } + \frac{{\partial \sigma_{\phi z} }}{\partial z} + \frac{2}{r}\sigma_{r\phi } + F_{\phi } = \rho \frac{{\partial v_{\phi } }}{\partial t}$$
(8)
$$\frac{{\partial \sigma_{{r{\text{z}}}} }}{\partial r} + \frac{1}{r}\frac{{\partial \sigma_{{\phi {\text{z}}}} }}{\partial \phi } + \frac{{\partial \sigma_{zz} }}{\partial z} + \frac{1}{r}\sigma_{rz} + F_{z} = \rho \frac{{\partial v_{z} }}{\partial t}.$$
(9)

We further simplify by recognizing that the only body force, the applied ultrasound radiation force, acts only in the \(z\) direction. We set \(F_{\phi } = F_{r} = 0\). Since our test chamber and the applied radiation force are entirely axisymmetric, we set all dependencies upon \(\phi\) equal to zero. Applying these simplifications to Eq. (7) yields:

$$\frac{{\partial \sigma_{rz} }}{\partial z} = \rho \frac{{\partial^{2} u_{r} }}{{\partial t^{2} }}$$
(10)
$$0 = 0$$
(11)
$$\frac{{\partial \sigma_{{r{\text{z}}}} }}{\partial r} + \frac{1}{r}\sigma_{rz} + F_{z} = \rho \frac{{\partial^{2} u_{z} }}{{\partial t^{2} }}.$$
(12)

We follow a similar strategy to expand Eq. (6), to yield:

$$\frac{{\partial \sigma_{{r{\text{z}}}} }}{\partial t} = \left( {\mu + \eta \frac{\partial }{\partial t}} \right)\left( {\frac{{\partial v_{r} }}{\partial z} + \frac{{\partial v_{z} }}{\partial r}} \right).$$
(13)

Collecting Eqs. (10)–(12) and Eq. (13) yields:

$$\rho \frac{{\partial v_{r} }}{\partial t} = \frac{{\partial \sigma_{rz} }}{\partial z}$$
(14)
$$\rho \frac{{\partial v_{z} }}{\partial t} = \frac{{\partial \sigma_{{r{\text{z}}}} }}{\partial r} + \frac{1}{r}\sigma_{rz} + F_{z}$$
(15)
$$\frac{{\partial \sigma_{{r{\text{z}}}} }}{\partial t} = \left( {\mu + \eta \frac{\partial }{\partial t}} \right)\left( {\rho \frac{{\partial v_{r} }}{\partial t} + \rho \frac{{\partial v_{z} }}{\partial t}} \right).$$
(16)

Equations (14)–(16) form a system of partial differential equations that predict how radiation force will induce shear waves and how those induced shear waves will interact. This system of equations is amenable to a finite difference solution using a staggered grid approach, similar to the Yee method.35 The finite difference representation of Eqs. (14)–(16) are:

$$v_{{r_{j,k}^{i + 1} }} = v_{{r_{j,k}^{i} }} + \frac{\Delta t}{\rho \, \Delta z}\left( {\sigma_{{rz_{{j,k + \frac{1}{2}}}^{i} }} - \sigma_{{rz_{{j,k - \frac{1}{2}}}^{i} }} } \right)$$
(17)
$$v_{{z_{{j + \frac{1}{2},k + \frac{1}{2}}}^{i + 1} }} = v_{{z_{{j + \frac{1}{2},k + \frac{1}{2}}}^{i} }} + \frac{\Delta t}{\rho }\left( {\frac{1}{\Delta r}\left( {\sigma_{{rz_{{j + 1,k + \frac{1}{2}}}^{i} }} - \sigma_{{rz_{{j,k + \frac{1}{2}}}^{i} }} } \right) + \frac{1}{{2\left( {j + \frac{1}{2}} \right)\Delta r}}\left( {\sigma_{{rz_{{j,k + \frac{1}{2}}}^{i} }} - \sigma_{{rz_{{j + 1,k + \frac{1}{2}}}^{i} }} } \right) + F_{{z_{{j + + \frac{1}{2},k + \frac{1}{2}}}^{i} }} } \right)$$
(18)
$$\sigma_{{rz_{{j,k + \frac{1}{2}}}^{i + 1} }} = \sigma_{{rz_{{j,k + \frac{1}{2}}}^{i} }} + \frac{\mu \, \Delta t}{\Delta z}\left( {v_{{r_{j,k + 1}^{i} }} - v_{{r_{j,k}^{i} }} } \right) + \frac{\mu \, \Delta t}{\Delta z}\left( {v_{{r_{{j + \frac{1}{2},k + + \frac{1}{2}}}^{i} }} - v_{{r_{{j - \frac{1}{2},k + \frac{1}{2}}}^{i} }} } \right) + \frac{\eta }{\Delta z}\left( {v_{{r_{j,k + 1}^{i + 1} }} - v_{{r_{j,k + 1}^{i} }} - v_{{r_{j,k}^{i + 1} }} + v_{{r_{j,k}^{i} }} } \right) + \frac{\eta }{\Delta z}\left( {v_{{z_{{j + \frac{1}{2},k + \frac{1}{2}}}^{i + 1} }} - v_{{z_{{j + \frac{1}{2},k + \frac{1}{2}}}^{i} }} - v_{{z_{{j - \frac{1}{2},k + \frac{1}{2}}}^{i + 1} }} + v_{{z_{{j - \frac{1}{2},k + \frac{1}{2}}}^{i} }} } \right).$$
(19)

This formulation was implemented in MATLAB (Natick, MA) using the staggered grid approach shown in Fig. 12 with a default sampling interval of 100 \(\mu\)m. Each computational cycle consisted of two steps. In the first, the velocity components were computed via Eqs. (17) and (18). In the next step shear was computed using Eq. (19). time–dependent velocities were integrated numerically using cumulative trapezoidal integration ( cumtrapz in MATLAB) to determine the time–varying displacement.

Figure 12
figure 12

FDTD model staggered grid geometry.

Appendix 2: Analytical Model

This analysis is based upon the Cauchy-Navier Equation for linear elasticity. We utilize a formulation incorporating viscoelastic material properties.

$$\left( {\left( {\lambda + \mu + \left( {\eta_{p} + \eta_{s} } \right)\frac{d}{dt} \cdot } \right)\nabla \left( {\nabla \cdot \vec{u}\left( {\vec{x},t} \right)} \right)} \right) + .\left( {\left( {\mu + \eta_{s} \frac{d}{dt} \cdot } \right)\nabla^{2} \vec{u}\left( {\vec{x},t} \right)} \right) - \rho \frac{{d^{2} }}{{dt^{2} }}\vec{u}\left( {\vec{x},t} \right) = \vec{F}\left( {\vec{x},t} \right).$$
(20)

To simplify analysis we assume an infinitely long cylindrical test chamber with radiation force applied to the central axis. We apply these simplifications to the components of Eq. (20):

$$\nabla \left( {\nabla \cdot \vec{u}\left( {\vec{x},t} \right)} \right) = \nabla \left( {\frac{1}{r}\frac{d}{dr}\left( {ru_{r} } \right) + \frac{1}{r}\frac{d}{d\phi }u_{\phi } + \frac{d}{dz}u_{z} } \right)$$
(21)
$$\nabla^{2} \vec{u}\left( {\vec{x},t} \right) = \left[ {\nabla^{2} u_{r} - \frac{1}{{r^{2} }}u_{r} - \frac{2}{{r^{2} }}\frac{d}{d\phi }u_{\phi } } \right]\hat{r} + \left[ {\nabla^{2} u_{\phi } - \frac{1}{{\phi^{2} }}u_{\phi } - \frac{2}{{r^{2} }}\frac{d}{d\phi }u_{r} } \right]\hat{\phi } + \left[ {\nabla^{2} u_{z} } \right]\hat{z}.$$
(22)

Because of cylindrical symmetry and the assumed infinite length of the chamber we can set all derivatives with respect to \(\phi\) and \(z\) equal to zero.

$$\nabla \left( {\nabla \cdot \vec{u}\left( {\vec{x},t} \right)} \right) = \nabla \left( {\frac{1}{r}\frac{d}{dr}\left( {ru_{r} } \right)} \right)\hat{r}$$
(23)
$$\nabla^{2} \vec{u}\left( {\vec{x},t} \right) = \left[ {\nabla^{2} u_{r} - \frac{1}{{r^{2} }}u_{r} } \right]\hat{r} + 0\hat{\phi } + \left[ {\nabla^{2} u_{z} } \right]\hat{z}.$$
(24)

Expanding gradients and Laplacians as appropriate, and setting derivatives with respect to \(\phi\) and \(z\) equal to zero yields:

$$\nabla \left( {\nabla \cdot \vec{u}\left( {\vec{x},t} \right)} \right) = \left[ {\frac{d}{dr}\left( {\frac{1}{r}\frac{d}{dr}\left( {ru_{r} } \right)} \right)} \right]\hat{r} + 0\hat{\phi } + 0\hat{z}$$
(25)
$$\nabla^{2} \vec{u}\left( {\vec{x},t} \right) = \left[ {\frac{1}{r}\frac{d}{dr}\left( {r\frac{d}{dr}u_{r} } \right) + \frac{1}{{r^{2} }}u_{r} } \right]\hat{r} + 0\hat{\phi } + \left[ {\frac{1}{r}\frac{d}{dr}\left( {r\frac{d}{dr}u_{z} } \right)} \right]\hat{z}.$$
(26)

This result shows that there are no cross terms between the radial and axial displacements. Our system measures axial displacement, so we discard the radial displacement term after substituting Eqs. (10) and (11) into Eq. (20), and assuming that the force is applied impulsively at time zero and is a spatial impulse at the center of the measurement cylinder.

$$\left( {\left( {\mu + \eta_{s} \frac{d}{dt} \cdot } \right)\left( {\frac{{d^{2} }}{{dr^{2} }}u_{z} \left( {r,t} \right) + \frac{1}{r}\frac{d}{dr}u_{z} \left( {r,t} \right)} \right)} \right) - \rho \frac{{d^{2} }}{{dt^{2} }}u_{z} \left( {r,t} \right) = \delta \left( t \right)\delta \left( r \right)/r.$$
(27)

We now utilize the Hankel Transform as defined in Eq. (28)26 to solve this partial differential equation.

$$\tilde{u}\left( {\xi_{i} ,t} \right) = \mathop \smallint \limits_{0}^{{R_{c} }} ru\left( {r,t} \right)J_{0} \left( {\xi_{i} r} \right)dr,$$
(28)

where \(R_{c}\) is the radius of the cylinder, \(J_{0} \left( x \right)\) is 0th order Bessel function of the first kind, and \(\xi_{i} = \lambda_{i} /R_{c}\) where \(\varvec{\lambda}_{\varvec{i}}\) is the ith root of the 0th order Bessel function of the first kind. Performing the Hankel Transform on Eq. (27) yields:

$$- \mu \xi_{i}^{2} \tilde{u}\left( {\xi_{i} ,t} \right) - \eta_{s} \xi_{i}^{2} \frac{d}{dt}\tilde{u}\left( {\xi_{i} ,t} \right) - \rho \frac{{d^{2} }}{{dt^{2} }}\tilde{u}\left( {\xi_{i} ,t} \right) = \delta \left( t \right).$$
(29)

The application of the Hankel Transform has reduced a partial differential equation to an ordinary differential equation. Assuming that \(\hat{u}\left( {\xi_{i} ,0} \right) = 0\) and that \(\frac{d}{dt}\tilde{u}\left( {\xi_{i} ,t} \right) = 0\) then:

$$\tilde{u}\left( {\xi_{i} ,t} \right) = \frac{1}{\sqrt A }e^{{ - t\eta_{s} \xi_{i}^{2} /2\rho }} \left( {e^{ - t\sqrt A /2\rho } - e^{t\sqrt A /2\rho } } \right).$$
(30)

where \(A = \xi_{i}^{2} \left( {\eta_{s}^{2} \xi_{i}^{2} - 4\mu \rho } \right)\). The finite Hankel Transform can be inverted using Eq. (31) [Equation 45 on page 83 of Ref. 26].

$$f\left( x \right) = \frac{2}{{R_{c}^{2} }}\mathop \sum \limits_{i = 1}^{\infty } \tilde{f}\left( {\xi_{i} } \right)\frac{{J_{0} \left( {\xi_{i} x} \right)}}{{\left[ {J_{1} \left( {\xi_{i} R_{c} } \right)} \right]^{2} }}.$$
(31)

Applying the inversion formula of Eq. (31) to Eq. (30) yields an expression for the time–dependent displacement.

$$u\left( {r,t} \right) = \frac{2}{{R_{c}^{2} }}\mathop \sum \limits_{i = 1}^{\infty } \frac{{J_{0} \left( {\xi_{i} r} \right)}}{{\sqrt A \left[ {J_{1} \left( {\xi_{i} R_{c} } \right)} \right]^{2} }}e^{{ - t\eta_{s} \xi_{i}^{2} /2\rho }} \times \left( {e^{ - t\sqrt A /2\rho } - e^{t\sqrt A /2\rho } } \right).$$
(32)

In comparing Eq. (33) to results from a simplified FDTD model we see that analytical model consistently exhibits high frequency content that is not visible in the FDTD model. We believe that this discrepancy is a result of the finite spatial sampling interval of the FDTD model. We have found empirically that the analytical model is in good agreement with the FDTD model when we include only terms up to those for which there is at least 1.5 samples per half cycle of the Bessel function. This truncated model, shown as Eq. (33) shows excellent agreement with the simplified FDTD model.

$$\hat{u}\left( {r,t} \right) = \frac{2}{{R_{c}^{2} }}\mathop \sum \limits_{i = 1}^{{2R_{c} /3\Delta x}} \frac{{J_{0} \left( {\xi_{i} r} \right)}}{{\sqrt A \left[ {J_{1} \left( {\xi_{i} R_{c} } \right)} \right]^{2} }}e^{{ - t\eta_{s} \xi_{i}^{2} /2\rho }} \times \left( {e^{ - t\sqrt A /2\rho } - e^{t\sqrt A /2\rho } } \right),$$
(33)

where \(\Delta x\) is the spatial sampling interval of the FDTD model. The terms of the sum take on two different forms depending upon whether \(A > 0\), in which case the system is overdamped, or \(A < 0\), in which case the system is underdamped. In the underdamped case the first term of the series is a sine wave under a decaying exponential envelope. Assuming that the system has very low damping, the resonant frequency is:

$$f_{0} = {\raise0.7ex\hbox{${\sqrt A }$} \!\mathord{\left/ {\vphantom {{\sqrt A } {2\pi 2\rho }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${2\pi 2\rho }$}} = {\raise0.7ex\hbox{${ - i\sqrt {\xi_{1}^{2} \left( {\eta_{s}^{2} \xi_{1}^{2} - 4\mu \rho } \right)} }$} \!\mathord{\left/ {\vphantom {{ - i\sqrt {\xi_{1}^{2} \left( {\eta_{s}^{2} \xi_{1}^{2} - 4\mu \rho } \right)} } {4\pi \rho }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${4\pi \rho }$}} \approx {\raise0.7ex\hbox{${ - i\sqrt {\xi_{1}^{2} \left( { - 4\mu \rho } \right)} }$} \!\mathord{\left/ {\vphantom {{ - i\sqrt {\xi_{1}^{2} \left( { - 4\mu \rho } \right)} } {4\pi \rho }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${4\pi \rho }$}} = \frac{{\xi_{1} }}{2\pi }\sqrt {\frac{\mu }{\rho }} = \frac{2.4048}{{2\pi R_{c} }}\sqrt {\frac{\mu }{\rho }} ,$$
(34)

where the factor of \(2\varvec{\pi}\) is required to convert from radial to temporal frequency. The shear modulus can be estimated from the resonant frequency:

$$\mu = \rho \left( {\frac{{2\pi f_{0} R_{c} }}{2.4048}} \right)^{2} .$$
(35)

The concept of resonant frequency is only meaningful when the resonant chamber is under-damped. The system is under-damped when \(A < 0\), or equivalently when:

$$\mu > \frac{{\eta_{s}^{2} \xi_{i}^{2} }}{4\rho }.$$
(36)

For our geometry assuming a blood density of 1060 kg/m3 and a shear viscosity of approximately 0.25 Pa s, the system will be under-damped for moduli greater than 19.3 Pa.

Appendix 3: Cramér–Rao Lower Bound on SEER

We begin our analysis from the formulation derived by Scharf in Ref. 22 The lower bound on the variance \(C_{ii}\) of estimates of the parameter \(\theta_{i}\) is:

$$C_{ii} \ge \left( {J^{ - 1} } \right)_{ii}$$
(37)

where \(J\) is the Fisher information matrix. The elements of the Fisher information matrix, for a single realization of data, are defined by Equation 6.134 of Ref. 22:

$$J_{ij} = \frac{1}{2}tr\left[ {R^{ - 1} \frac{\partial R}{{\partial \theta_{i} }}R^{ - 1} \frac{\partial R}{{\partial \theta_{j} }}} \right] + \frac{{\partial m^{t} }}{{\partial \theta_{i} }}R^{ - 1} \frac{\partial m}{{\partial \theta_{j} }},$$
(38)

where \(tr\left[ X \right]\) is the trace of matrix \(X\), \(R\) is the expected covariance matrix of the observed data, and \(m\) is the expected mean of the observed data. In SEER the observed data is the vector of time–displacement estimates. We assume that the noise in the displacement estimates is generally uncorrelated so that the covariance matrix can be modeled as an identity matrix scaled by the variance of the displacement estimates.

$$R = {\sigma_{x}^{2} }I,$$
(39)

where \(\sigma_{x}^{2}\) is the variance of the displacement estimates. Assuming that the error covariance does not depend upon the parameter (either \(\theta_{i}\) or \(\theta_{j}\)) and that the covariance can be represented as in (39), then (38) can be simplified as:

$$J_{{ij}} = \frac{{\partial m^{t} }}{{\partial \theta _{i} }}R^{{ - 1}} \frac{{\partial m}}{{\partial \theta _{j} }} = \frac{{\partial m^{t} }}{{\partial \theta _{i} }}\left[ {{{\sigma _{x}^{2} }}I} \right]^{{ - 1}} \frac{{\partial m}}{{\partial \theta _{j} }} = \frac{1}{{\sigma _{x}^{2} }}\frac{{\partial m^{t} }}{{\partial \theta _{i} }}\frac{{\partial m}}{{\partial \theta _{j} }}.$$
(40)

For the free-viscosity estimator we estimate both the viscosity and the modulus, so the Fisher Information Matrix is:

$$J = \left[ {\begin{array}{*{20}c} {\frac{1}{{\sigma_{x}^{2} }}\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\mu } }}} & {\frac{1}{{\sigma_{x}^{2} }}\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\eta } }}} \\ {\frac{1}{{\sigma_{x}^{2} }}\frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\mu } }}} & {\frac{1}{{\sigma_{x}^{2} }}\frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\eta } }}} \\ \end{array} } \right] = \frac{1}{{\sigma_{x}^{2} }}\left[ {\begin{array}{*{20}c} {\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\mu } }}} & {\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\eta } }}} \\ {\frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\mu } }}} & {\frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\eta } }}} \\ \end{array} } \right].$$
(41)

Substituting Eqs. (41) into (37) yields the following lower bounds on estimation in the free-viscosity algorithm:

$$C_{\mu \mu } \ge \frac{{\sigma_{x}^{2} }}{{\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\mu } }}\frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\eta } }} - \frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\eta } }}\frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\mu } }}}} \frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\eta } }}$$
(42)
$$C_{\eta \eta } \ge \frac{{\sigma_{x}^{2} }}{{\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\mu } }}\frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\eta } }} - \frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\eta } }}\frac{{\partial m^{t} }}{{\partial \theta_{\eta } }}\frac{\partial m}{{\partial \theta_{\mu } }}}}\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\mu } }},$$
(43)

where \(\mu\) is modulus and \(\eta\) is viscosity. The mean signal value is the expected displacement for a given combination of modulus and viscosity. These bounds offer a lower limit on performance for any unbiased estimator. It is an optimistic estimate because it does not fully account for signal decorrelation associated with radiation force induced displacement.

The performance bound for the fixed-viscosity estimator is simpler as it has only one free parameter, yielding:

$$C_{\mu \mu } \ge {\raise0.7ex\hbox{${\sigma_{x}^{2} }$} \!\mathord{\left/ {\vphantom {{\sigma_{x}^{2} } {\left( {\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\mu } }}} \right)}}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\left( {\frac{{\partial m^{t} }}{{\partial \theta_{\mu } }}\frac{\partial m}{{\partial \theta_{\mu } }}} \right)}$}}.$$
(44)

The lower bound on the variance of modulus estimates for the fixed-viscosity algorithm is smaller than that for the free-viscosity estimator. This is not surprising because the fixed-viscosity estimator has only one degree of freedom. It should be noted that the total error for the fixed-viscosity algorithm may actually be higher however because of bias resulting from differences between the actual and assumed viscosity.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Corey, F.S., Walker, W.F. Sonic Estimation of Elasticity via Resonance: A New Method of Assessing Hemostasis. Ann Biomed Eng 44, 1405–1424 (2016). https://doi.org/10.1007/s10439-015-1460-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s10439-015-1460-y

Keywords

Navigation