Determining mechanical properties via ultrasound-induced resonance

Abstract

A device for estimating a mechanical property of a sample is disclosed herein. The device may include a chamber configured to hold the sample; a transmitter configured to transmit a plurality of waveforms, including at least one forcing waveform; and a transducer assembly operatively connected to the transmitter and configured to transform the transmit waveforms into ultrasound waveforms. The transducer assembly can also transmit and receive ultrasound waveforms into and out of the chamber, as well as transform at least two received ultrasound waveforms into received electrical waveforms. The device also includes a data processor that can receive the received electrical waveforms; estimate a difference in the received electrical waveforms that results at least partially from movement of the sample; and estimate a mechanical property of the sample by comparing at least one feature of the estimated difference to at least one predicted feature, wherein the at least one predicted feature is based on a model of an effect of the chamber wall. Finally, the device can also include a controller configured to control the timing of the ultrasound transmitter and data processor.

Description

BACKGROUND

Hemostasis is the physiological process that stops bleeding. Functional hemostasis requires the balanced engagement of plasma coagulation factors to initiate clotting, adequate fibrinogen to form the fibrin mesh, platelets to modulate factor function and to mechanically stiffen the fibrin mesh, and fibrinolytic enzymes to dissolve the clot when its useful life is over. Perturbation of any of these subsystems can disrupt hemostasis; either by impeding the cessation bleeding or by initiating clotting when it is unneeded. Disruptions of hemostasis contribute significantly to morbidity and mortality in patients suffering from heart disease, stroke, traumatic injury, cancer, and sepsis.

While hemostatic dysfunction impacts a broad range of medical conditions, it has been studied with particular intensity in cardiac surgery. Cardiac bypass surgery is associated with significant post-operative bleeding. This is caused by a combination of platelet damage by the bypass pump, factor and fibrinogen consumption associated with surgical trauma, and the occasional presence of residual anti-coagulant. A number of strategies are currently used to manage this dysfunction. The crudest strategy is a “shotgun therapy” approach; transfusing varying combinations of fresh frozen plasma, cryoprecipitate or fibrinogen concentrate, and platelet concentrate. This approach is often successful in controlling bleeding, however unnecessary transfusion carries a significant financial cost and increases patient morbidity and mortality. Recognition of the risks associated with excessive transfusion have led to increasingly specific and detailed guidelines to manage transfusion. These guidelines call for transfusion to be guided by point-of-care tests that assess hemostatic function in a timely and accurate manner.

A broad range of approaches have been proposed to fill the need for point-of-care hemostasis testing. These technologies can be separated into broad categories: clot-time assays, platelet-only tests, and viscoelastic tests. Clot-time assays can be implemented in simple systems, however a rapidly forming clot may not be a physiologically useful dot, so clot-time results may have limited clinical value. Moreover, clot-time assays generally operate on plasma, rather than whole blood, and therefore typically neglect the important interactions between plasma coagulation factors and platelets. Platelet-only tests provide useful information, but are also limited in that they neglect interactions between platelets and plasma coagulation factors. Viscoelastic tests have been shown to provide highly useful data. However, their operational complexity has traditionally limited their point-of-care utility. None of the currently available point-of-care tests adequately assess hemostatic function in a timely and accurate manner. A fast, accurate test is therefore still needed to fill this gap.

Other systems, methods, features and/or advantages will be or may become apparent to one with skill in the art upon examination of the following drawings and detailed description. It is intended that all such additional systems, methods, features and/or advantages be included within this description and be protected by the accompanying claims.

SUMMARY

A device for estimating a mechanical property of a sample is disclosed herein. The device may include a chamber configured to hold the sample; a transmitter configured to transmit a plurality of waveforms, including at least one forcing waveform; and a transducer assembly operatively connected to the transmitter and configured to transform the transmit waveforms into ultrasound waveforms. The transducer assembly can also transmit and receive ultrasound waveforms into and out of the chamber, as well as transform at least two received ultrasound waveforms into received electrical waveforms. The device also includes a data processor that can receive the received electrical waveforms; estimate a difference in the received electrical waveforms that results at least partially from movement of the sample; and estimate a mechanical property of the sample by comparing at least one feature of the estimated difference to at least one predicted feature, wherein the at least one predicted feature is based on a model of an effect of the chamber wall. Finally, the device can also include a controller configured to control the timing of the ultrasound transmitter and data processor.

In one implementation, the at least one predicted feature predicted by the data processor is based on a model of an aspect of the induced movement that is caused at least in part by a boundary effect of the chamber wall.

In another implementation, the at least one predicted feature predicted by the data processor is based on a model of an aspect of the induced movement that is caused at least in part by resonance within the chamber.

In yet another implementation, the at least one predicted feature predicted by the data processor is based on a model of an aspect of the induced movement that is caused at least in part by reflection of an induced shear wave from the chamber wall.

Also disclosed is a device for estimating a mechanical property of a sample. The device may include a chamber configured to hold the sample; a transmitter configured to transmit a plurality of waveforms, including at least one forcing waveform; and a transducer assembly operatively connected to the transmitter and configured to transform the transmit waveforms into ultrasound waveforms. The transducer assembly can also transmit and receive ultrasound waveforms into and out of the chamber, as well as transform at least two received ultrasound waveforms into received electrical waveforms. The device also includes a data processor that can receive the received electrical waveforms; estimate a difference in the received electrical waveforms that results at least partially from resonance of the sample; and characterize a mechanical property of the sample from at least one feature of the estimated difference. Finally, the device can also include a controller configured to control the timing of the ultrasound transmitter and data processor.

Additional embodiments, implementations, and/or examples are provided below as well.

BRIEF DESCRIPTION OF THE DRAWINGS

The following detailed description will be better understood when read in conjunction with the appended drawings, in which there is shown one or more of the multiple embodiments of the present disclosure. It should be understood, however, that the various embodiments of the present disclosure are not limited to the precise arrangements and instrumentalities shown in the drawings.

FIG. 1 is an example schematic of a device for measuring a mechanical property of a sample.

FIG. 2 is an example schematic of the signal processing and data acquisition process using the implementation of FIG. 1.

FIG. 3 shows representative experimental time-displacement curves from clotting human blood.

FIG. 4 shows representative time-modulus and time-viscosity curves estimated utilizing a computer model.

FIG. 5 shows an example contour plot of correlations between experimental time-displacement data and computer models across a range of moduli and viscosities.

FIG. 6 shows a representative time-modulus curves estimated utilizing a computer model with the viscosity held constant.

FIG. 7 shows an example of the use of correlation masking to remove outlying modulus estimates.

FIG. 8 depicts an example embodiment of the present disclosure wherein a focusing element is associated with the transducer assembly

FIG. 9 depicts an example embodiment of the present disclosure wherein an acoustic couplant is interposed between the transducer assembly and the test chamber.

FIG. 10 depicts an example embodiment of the present disclosure wherein a focusing element is associated with the transducer assembly and an acoustic couplant is interposed between the transducer assembly and the test chamber.

FIG. 11 depicts an example embodiment of the present disclosure wherein a focusing element is associated with the test chamber and an acoustic couplant is interposed between the transducer assembly and the test chamber.

FIG. 12 depicts an example embodiment of the present disclosure incorporating numerous additional elements.

FIG. 13 depicts an example embodiment of the present disclosure wherein the transducer assembly comprises two separate transducer elements.

FIG. 14 depicts an embodiment of the present disclosure wherein the transducer assembly comprises two separate transducer elements with one dedicated to waveform transmission and the other dedicated to waveform reception.

FIG. 15 shows the computational grid used for the Finite Difference Time Domain method described in the present disclosure.

FIG. 16 shows a representative time-modulus curve estimated utilizing an analytical model.

DETAILED DESCRIPTION

Assessing hemostatic function in a timely and accurate manner may be accomplished by measuring mechanical properties of a blood sample as it clots. For example, the shear modulus of a blood sample may be measured over time during the process of coagulation. Throughout this application, “shear modulus” is referred to as “modulus” interchangeably. The terms “stiff” and “stiffness” are likewise indicative of modulus.

In embodiments described in the current disclosure, a forcing ultrasound waveform is applied to a sample within a test chamber. This forcing waveform applies acoustic radiation force to the sample, thus inducing motion. This motion is impacted by the presence of the test chamber walls. Ultrasound sensing pulses are applied to the sample and differences in their echoes provide information about movement of the sample. Such differences may include phase changes or time shifts, either of which can be related to displacement. Finally, the features of these differences, such as the period of oscillation of sample motion, are compared to analytical or computational models to estimate mechanical properties of the sample.

Embodiments described in this disclosure may comprise elements shown in FIG. 1. A transmitter 112 emits electrical waveforms including at least one forcing waveform. These electrical waveforms are converted by transducer assembly 114 into ultrasound waveforms. The forcing waveform induces sample motion within test chamber 116. Returned ultrasound echoes are converted to electrical waveforms by transducer assembly 114. These electrical waveforms are analyzed by data processor 118 to estimate sample mechanical properties. The details of implementation of this and other embodiments may be understood through the detailed description that follows. The timing of transmission, reception, and data processing are controlled by the controller 110.

The present disclosure makes use of the phenomena of ultrasound radiation force, which may sometimes be referred to as acoustic radiation force. Ultrasound radiation force is a body force that acts in the same direction as the direction of propagation of the ultrasound wave. It is the result of a transfer of momentum that occurs as the travelling ultrasound wave is absorbed or reflected. Ultrasound radiation force is localized to the ultrasound beam; with the force magnitude proportional to the intensity of the ultrasound beam. Thus a focused ultrasound beam may be used to apply a localized acoustic radiation force field. In the present disclosure the applied radiation force field is generally smaller than the test chamber to which it is applied. The sample motion induced by the radiation force field will initially be localized to the region of the force field. Over time however the displacement field will propagate outward from the region of force application.

In embodiments described in the current disclosure, ultrasound radiation force is used to generate a shear wave within a blood sample. A shear wave is a mechanical wave in which the direction of particle displacement is perpendicular to the direction of wave motion. The shear wave of the present disclosure may be generated by directing an ultrasonic waveform of sufficient magnitude into the sample via a “forcing waveform.” The forcing waveform carries sufficient energy so that its absorption and reflection within the test sample generates acoustic radiation force. Acoustic radiation force is induced along the direction of ultrasound wave propagation and can be considered a body force acting upon some volume of the medium constrained by the dimensions of the ultrasound beam.

The induced shear wave will travel within the test chamber, reflecting from one or more of the walls. In some embodiments a single reflected shear wave may be detected and its arrival time at the ultrasound sensing beam may form the basis for estimating the modulus of the sample. For some combinations of chamber geometry and sample mechanical properties, it may not be possible to detect a reflected shear wave. For example the viscosity of the sample may be so high that the shear wave, while generated, has no easily discernable features and therefore no shear wave reflection from the wall can be readily measured. Nonetheless the sample motion resulting from acoustic radiation force will exhibit disturbances that result from the interaction of an induced shear wave with the chamber wall. The characteristics of this sample motion may form the basis for estimating the mechanical properties of the sample, even when the shear wave itself is not clearly resolved in the induced sample motion. In yet other cases the shear wave may reflect repeatedly within the chamber, generating a resonance that can form the basis for estimating the modulus of the sample.

Throughout this specification the terms test chamber, resonant chamber, resonance chamber, and chamber are used interchangeably with no loss in generality.

In the present disclosure the terms forcing waveform and forcing pulse will be used interchangeably with no loss in generality. Likewise the terms sensing pulse and sensing waveform will be used interchangeably with no loss in generality.

The motion of the induced shear wave, including its perturbations associated with reflections and/or reverberation, can be estimated by considering differences in the echoes from sensing waveforms. Said differences are not intended to imply a subtraction; instead, the term is used in a broad sense to mean any aspects that are not identical between waveforms. The term sensing waveform is used here to indicate an ultrasound waveform which has a magnitude too small to generate significant acoustic radiation force, and therefore too small to induce significant shear waves, but large enough to return ultrasound echoes for difference analysis. In an alternative embodiment the same waveform(s) may be used for forcing and sensing.

Shear modulus can be related to other measures of mechanical properties such as Young's modulus and Lamé constants. Thus while the present disclosure is focused on the measurement of shear modulus, these estimates can be transformed to provide estimates of other mechanical properties.

In one embodiment, a device is provided for emitting a forcing waveform and a plurality of sensing waveforms into a sample within a resonant chamber, and processing echoes returned from the sensing waveforms. The device discerns mechanical characteristics of the sample from the characteristics of the resonance. For example, such a device may include at least a controller, a transmitter, a transducer assembly, a resonant chamber, and a data processor.

FIG. 1 shows a high-level block diagram of an example embodiment of the present disclosure. In this embodiment, a controller 110 is provided to manage the timing of various aspects of the testing procedure. For example, controller 110 may control the timing of transmission data digitization and data processing. Controller 110 may be a general purpose computer. In other embodiments, controller 110 may be a dedicated controller such as, for example, a field programmable gate array (FPGA). In one particular embodiment, a Xilinx Spartan 6 FPGA may be utilized. Alternatively an embedded processor or DSP chip might be used as the controller.

Controller 110 may control the timing of, among other things, a transmitter 112. Transmitter 112 may be used to transmit voltage waveforms. The controller 110 can direct the transmitter 112 to connect and disconnect from power supplies to the transducer assembly at specific time intervals. In one example embodiment, the transmitter 112 may transmit desired waveforms comprising positive, negative, and/or neutral voltages with specific time intervals for the transmissions between voltage levels. In other embodiments the transmitter may be capable of multiple voltage amplitude levels, enabling generation of a broader range of waveform shapes. In one embodiment, transmitter 112 comprises a Supertex MD1810 level shifter to control Supertex TC6320 MOSFETs to switch +/−100V supplies, and Supertex TC2320 to clamp transmit waveforms to ground. Various hardware devices, firmware, and/or software, or combinations thereof, may be used as well. The input signals to transmitter 112 may originate from the controller 110.

In the embodiment of FIG. 1, transmitter 112 transmits voltage waveforms to the transducer assembly 114. In this embodiment, the transducer assembly 114 is an ultrasound transducer. The transducer assembly 114 can convert transmit voltage waveforms into ultrasound waveforms and converts ultrasound echoes into receive voltage waveforms. In one example embodiment, the ultrasound transducer is a single-element composite piston transducer. However, other types of ultrasound transducers may be used, and may include hardware, firmware, and/or software, or combinations thereof. In alternative implementations the transducer may include a piezoelectric material (including single crystal materials), CMUT (capacitive micro-machined ultrasound transducer), a relaxor-ferroelectric transducer, thermo-acoustic sources, or voice coils, other transducer technologies. In another alternative embodiment ultrasound transmission is performed using a thermoacoustic method where rapid heating causes thermal expansion that in turn generates an ultrasound wave. Transducer assembly 114 may further comprise an active transducer element—for example, a piezoelectric material—mounted to a single acoustic matching layer, which in turn may be mounted to a polymeric standoff. In one embodiment the transducer is air backed to enhance electro-mechanical efficiency. In one example embodiment, the transducer element of transducer assembly 114 is broad bandwidth and has sensitivity between approximately 5 and 12 MHz. In some embodiments, a series matching inductor is placed between the transducer assembly 114 and the transmitter 112 to tune out an imaginary component of the electrical impedance of the transducer. Other electrical circuits may prove advantageous for matching the electrical impedance of the transmitter, transducer, and receiver.

In an example embodiment, a test sample is placed within chamber 116 for testing. The transducer assembly 114 directs ultrasound energy through the test sample held within chamber 116. In some embodiments, chamber 116 is axisymmetric and has a major axis that is co-linear with the propagation vector of the ultrasound beam. In an alternative embodiment only a portion of the chamber 116 is axisymmetric, while other portions have arbitrary geometries as needed to support filling of the sample and to avoid blocking the ultrasound beam.

In some embodiments chamber 116 is fabricated from a material that is substantially stiffer (higher in shear modulus) than the material being characterized. Thus a resonant chamber fabricated of polystyrene, or similarly rigid material, can be effectively considered infinitely rigid for the purposes of analyzing blood clots. For example, chamber 116 may be several thousand, several hundred thousand, or even over a million times stiffer than the sample within chamber 116. For example, blood clots typically have a shear modulus of a few kiloPascals. Thermoplastics such as polystyrene have a shear modulus of roughly one gigaPascal.

In some embodiments, a “forcing waveform” and a “sensing waveform” may be directed into the test sample in chamber 116. The forcing waveform may be an ultrasonic waveform capable of inducing shear waves in the sample via acoustic radiation force, while the sensing waveform may be a lower energy waveform used to sense various aspects of the sample at a given point in time. These waveforms and their uses are described in more detail below. The modulus of the sample can be estimated by analyzing the resonance of ultrasound-induced shear waves within chamber 116.

In some embodiments the data processor 118 incorporates a number of functions to enable analysis of received echoes. For example, data processor 118 may incorporate a receiver and a digitizer that together provide digital data to a general purpose processor for data analysis. In this embodiment the receiver of the data processor 118 receives and amplifies electrical signals corresponding to ultrasound echoes within chamber 116. The receiver would be operatively coupled to the transducer in this embodiment. The receiver may also include a protection circuit that prevents high-voltage waveforms from overwhelming the one or more amplifiers of the receiver. An example of such a protection circuit is the Supertex MD0100. In some embodiments the input to the protection circuit is tied to the transducer while the output of the protection circuit is coupled to a low-noise amplifier and then to a variable-gain amplifier. Filtering stages also may be interposed to eliminate out-of-band noise. For example, in one embodiment the Analog Devices AD8334 LNA/VGA combination is utilized to amplify the incoming signal.

In one embodiment the receiver may be operatively coupled to a digitizer. Specifically the output of the amplifier might form the input to a digitizer. The digitizer transforms an analog signal to a digital signal. In one example embodiment, a 12-bit analog-digital converter (ADC) such as the Analog Devices AD9238 is utilized.

In the example embodiment of FIG. 1, receive echo data may be stored in a memory within data processor 118. Such a memory could capture digital output from the digitizer. Data processor 118 may include an FPGA, a general purpose processor, a dedicated DSP processor, or some combination of these items. For example, data processor 118 may include an FPGA storage unit in which echo data is temporarily buffered before it is transferred to an embedded processor. In that case the data is again buffered in the embedded processor before being transferred to an embedded PC for processing and modulus estimation. In an example embodiment, the data processor 118 estimates modulus through two distinct and interconnected steps. First, data processor 118 analyzes incoming echo signals to determine the displacement between the echoes returned from various sensing waveforms. In the second step the data processor 118 compares features of the measured displacements to the predicted features of either analytical or computer model predictions for the given chamber 116 geometry to estimate the modulus of the sample within the chamber 116.

FIG. 2 depicts an example of an estimation process, including a data acquisition step 202, a motion estimation step 204, and a modulus estimation step 206. FIG. 2 also shows a plot 224 of modulus over time, including data point 222, the estimation of which is represented in more detail in panels 202, 204, and 206. The data acquisition panel 202 shows the data acquisition process. A series of ultrasound transmission are transmitted into the test chamber. The series of ultrasound waveforms are: (1) sensing pulse 210, (2) forcing pulse 211, (3) sensing pulse 212, (4) sensing pulse 213, (5) sensing pulse 214, and (6) sensing pulse 215. These waveforms are merely examples and the disclosure is not limited to the particular number or order of waveforms shown.

In some example embodiments, the sensing pulses are designed to impart minimal acoustic radiation force on the sample while returning echoes with high signal-to-noise ratio and bandwidth. Sensing pulse 210 may be used to establish a baseline echo for the sample before application of a forcing waveform. The forcing pulse 211, on the other hand, is designed to impart a substantial acoustic radiation force field. Following forcing pulse 211, a series of low intensity sensing pulses (212-215) are transmitted into the sample. The timing between the various sensing pulses is controlled to maintain precision in downstream signal processing steps.

In an alternate embodiment all of the waveforms are of sufficient energy to impart acoustic radiation force. In this particular embodiment these forcing waveforms also act effectively as sensing waveforms. The received echoes from any or all of these waveforms may be processed to estimate modulus using the methods and devices of the present disclosure. The present disclosure is further anticipated to incorporate arbitrary combinations of forcing waveforms, sensing waveforms, and combined use forcing/sensing waveforms.

Undisturbed coagulation of a healthy blood sample leads to the formation of a stiff clot. If however that same sample is mechanically stressed during coagulation then the forming fibrin mesh may be damaged, leading to a soft clot. Any measurement of mechanical properties during coagulation that applies significant mechanical thus has the potential to corrupt the evolution of the mechanical properties being measured. This biasing effect of viscoelastic clot measurement is particularly likely to corrupt the measurement of soft clots formed from the blood of patients with dysfunctional hemostasis. This problem has been previously addressed by the adaptive force measurement method disclosed in Patent Application PCT/US2010/049342. This application discloses a coagulation measurement system wherein the magnitude of the applied force is adjusted to constrain the magnitude of the induced displacement. In times where the clot is soft, the applied force is reduced to avoid damaging the clot. In times where the clot is stiff, the applied force is increased to maximize the sensitivity of the mechanical property estimation. This adaptive force approach is equally applicable to the present disclosure and combination with it is anticipated.

In one example embodiment, sensing waveforms are transmitted at an interval of approximately 122 microseconds, providing a sampling frequency of about 8.2 kHz. Other interrogation frequencies may be used as well. Generally speaking, the interrogation of stiffer materials requires higher interrogation frequency because the shear wave resonance has a high frequency (assuming a constant resonant chamber geometry). More accurate results may be achieved in softer materials by using a lower interrogation frequency so as to minimize the accumulated radiation force from the sensing waveforms and so as to enable the acquisition of data over a longer period of time for a fixed data memory size. With respect to blood samples, suitable interrogation frequencies range from about 2 kHz to about 16 kHz, for example.

A grouping of a single forcing waveform and a plurality of sensing waveforms may be referred to as an “ensemble.” In one embodiment, an ensemble comprises about 500 sensing waveforms. In other embodiments, however, the ensemble may comprise between about 16 and about 2048 sensing waveforms. Other ensemble sizes may be used for measuring materials having higher or lower stiffnesses. One ensemble is processed to yield a single modulus estimate.

In some embodiments, the acquisition time for a single ensemble is approximately 62 milliseconds. However, the acquisition time may be lower or higher. For example, accurate results may be obtained by using an acquisition time of approximately 20-30 milliseconds for a single ensemble. Even lower acquisition times, such as 10 milliseconds, for example, may be used as well. A longer ensemble time period enables accurate measurement of a broader range of moduli. In some embodiments, ensembles may be repeated at a rate of approximately 16 Hz to measure rapidly changing moduli. In other embodiments the physical process under examination (coagulation) is slow enough that ensembles may be repeated at a rate of only once every six seconds and still provide data that accurately reflects the changes in modulus.

In some embodiments it may be advantageous to limit the range of mechanical properties being considered in a given ensemble based upon previously measured mechanical properties for the same test chamber. During coagulation, for example, it is expected that the shear modulus will vary smoothly with time, assuming that measurements are made at a small enough time interval. If, for example, the modulus is 1.0 kPa in a given measurement, then it may be advantageous to limit the range of possible moduli in successive measurements to a range between 0.5 kPa and 2.0 kPa. Even in cases where the modulus range is not explicitly limited, it may be advantageous to smooth modulus estimates across time by either linear filtering (convolution with a filter kernel) or non-linear filtering methods such as median filtering or combinations of both.

Each of the transmitted waveforms travels from left to right along an hourglass beam shape as shown by the series of vertical lines within the test chambers of panel 202 of FIG. 2. The sample is at rest when interrogated by sensing waveform 210. The sample remains at rest when impacted by the forcing waveform 211 because of the inertia of the sample material. Immediately after the forcing waveform passes however, the acoustic radiation force imparted by the forcing waveform 211 causes the sample material along the beam to move in the direction of propagation of the forcing waveform 211. In the embodiment of FIG. 2, this motion is first visible when the sample is interrogated by sensing waveform 212. This movement is shown in the figure as a shaded region lying under the ultrasound beam.

As sensing waveform 213 is transmitted, the radiation-force-induced displacement is beginning to propagate outward from the acoustic beam towards the test chamber wall. This displacement propagates predominantly in the form of a shear wave. Over time, the shear wave reflects from the wall, travels back through the acoustic beam, and reflects from the wall again. This repeated reflection represents resonance of the shear wave within the test chamber. Reverberation eventually settles down as the viscous losses within the sample and other losses damp out the propagating shear wave. Note that for some combinations of modulus, viscosity, and resonant chamber velocity the induced shear wave may reach the chamber wall so quickly that an observer would not see it propagate. Rather it would seem that the entire contents chamber were oscillating in time. Although qualitatively different than a propagating shear wave, this standing wave pattern is nonetheless an example of resonance and is anticipated by the current disclosure. In yet other combinations of modulus, viscosity, and chamber geometry, the excited shear wave might be quite discrete in time and space, and the shear wave reflected from the wall is a waveform that is discrete in space and time. Measurement of modulus by examining the arrival times of such a distinct pulse is one possible embodiment of the present disclosure. In this embodiment the shear wave echo reflection arrival time is an aspect of the differences in the received waveforms that is compared to the same aspect (shear wave reflection arrival time) of the modeled differences.

In yet other implementations, the chamber 116 is so small relative to the shear wavelength that no true shear wave can be generated. In this case it can be said that resonance does not exist for this combination of modulus and chamber geometry. It is still true however that the induced displacement is impacted by the presence of the chamber wall. This boundary effect will act to alter the induced displacement, compared to what would be expected if the same force were applied to an infinite or semi-infinite medium. This variation may take the form of a change in the time dependent displacement from that predicted for a semi-infinite medium. This variation is an aspect of the estimated time-displacement that is caused at least in part by a boundary effect of the chamber wall. In such a case the time course of the induced displacement can be considered in combination with the chamber geometry to estimate the sample modulus. In one embodiment this modulus estimation is performed by comparing the measured displacements to those predicted by a series of computer models, as described in more detail below.

The presence or absence of resonance may be determined by analyzing the differences in received echoes in general or by analyzing estimated time-displacements specifically. In one embodiment the time-displacement curve is analyzed for the presence of a trough (negative peak displacement). If such trough is found then it may be concluded that resonance is occurring. Should no trough (negative peak) in time displacement be observed then it may be concluded that resonance is not present. This conclusion about the presence or absence of resonance may be represented as a parameter indicating the strength of resonance. It this simple example the parameter would hold a value of 1 when a time-displacement trough is detected and a value of 0 when no such trough is detected. This concept can be extended further by considering whether a displacement peak is detected in conjunction with a displacement trough, as this would indicate even stronger resonance. An alternative parameter indicating the strength of resonance would be the mean value of the time-displacement curve. When resonance is absent, the time-displacement curve will be primarily unipolar and will therefore have a high mean value. Alternatively, when resonance is strong, the time-displacement curve will exhibit a strong oscillation around zero and will therefore have a low mean value. An alternative parameter indicating strength of resonance is the ratio of the mean displacement to the peak displacement over the measurement ensemble. A high value of this parameter indicates weak resonance. Other parameters can be computed to indicate the strength of resonance.

Each of the sensing waveforms returns an echo from inhomogeneities (acoustic scatterers) within the sample. In the case of whole blood, these inhomogeneities are primarily red blood cells. The present disclosure can also be used to measure homogenous materials, such as blood plasma, by the addition of polystyrene microspheres, or other agents, which act as acoustic scatterers.

As acoustic scatterers move away from the ultrasound transducer, the acoustic path length between the transducer and the scatterers lengthens. Assuming that the speed of sound is constant, this causes the echoes to arrive at a later time when the targets are pushed farther away by the resonating shear waves. Likewise, the echoes will arrive earlier in cases where the scatterers have moved closer to the ultrasound transducer. These changes in echo arrival time are differences between these waveforms that are indicative of underlying movement of the sample. If the speed of ultrasound propagation (speed of sound c) is known in the sample, or can be measured, then the measured time delays can be related to the underlying physical displacement through the well known relation dx=c dt/2.

In the preceding expression dt is the measured time shift between echoes, c is the speed of sound (ultrasound, not shear waves) and dx is the estimated relative displacement. Note however that the present disclosure does not require knowledge of the speed of sound as the displacement characteristics used to estimate modulus do not necessarily include absolute displacements. In one embodiment of the present disclosure, phase shifts are measured between the various sensing waveform echoes. These phase shifts are differences in the received waveforms that result from the underlying movement of the sample. For a resonating sample these phase shifts will exhibit the feature of oscillation with a frequency related to the modulus of the sample. The observed feature (oscillation frequency) can be compared to a predicted feature (oscillation frequency predicted by theory) to estimate the modulus of the sample.

It is well known that the speed of sound in blood changes as it clots. This evolving speed of sound has little or no effect upon the current measurement, however, for at least two reasons. First, as stated above, many of the algorithms used to relate measured displacement to modulus do not need to know the true displacement; relative displacements are all that are needed. Second, the speed of sound change occurs over a matter of minutes, while the measurement ensembles described herein occur over a matter of milliseconds. Thus the slowly evolving speed of sound has an imperceptible impact on the time-delay estimates from any single ensemble.

The echo data from the ensemble is processed to find differences in the received waveforms that are indicative of motion of targets along the ultrasound beam. One process of analyzing these waveform differences is known as “motion estimation,” a process shown conceptually in the motion estimation panel 204 of FIG. 2. Each echo produced by, for example, sensing pulses 212-215, is compared to the reference echo produced by sensing pulse 210 in order to find the time delay between them. By using a measured or assumed speed of sound, the time delays between various echoes can be converted to displacements. All of the displacement estimates from a single ensemble are combined to form a time-displacement curve 220, as shown vertically on the right of motion estimation panel 204. This time-displacement curve is a feature indicative of the modulus of the sample. Note that the example time-displacement curve depicts oscillations associated with shear-wave resonance and a decay of those oscillations associated with the chamber geometry and the intrinsic viscous damping of the medium.

The motion estimation algorithm used to calculate differences in the receive waveforms may be an algorithm known in the art. Example algorithms include those presented by Kasai (C. Kasai, K. Natnekawa, A. Koyano, and R. Omoto, “Real-Time Two-Dimensional Blood Flow Imaging Using an Autocorrelation Technique” IEEE Trans. Sonics Ultras., vol. SU-32, pp. 458-464, 1985), Loupas (Loupas et al., “Experimental evaluation of velocity and power estimation for ultrasound blood flow imaging, by means of a two-dimensional autocorrelation approach”, IEEE Trans Ultrason Ferroelect Freq Contr. 42:689-699, 1995), and Walker (U.S. Pate. No. 8,306,293).

Alternatively, waveform differences may be analyzed to estimate motion by finding the time-delay corresponding to the peak of the correlation function between different various received echo signals. Direct measurement of time delays can be extended to the measurement of the delay envelope of the waveform after demodulation. As another alternative, relative phase shifts between various received echo waveforms are differences representative of motion in the sample. These phase shifts can be computed digitally by comparing the complex Hilbert Transform of the waveforms associated with different transmissions. As yet another alternative, the receive waveforms can be sampled digitally at an interval approximately ¼ of a period apart so as to approximate In-phase and Quadrature (IQ) signals. This so-called Direct Sampled In-phase and Quadrature (DSIQ) sampling scheme has previously been used to simplify ultrasound beamformer design (US20070016022 A1) and can be applied to calculate waveform differences representative of motion in the present disclosure. In another embodiment the received ultrasound waveforms are processed via quadrature demodulation to yield a complex waveform with an angle between its real and imaginary components indicating the phase of the received signal. This phase is a difference indicative of motion of the sample.

FIG. 3 shows experimental time-displacement curves obtained by the present disclosure. Waveform differences were analyzed to yield displacement estimates over a series of 300 ensembles. Each ensemble consisted of a single forcing waveform and 512 sensing waveforms transmitted at a pulse repetition frequency of 8,206 Hz. Larger displacements and lower frequency oscillations are seen earlier in the clotting process. The first time-displacement curve 302 shows no oscillation, while oscillations become clearer as the clot formed higher modulus. A simple algorithm based upon only the frequency of mechanical resonance would fail with the data of curve 302. This curve corresponds to a modulus for which this chamber geometry cannot support resonance. The present disclosure, however, can estimate the modulus even when no clear oscillation is visible, such as for time-displacement curve 302. Modulus in this case can be estimated by comparing the measured displacements (features) to computer or analytical models of the dynamic sample motion induced by the forcing waveform (predicted features).

The modulus estimation panel 206 of FIG. 2 depicts an example process of estimating the shear modulus of the sample from the experimentally determined time-displacement curve 220. The shape of time-displacement curve 220 is a feature of the estimated displacements that can, in turn, be described as waveform differences. In one example embodiment, a computer model has been used to generate a set of reference models, where each reference model is the predicted time-displacement curve for a given shear modulus and viscosity for the specific test-chamber geometry. These reference models incorporate predicted features associated with the modulus of the computer model. Such a computer model may make use of the Finite Difference Time Domain method, as described below. Alternatively finite element or boundary element computational models might be used. The processor searches through a library of reference models (predicted features) to find the one that most closely matches the experimentally measured time-displacement curve 220 (feature). In this embodiment the library of reference models can be formed offline using a Finite Difference Time-Domain (FDTD) model, as discussed below. Alternatively the reference models may be computed using a finite element or boundary element model. The reference models (predicted features) are shown in modulus estimation panel 206 as gray curves that each correspond to a different shear modulus. Each reference model of modulus estimation panel 206 shows the time-displacement curve 220 (features) overlaid in black. In this particular example embodiment, the reference model for 3 kPa shear modulus most closely matches the time-displacement curve 220. The modulus and viscosity used to form this reference model are the estimate of the modulus and viscosity of the sample. In other embodiments, the computer model calculates a range of models corresponding to potential modulus and viscosities for a given time-displacement curve immediately after the time-displacement curve is generated. However, this dynamic modeling approach may be excessively computationally burdensome relative to searching through predefined reference models. Improvements in computer power or reductions in algorithm complexity will enable dynamic computation of reference models. Such an approach would allow more precise estimation of modulus and/or viscosity. Such an approach could also be coupled with a coarsely sampled (in the modulus dimension) library of reference models to trade-off between computational complexity and storage requirements.

Note that the term “analytical model” may refer to a something as simple as an expression that relates resonant frequency to modulus, or as complex as a full time-displacement waveform predicted by an analytical expression. In cases where a full time-displacement waveform is used as a model it may be advantageous to evaluate that analytical expression and build a set of reference models, analogous to the approach described above for computational models.

In an alternative embodiment, no explicit reference curves are used to estimate modulus. Rather a feature of the time displacement curve, such as its period of oscillation, is computed and used with the period of oscillation predicted by an analytical model (predicted feature) to estimate the modulus. An example analytical model is derived below. That model shows that the resonant frequency (inverse of period) is related to the resonant chamber radius and material modulus and density through the following expression:

$f=\frac{2.4048}{2\pi R}\frac{\sqrt{G}}{\sqrt{\rho }}$

This expression can be rearranged so that the modulus can be estimated directly from the measured period:

$G={\rho \left(\frac{2\pi R}{2.4048T}\right)}^2$

where T is the period of resonance. Thus the period of oscillation is a feature that can be compared to a predicated feature (the period of oscillation of an analytical model) to estimate modulus. Likewise, the period of oscillation is an aspect of induced movement that is caused at least in part by resonance within the chamber. Note that in the case of an analytical model it is not necessary to test a series of model predictions, rather the step of comparing predicted and experimental features can be accomplished through a simple mathematical expression. Such an approach is anticipated by the current disclosure. Alternatively resonant frequency could be used with the above expression to estimate modulus.

In some applications it may not be necessary to estimate a mechanical property. Rather, it may prove useful to instead characterize a mechanical property. In this context we intend the concept “estimating a mechanical property” to refer to making a quantified estimate of a well-known mechanical property such as shear modulus. Alternatively we consider the concept of “characterizing a mechanical property” to refer to determining something about the mechanical property, but not necessarily putting a unit on it or even determining it in a directly proportional way. For example, estimating the resonant frequency of a sample characterizes a mechanical property of that sample without taking the extra step of referencing to the true modulus. Rather than producing a graph of modulus versus time, it might be useful to produce a graph of “resonant frequency” versus time. While “resonant frequency” is not identical to modulus, nor even proportional, tracking resonant frequency would allow one to obtain significant information about coagulation. We thus recognize the utility of such characterization without the need to connect it to a fundamental mechanical property such as shear modulus.

The period of resonance is a feature that can be used as a basis for estimating mechanical properties or characterizing mechanical properties. Given that period is simply the inverse of frequency, estimation of period is equivalent to estimation of frequency, and vice versa. The period of resonance can be estimated in a variety of ways from differences between received waveforms. For the purpose of this discussion we will restrict our consideration of differences between received waveforms to motion estimates made from received waveforms, although other approaches are envisioned particularly including phase changes estimated from the received waveforms. We begin by considering a time-displacement curve, like one of those depicted in FIG. 3. The period of oscillation of such a curve may be estimated by taking a Fast Fourier Transform (FFT) and estimating the frequency at which the energy is at its peak. Alternatively one could subtract the mean from the time-displacement curve and then utilize Burg's method for autoregressive power spectral density estimation. The frequency of peak energy corresponds to the frequency of oscillation. The ROOT-MUSIC algorithm can also be applied to estimate the frequency of oscillation. Other spectral estimation techniques could similarly be employed.

The frequency of oscillation may also be estimated through other means. Specific algorithms are known for estimating the frequency of a decaying sinusoid. One such algorithm is described in the paper “Parameter estimation of exponentially damped sinusoids using a higher order correlation-based approach,” by D .P Ruiz et. al. and published in IEEE Transactions on Signal Processing, vo. 43, no. 11, November 1995. A review of algorithms dedicated to this problem is presented in the paper “Frequency and damping estimation methods—an overview” by T. P. Zielinski and K. Duda and published in Metrology and Measurement Systems, vol. 18, no. 4, 2011. Other approaches are also available.

The period of oscillation may be estimated directly from the time-displacement signal. In one approach, the time of the first trough in displacement (maximum negative displacement) is used as an estimate of ½ of the oscillation period. Since the time-displacement signal is discretely sampled in time but the period can take on continuous values, it is advantageous to apply interpolation. In one embodiment the time-displacement signal is directly interpolated to a higher sampling frequency prior to locating the time of the trough. Such interpolation may be performed by resampling via FFT, piecewise cubic spline interpolation, or other known methods. Alternatively the discrete sample at which the trough occurs can be identified and the true location of the trough found via an analytical interpolation scheme. In one embodiment a parabola is fit to the discrete trough value and its two nearest neighbors, and the time of minimum of that parabola is used as the estimated time of the trough. Higher order functions including piecewise cubic splines can alternatively be used for such interpolation.

The period of displacement oscillation may be estimated through other means as well. In an alternative approach the locations of the first trough and the second peak are determined and the time interval between them is the estimate of ½ of the oscillation frequency. Note that the use of the second peak is generally preferred as the timing of the first peak is distorted by the application of the forcing pulse. This peak-finding strategy can be extended further by estimating the locations of multiple positive and negative peaks and combining those locations to estimate the period. In one embodiment the first and second troughs, and the second and third peaks are identified. Consider the time of the first trough as t_n1, the time of the second peak as t_p2, the time of the second trough as t_n2, and the time of the third peak as t_p3. In this case the period can be related to each of the peaks and troughs as follows:t_n1=T/2+e₁t_p2=T+e₂t_n2=3T/2+e₃t_p3=2T+e₄

where e₁, e₂, e₃, and e₄ represent error terms resulting from noise in the peak/trough location estimates. The period of oscillation may be measured directly from the estimated times of these peaks and troughs through the following expression:{circumflex over (T)}=(2t_n1+t_p2+2/3t_n2+1/2t_p3)/4

This expression has the advantage of reducing the weighting of later peaks and troughs that will have lower relative amplitude and therefore will be more susceptible to noise. Alternatively the period may be estimated from the same peak and trough times through the following expression:{circumflex over (T)}=(t_n2−t_n1+t_p3−t_p2)/2

The present disclosure anticipates other variants of this approach including more or fewer peaks and troughs and alternative expressions for estimating the period from the times of the peaks and troughs.

The period of oscillation may alternatively be estimated from locations of the zero-crossings of the time-displacement curves. While the peak to trough time interval is ½ period, the zero-crossing interval is approximately half as long at ¼ period. The above described methods of combining peak and trough times to estimate oscillation period may be readily modified to combine zero-crossing times to estimate oscillation period. Numerous algorithms for estimating the time of a zero crossing are known. One approach fits a line to the data points near the zero crossing and the finds the time at which the fit line equals zero. Higher order approaches using polynomials or splines are also envisioned.

In some embodiments, the shear-modulus estimate for each ensemble is plotted as a single point of a curve, as shown in FIG. 2. In this particular embodiment, the example calculations represented by panels 202, 204, and 206 result in data point 222. Repeated calculations, for a plurality of ensembles over a period of time, result in the shear-modulus curve 224.

In one example embodiment, the similarity between reference models and the experimentally determined time-displacement curve 220 is quantified using a normalized correlation coefficient. A normalized correlation coefficient between two different signals a[n] and b[n] is given by the following equation:

${\rho }_{\mathrm{ab}}=\frac{\sum _{i=1}^N\left(a\left[i\right]-\stackrel{\_}{a}\right)\left(b\left[i\right]-\stackrel{\_}{b}\right)}{\sqrt{\sum _{i=1}^N{\left(a\left[i\right]-\stackrel{\_}{a}\right)}^2\sum _{i=1}^N{\left(b\left[i\right]-\stackrel{\_}{b}\right)}^2}}$ $\mathrm{where}\stackrel{\_}{a}=\frac{1}{N}\sum _{i=1}^{N}a\left[i\right]\mathrm{and}\stackrel{\_}{b}=\frac{1}{N}\sum _{i=1}^{N}b\left[i\right]$

FIG. 4 depicts the results of modulus estimation by the present disclosure as applied to clotting human blood. A small amount of kaolin was added to a sample of human whole blood. The top panel shows the modulus estimated over time, with the reagent (kaolin) added to the blood sample at time zero. The middle panel depicts the estimated viscosity over time. The lowest panel depicts the normalized correlation between the experimental time-displacement curves and those predicted by the Finite Difference Time Domain model described below. The model used for this figure considered moduli between 10 and 10,000 Pa and viscosities between 0.025 and 0.8 Pa s. A total of 16,032 combinations of modulus and viscosity were modeled. All models assumed a density of 1.06 g/cm³. The geometry of our experimental test chamber was modeled on a spatial grid sampled at 100 um in each dimension. The resonant portion of the test chamber consisted of a cylindrical region 4.2 mm in diameter and approximately 1.5 mm in length with a hemispherical cap the same radius as the cylinder. The first 18 time-modulus estimates were removed from the experimental data and the models prior to computing the correlation between them. A best-fit line was removed from both the experimental time-displacements and computer model time-displacement predictions prior to computing the correlation. The erroneous modulus estimates in the early stages of clot formation occur because the model fails to allow for the low modulus of liquid blood or blood in the very earliest stages of clot formation.

The results shown in FIG. 4 are notable, in part, because of the relatively high correlation between the experimentally determined time-displacement curves and the predictions of the FDTD model. We see from the bottom panel of FIG. 4 that for a well formed clot the correlation between the FDTD model and the experimentally measured time-displacement curves is well over 0.95. This provides a strong indication that the computation model is an appropriate one for these experimental conditions. The correlation between the model and the experimental data is much worse for the liquid blood phase. This is not surprising given that the softest model in that case corresponds to a modulus of 10 Pa, which is far from the liquid state, which has a modulus of zero. Similar correlation tests were performed with time-displacement waveforms formed by numerical evaluation of the analytical model derived below. Although not shown, these results also correlated well with experiment; confirming the utility of the analytical model. The correlations for the analytical model were somewhat worse that the FDTD model. This is likely because the analytical model assumes an infinitely long cylinder, while the FDTD assumes a capped cylinder like the real experimental test chamber.

The correlation between an experimental time-displacement curve and a range of reference models can be plotted as a two-dimensional function of modulus and viscosity, as shown in the contour plot of FIG. 5. This figure depicts the normalized correlation between the experimental time-displacement curves and those predicted by the Finite Difference Time Domain computational model described below. The model used for this figure considered moduli between 10 and 10,000 Pa and viscosities between 0.025 and 0.8 Pa s. A best-fit line was removed from both the experimental time-displacements and computer model time-displacement predictions prior to computing the correlation. The peak correlation lies within the smallest ellipse labeled 0.95 in this figure. The location of the peak corresponds to the estimated modulus and viscosity for a particular example ensemble. The smallest ellipse (highest correlation) is oblong, spread across a large range of viscosities and smaller range of moduli. Small amounts of noise may cause successive estimates to wander around the peak correlation, although it will likely stay confined to the high correlation contour. The shape of the 0.95 correlation contour indicates that the estimates of viscosity may include significant variability, as evidenced by the breadth of the contour in the viscosity dimension. The variability of modulus estimates will be relatively smaller as a result of the narrower contour in this dimension. However, the correlation function is not separable in the viscosity and modulus dimensions. Rather, an error that shifts the viscosity must also shift the modulus. This observation suggests that, in cases where the viscosity can be assumed to hold a certain value, we can fix that parameter and obtain more precise and repeatable estimates of shear modulus. This approach of fixing the viscosity has the added bonus of limiting what was a two-dimensional search (modulus and viscosity) to a one-dimensional search (modulus), which improves computational efficiency.

FIG. 6 depicts the result of modulus estimation with the viscosity held constant. The top panel shows the modulus estimates. The middle panel shows the viscosity, which was fixed at a value of 0.25 Pa s. The bottom panel shows the correlation between the model and the experimental time-displacements. The model used for these estimates was in many ways identical to that used to obtain the results of FIG. 4. There were two notable differences. First, the modulus varied between 0.01 and 10,000 Pa. Second, the viscosity was held constant at a value of 0.25 Pa s. A total of 1,167 models were created. The results of FIG. 6 show a significant improvement in modulus estimates before a firm clot has formed. Although the correlations for these early estimates are fairly low, they are significantly higher than they were for FIG. 4. This result is encouraging as it shows that adjusting the range of the model to incorporate lower moduli yields a closer estimate of mechanical properties when the modulus of the sample is lower. This is another confirmation of the validity of the computational model. The results of FIG. 6 required less than 1/10^th the computation of the results in FIG. 4, with no loss in the quality of the result.

The computational and analytical models described in the present disclosure have a number of adjustable parameters that might be estimated or alternatively held constant, depending upon the specific sample and chamber geometry. In the preceding example it was show to be advantageous to hold the sample viscosity and density constant and evaluate the model for a range of shear moduli. In human blood it is reasonable to hold density at a constant value as the variation among individuals is quite low. While viscosity can vary somewhat more, the advantage of limiting that degree of freedom in estimation may outweigh any absolute error in viscosity. For different samples it may be advantageous to hold modulus constant while allowing either viscosity or density to vary. In yet other embodiments the chamber geometry may not be known precisely and models may allow variation of certain aspects of the chamber geometry among models.

In some embodiments, the value of the normalized correlation coefficient can be used to reject modulus and viscosity estimates that are unreliable and are unlikely to be correct. For example, the normalized correlation coefficient can be used to reject modulus and viscosity estimates where a bubble in the acoustic path corrupts the underlying ultrasound echoes. Processing this corrupted echo data will yield a time-displacement curve that is also excessively noisy. While a typical correlation between the best reference model and the experimental time-displacement curve may be 0.98, the correlation for the corrupted time-displacement data might fall to 0.40, for example. This significant reduction in correlation would clearly indicate that such a measurement is not reliable and should be rejected. This process of rejecting estimates with low correlation may be referred to as “masking.” The appropriate threshold value can be determined empirically from experimental trials, however a reasonable threshold may be around 0.9.

The masking of shear-modulus estimates according to their peak correlation with the reference models is particularly valuable for cleaning up noisy estimates formed before the blood has begun to coagulate. For liquid blood, the time-displacement estimates show very large displacements that are often irregular. The irregular nature of these curves makes it possible for them to fit almost any reference model, although this fit is merely coincidence rather than an actual match. Moreover, the reference model may not actually be applicable until coagulation begins forming a solid “clot.” For example a model based upon a viscoelastic solid is inappropriate for a liquid sample. In these cases, while there is a fit, the quality of match and therefore correlation coefficient is quite low. By simply removing modulus estimates for which the peak correlation is below a reasonable threshold, it is straightforward to remove erroneous modulus estimates associated with liquid blood.

The concept of “correlation masking” was applied to the experimental results of FIG. 6, with the results shown in FIG. 7. The top panel shows the raw modulus estimates; including two data points at the beginning of the measurement. These data points are clearly in error. The central panel shows the correlation between the experimental data and the best-fit models. We select a correlation threshold of 0.6, which was selected to remove clearly erroneous data points. The lowest panel of FIG. 7 shows the remaining modulus estimates after the removal of the moduli with correlation below 0.6. The two erroneous modulus estimates have been eliminated.

The correlation between a model and experimental time-displacement estimates may also be used for other purposes. In one alternative embodiment, multiple different models are formed with variations in both mechanical properties and in underlying assumptions. For example, one family of models could be developed in which the sample is assumed to be a viscoelastic solid, while a second family of models assumes a liquid sample. Within each family the mechanical properties are allowed to vary. In this embodiment the data processor would test the correlation between each model in both of the two families against the estimated time-displacement curve. The specific model with the best correlation to the time-displacement curve would indicate both the value of the mechanical properties and the type of material that best describes the sample. In this specific example the material type would be either a viscoelastic solid or a liquid.

The library of reference models necessarily contains a finite number of reference models. In some embodiments, the modulus estimate may fall between reference models. This limitation can be overcome with little computational cost by interpolating the correlation coefficient around the measured peak to locate the shear modulus corresponding to the best correlation, whether a reference model was computed for that modulus or not. In the simplest case, where the viscosity is held constant and the modulus alone is allowed to vary, the modulus estimate can be determined by locating the peak of a parabola fit through the correlation of the best fit reference model and the correlations of its two nearest neighbors. Other interpolation schemes, including spline based and higher order polynomials, may yield even more precise results. Further, the interpolation described here can be applied to two-dimensional (modulus—viscosity) estimates by using a two-dimensional interpolation scheme. The experimental results presented here all made use of interpolation to form modulus estimates.

The experimentally determined time-displacement curve may be corrupted by a variety of physical and electronic effects. For a finite-volume test chamber, a forcing waveform may generate significant reverberant ultrasound echoes that may extend into the acquisition periods of the following sensing waveforms. In one example embodiment, the impact of such reverberation may be mitigated by designing the frequencies of the forcing waveform and the sensing waveform to lie within different frequency bands, and then using analog or digital filtering so suppress the reverberation associated with the forcing waveform. In an alternate embodiment, sensing waveforms sent shortly after a forcing waveform can be removed from the time-displacement estimation curve prior to searching the reference models for a best fit. This approach of rejecting erroneous time-displacement estimates can either be performed statically (e.g., the first 6 estimates from every time-displacement curve are assumed erroneous) or performed dynamically (e.g., a quality metric is computed to for the time-displacement estimates and estimates with a quality metric below a certain threshold are discarded).

FIG. 16 depicts an experimentally determined time-modulus curve using an analytical expression to relate the measured period of displacement oscillation to the shear modulus. This approach used equation 58 to relate resonant frequency to modulus. The period of resonance was estimated from the time difference between the first trough in displacement to the second peak of displacement. This approach yields a robust estimate of clot modulus, although it is challenged by the early liquid phase before coagulation begins. As described elsewhere in this disclosure, such erroneous estimates could be masked depending upon a parameter indicative of the strength of resonance. This analytical estimate of modulus is quantitatively different from the estimated based on the computational model. This is likely due to the different geometric assumptions of the two approaches. The analytical model assumes an infinitely long cylinder while the computational model assumes a capped cylinder like the experimental test chamber.

The present disclosure anticipates estimating mechanical properties based upon at least three different but interrelated modes of mechanical behavior. In the first mode a forcing pulse is applied to the chamber and has the potential to excite a mechanical resonance of the sample, which is bound by the test chamber. In one embodiment, the frequency of resonance is the basis for estimating the shear modulus of the sample. Note however that for some combinations of modulus, viscosity, density, and chamber geometry, it may not be possible to induce true resonance. For example, if the sample viscosity is high relative to the modulus then the system will be under-damped and no resonance will occur. As the clot stiffens, however, the modulus will increase while the viscosity remains relatively constant. Thus as the clot stiffens, an example device disclosed herein will be able to detect the initiation of resonance. Such a transition might offer a useful and simple indicator of clot formation.

In another embodiment, the viscosity of the sample is low enough relative to the modulus that a well-defined (in space and time) shear wave can be initiated within the sample. In cases where the shear wave is constrained in time and space its reflection from the chamber walls will be clear and easily detected from the differences between echoes received from scatterers within the sample. Knowledge of the chamber geometry (shear wave path length) enables direct estimation of shear wave velocity (and therefore modulus) from the arrival time of shear wave echoes from the chamber wall. As with the prior example of resonance, there are certain combinations of mechanical properties and geometries for which clearly measurable shear waves cannot be generated. As with the resonance example, a transition between states where clearly detectable shear waves cannot be generated to a state where they can be might offer a good proxy for the transition from liquid blood to a formed clot. The timing of such a transition might therefore present a useful measure of clot formation time.

In a third embodiment, neither clear resonance nor propagating shear waves are readily measured. That does not, however, preclude estimation of mechanical properties within the test chamber. For any sample within a finite chamber, the walls of the chamber will change the apparent mechanical impedance of the sample. In this context mechanical impedance describes the relationship between and applied force and the resultant displacement. A very soft sample within a very large test chamber will exhibit a mechanical impedance much like that sample would in an infinitely large test chamber. Alternatively, a stiff sample within a small test chamber will exhibit a mechanical impedance quite different from that observed in a larger test chamber. In the context of dynamic force application, the magnitude of the sample viscosity will further impact the perceived mechanical impedance. The present disclosure anticipates comparing models of the sample within test chamber to experimental measurements to estimate the mechanical properties of the sample. In this embodiment the impact of the wall (boundary effect) is explicitly considered. Since the impact of the boundary (wall) for a given chamber geometry varies depending upon the mechanical properties of the sample, a transformation from the boundary effect being insignificant to it being significant can be a proxy for the transition from a liquid blood sample to a formed clot.

In another embodiment, shown in FIG. 8, the transducer assembly 114 of FIG. 1 incorporates a focusing element 802. The addition of a focusing element increases the intensity of the applied forcing wave and increases the strength of the received echoes. In addition the use of the focusing element 802 enables use of a larger transducer element within transducer assembly 114. The substitution of such a device is wholly contemplated herein by Applicant.

In an example embodiment, shown in FIG. 9, a separable test chamber 116 is provided with a couplant 902. In many commercial applications, the separable test chamber will be particularly useful as a consumable component. Throughout this application we refer to the separable test chamber as “consumable,” while recognizing that it could be a reusable component. In the context of this embodiment, we describe an acoustic path comprising of a transducer assembly 114 (instrument), a couplant 902 (consumable), and test chamber 116 (consumable). Additionally, a focusing element can be interposed between the transducer assembly and couplant, or between the couplant and test chamber. In one alternative embodiment (shown in FIG. 10) the focusing element 1004 is combined with the instrument, and a couplant 1002 is associated with either the instrument or the consumable. The acoustic path from the transducer assembly outward is: transducer assembly 114 (instrument), focusing element 1004 (instrument), couplant 1002 (instrument or consumable), and test chamber 116 (consumable). In yet another alternative embodiment the focusing element is incorporated with the consumable, as shown in FIG. 11. In this embodiment the couplant 1102 can be associated with either the instrument or the consumable, while the focusing element 1104 is associated with the consumable.

In the example embodiment of FIG. 11, the ultrasound waveform travels from couplant 1102 to focusing element 1104. Focusing element 1104 may be shaped such that it refracts incoming ultrasound energy and focuses it within the resonant test chamber as desired. In some embodiments focusing unit 1104 is a thermoplastic, although other suitable materials may be used instead of, or in addition to, a thermoplastic. In one possible embodiment the couplant 1102 may comprise a liquid such as water. In another possible embodiment the couplant 1102 may comprise a water-based or oil-based gel such as that commonly used for ultrasound imaging.

Focusing element 1104 focuses ultrasound energy into chamber 116. Chamber 116 may also be described interchangeably as a resonant chamber, resonance chamber, or testing chamber. The focusing element 1104 is designed so that its curvature and speed of sound act to refract the transmitted ultrasound waves into a focused beam. The focusing element 1104 may be designed to form a sharp or broadly focused beam. A sharp focus will tend to impart a higher magnitude of radiation force, however the sharp focus may make the overall system more sensitive to minor errors in aligning the ultrasound beam and the resonant chamber 116. Alternative the focusing element 1104 may be designed to more broadly focus the ultrasound beam. This may generate a lower magnitude of radiation force, but will make the system less sensitive to errors in the alignment of the ultrasound beam with the resonant chamber 116.

FIG. 12 depicts an alternative embodiment comprising a durable instrument and a consumable test cartridge. The controller 110, transmitter 112, transducer assembly 114, data processor 118, transducer alignment feature 1210, and clamping mechanism 1214 are incorporated in the instrument. The test chamber 116, couplant 1202, focusing element 1204, chamber alignment feature 1212, thermal control unit 1220, and surface modification 1230 are incorporated in the consumable.

In the embodiment of FIG. 12 the consumable component is brought into acoustic contact with the instrument by the action of a clamping mechanism 1214. In one embodiment the clamping mechanism is lead screw actuated by a stepper motor. Other embodiments may include, for example, a mechanically actuated lever or some other mechanical mechanism. In cases where the couplant 1202 is a fairly firm material, the clamping mechanism 1214 might require significant force to eliminate any air gaps and form good acoustic contact. In one embodiment the clamping mechanism 1214 applies about 110 Newtons of force for clamping a single resonant chamber 116. In an alternative embodiment, where a single consumable consists of multiple resonant chambers, a single clamping mechanism may be shared across channels. In this case the clamping mechanism may be required to apply more than 110 Newtons of force. For example, in an embodiment having four chambers, the clamping mechanism may apply about 440 Newtons of force.

In the embodiment of FIG. 12 alignment features 1210 and 1212 are associated with the transducer assembly 114 and the focusing component 1204. The alignment features 1210 and 1212 fit together so that the acoustic beam is precisely generated and accurately placed within the test chamber 116. In one embodiment the alignment feature 1210 incorporates four different components. Two of these are pins that each have an axis parallel with the ultrasound beam. One of these pins may fit into a hole as a component of transducer alignment feature 1212. The other pin may fit a slot as a component of chamber alignment feature 1212. In addition to the pin/hole and pin/slot alignment features, another embodiment may include at least one hard stops to control the degree of compression of the couplant 1202. In one embodiment the hard stops consist of rails above and below the couplant 1202 and associated with the focusing component 1204 that come into contact with the transducer assembly 114. Together these rails may force the focusing assembly 1204 to be parallel with the face of the transducer assembly 114.

In the embodiment of FIG. 12, the test chamber 116 is in thermal contact with a thermal control unit 1220. Since coagulation is a temperature-dependent process, the use of a thermal control unit serves to improve repeatability and increase the rate of coagulation so that clinicians receive results more rapidly. In one embodiment the temperature control unit is an aluminum block with an embedded thermistor and an electric heater. The thermistor acts as an input to a PID Controller (Proportional-Integral-Derivative Controller). In another embodiment the thermal control unit makes use of a Peltier device to enable both heating and cooling of the sample relative to its incoming temperature.

In the embodiment of FIG. 12, the interior surface of the test chamber 116 is formed or modified so as to have a surface 1230 conducive to clot adhesion. Such a treatment may be needed so that as the blood clots and platelets contract, the clot cannot pull away from the chamber wall. This may be advantageous because, should the clot pull away from chamber wall, the assumed properties upon which a computational or analytic model is built can be wrong. For example, instead of a uniform material extending all throughout the test chamber 116, a chamber with a retracted clot will include a serum layer filling a gap between the contracted clot and the resonant chamber wall. Such a gap will alter the resonant geometry of the resonant chamber, thereby biasing modulus measurements under these conditions. Even if clot retraction only opens a tiny gap between the clot and the chamber wall, the boundary conditions upon which our models are predicated will no longer be appropriate and therefore estimation of mechanical properties will be corrupted. Even when we restrict ourselves to characterizing the mechanical properties (not computing absolute properties) the characterization will be corrupted by the altered boundary conditions in the chamber. These errors can be eliminated by forming a chamber with a rough inner surface to which the clot can stay firmly adhered. Such a surface can be formed by patterning the injection mold used to form the chamber, or by roughening the interior surface after manufacture through a mechanism like bead blasting. Alternatively the interior surface can undergo plasma or corona treatment, which can both roughen the surface and modify its chemical properties to enhance clot adhesion. In another alternative embodiment the inner surface can be treated by adhering a protein like fibronectin, to which the clot can adhere.

Another alternative embodiment of the present disclosure is shown in FIG. 13. In this embodiment the transducer assembly 114 comprises a focusing element 802, and two separate transducer elements 1302 and 1304. The first transducer element 1302 is connected to the transmitter but not the data processor and is intended to transmit forcing waveforms. The second transducer element 1304 is connected to both the transmitter and the data processor and is intended to transmit sensing waveforms and receive their echoes. The configuration of FIG. 13 may prove advantageous because, for example, the ultrasound beam shape most desirable for applying forcing waveforms may be different than that which is optimal for sensing. Furthermore, as the forcing waveforms are generally higher in energy than the sensing waveforms, this configuration likely reduces the requirements on protection circuitry associated with the data processor and improves sensitivity.

Another alternative embodiment of the present disclosure is shown in FIG. 14. In this embodiment the transducer assembly 114 comprises a focusing element 802, and two separate transducer elements 1402 and 1404. The first transducer element 1402 is connected to the transmitter but not the data processor and is intended solely for waveform transmission. The second transducer element 1404 is connected to the data processor and is intended solely for echo reception. The configuration of FIG. 14 may prove advantageous because it will eliminate the need for protection circuitry associated with the data processor.

One advantage of the present disclosure over the prior art is a significant reduction in required blood volume for testing. For example, at least one prior-art method required a 60 ml sample of blood per test. In various embodiments described herein, only about 330 microliters per test chamber is required. Larger or smaller volumes of blood may be used as well, depending on the testing conditions. This reduced test chamber volume allows a patient to provide less blood while also benefiting from multiple tests in parallel within a single cartridge. In some embodiments, a single cartridge may comprise multiple test chambers that can be tested in parallel. In alternative embodiments in which the focusing element is incorporated in the transducer assembly the total blood volume may be reduced even further to as few as tens of microliters per test chamber.

In one embodiment, described above, the sensing waveform is designed so as to impart no significant acoustic radiation force upon the sample. In other embodiments, however, the sensing waveform may impart a small amount of radiation force. In these embodiments the measured displacements will include a superposition of resonant time-displacement curves. The dominant component of this superposition will be the resonant time-displacement curve resulting from the application of the forcing waveform. However, added to this curve will be smaller delayed responses from each of the sensing waveforms. In cases where the displacements due to the applied sensing waveforms have a significant effect, the ideal response can be reconstructed by deconvolving the inputs. The impact of small amounts of radiation force from the sensing waveforms might also be mitigated by lowering the rate of sensing waveform transmission.

The present disclosure has numerous applications beyond the measurement of coagulation in human blood. For example, in the veterinary sciences and in biomedical research it may be useful to quantify the process of blood clotting samples from non-human animals.

In the food industry there is strong interest in quantifying the mechanical properties of foods and food additive. In that domain it would be straightforward to fill a resonant chamber with the food to be tested. One might find it valuable to apply the present disclosure to monitor the aging of cheese or other dynamic processes. In one example, the present disclosure could be applied to measure aging of a cheese wheel. Large cheese wheels can be over 50 inches in diameter. Such a cheese wheel could be considered a resonant chamber in the context of the present disclosure. With dimensions approximately 300 times larger than the resonant test chamber considered in the detailed examples above, it may prove advantageous to scale the operating frequency similarly. This suggests that measurement of such a large object could require the use of lower frequency acoustic energy; possibly extending into the audible regime. The chamber of this application would be formed by the boundary of the cheese wheel itself with the surrounding air forming a free boundary, in contrast to the fixed boundary discussed above for blood.

In the art of tissue engineering there is an unmet need for measuring the mechanical changes that occur as cells grow, mature, and build extracellular matrices within certain tissue-engineered constructs. The present disclosure could be applied to that problem by growing cells within a resonant chamber. Alternatively the matrix upon which cells grow could be placed in a nutritive broth so the geometry of the matrix itself forms the resonant geometry. In this case the analytical or computational model would assume a free boundary, rather than the infinitely stiff boundary considered elsewhere in this application.

Example—Finite Difference Time Domain Model

The example model described below is a derivation of a Difference Time Domain model for radiation force induced shear waves in cylindrical geometry. This formulation can be used to model time-displacements in the resonant chamber.

Our analysis begins with a velocity-stress formulation of the shear wave equation. Note that while this formulation considers velocity, results computed with this approach can be numerically integrated to yield displacements, so as to be consistent with experimental results. We expand upon traditional formulations by including a term to account for viscous losses.

$\begin{array}{cc}\rho \frac{\partial }{\partial t}\overrightarrow{v}\left(\overrightarrow{x},t\right)=\nabla ·\sigma \left(\overrightarrow{x},t\right)+\overrightarrow{F}\left(\overrightarrow{x},t\right)& \left(1\right)\end{array}$ $\begin{array}{cc}\frac{\partial }{\partial t}\sigma \left(\overrightarrow{x},t\right)=\left(\mu +\eta \frac{\partial }{\partial t}\right)\left[\left(\nabla \overrightarrow{v}\left(\overrightarrow{x},t\right)\right)+{\left(\nabla \overrightarrow{v}\left(\overrightarrow{x},t\right)\right)}^T\right]& \left(2\right)\end{array}$

We expand equation 1 by expressing it in cylindrical coordinates and expanding the vector velocity into its constituent components.

$\begin{array}{cc}\frac{\partial {\sigma }_{\mathrm{rr}}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{r\theta }}{\partial \theta }+\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial z}+\frac{1}{r}\left({\sigma }_{\mathrm{rr}}-{\sigma }_{\mathrm{\theta \theta }}\right)+F_r=\rho \frac{\partial v_r}{\partial t}& \left(3\right)\end{array}$ $\begin{array}{cc}\frac{\partial {\sigma }_{r\theta }}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\theta \theta }}{\partial \theta }+\frac{\partial {\sigma }_{\mathrm{\theta z}}}{\partial z}+\frac{2}{r}{\sigma }_{r\theta }+F_{\theta }=\rho \frac{\partial v_{\theta }}{\partial t}& \left(4\right)\end{array}$ $\begin{array}{cc}\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\mathrm{\theta z}}}{\partial \theta }+\frac{\partial {\sigma }_{\mathrm{zz}}}{\partial z}+\frac{1}{r}{\sigma }_{\mathrm{rz}}+F_z=\rho \frac{\partial v_z}{\partial t}& \left(5\right)\end{array}$

Further consideration of our problem allows significant simplification. First, we recognize that the only body force is the applied ultrasound radiation force. Assuming that force is entirely in the z direction, we can set F_θ=F_r=0. Assuming that our test chamber and applied radiation force are entirely axisymmetric, we can set all dependencies upon θ equal to zero. Applying these simplifications to equations 3-5 yields:

$\begin{array}{cc}\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial z}=\rho \frac{{\partial }^2{u}_r}{\partial t^2}& \left(6\right)\end{array}$ $\begin{array}{cc}0=0& \left(7\right)\end{array}$ $\begin{array}{cc}\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial r}+\frac{1}{r}{\sigma }_{\mathrm{rz}}+F_z=\rho \frac{{\partial }^2{u}_z}{\partial t^2}& \left(8\right)\end{array}$

We follow a similar strategy to expand equation 2, to yield:

$\begin{array}{cc}\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial t}=\left(\mu +\eta \frac{\partial }{\partial t}\right)\left(\frac{\partial v_r}{\partial z}+\frac{\partial v_z}{\partial r}\right)& \left(9\right)\end{array}$

Collecting equations 6-9 we have:

$\begin{array}{cc}\rho \frac{\partial v_r}{\partial t}=\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial z}& \left(10\right)\end{array}$ $\begin{array}{cc}\rho \frac{\partial v_z}{\partial t}=\frac{\partial {\sigma }_{\mathrm{rz}}}{\partial r}+\frac{1}{r}{\sigma }_{\mathrm{rz}}+F_z& \left(11\right)\end{array}$ $\begin{array}{cc}\frac{\partial {\sigma }_{\mathrm{rz}}}{\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)& \left(12\right)\end{array}$

Equations 10-12 form a system of partial differential equations that can be solved together to predict how radiation force will induce shear waves and how those induced shear waves will interact. This system of equations is particularly amenable to a finite difference solution using a staggered grid approach, similar to the Yee method. The finite difference representation of equations 10-12 are:

$\begin{array}{cc}{v}_{r_{j,k}}^{i+1}=v_{r_{j,k}}^i+\frac{\Delta t}{\rho \Delta z}\left({\sigma }_{{\mathrm{rz}}_{j,k+\frac{1}{2}}}^i-{\sigma }_{{\mathrm{rz}}_{j,k-\frac{1}{2}}}^i\right)& \left(13\right)\end{array}$ $\begin{array}{cc}{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 }_{{\mathrm{rz}}_{j+1,k+\frac{1}{2}}}^i-{\sigma }_{{\mathrm{rz}}_{j,k+\frac{1}{2}}}^i\right)+\frac{1}{2\left(j+\frac{1}{2}\right)\Delta r}\left({\sigma }_{{\mathrm{rz}}_{j,k+\frac{1}{2}}}^i-{\sigma }_{{\mathrm{rz}}_{j+1,k+\frac{1}{2}}}^i\right)+F_{z_{j++\frac{1}{2},k+\frac{1}{2}}}^i\right)& \left(14\right)\end{array}$ $\begin{array}{cc}{\sigma }_{{\mathrm{rz}}_{j,k+\frac{1}{2}}}^{i+1}={\sigma }_{{\mathrm{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)& \left(15\right)\end{array}$

The mathematical formulation described above can be implemented numerically using a staggered grid approach as shown in FIG. 15. Each computational cycle consists of two steps. In the first step the velocity components are computed using the finite difference equations above. In the next step the shear is computed using the finite difference equations.

Example—Analytical Modeling of Impulse-Induced Resonance

The example analytical model derived below represents the mechanical model of blood clot displacement in a cylindrical resonant chamber. The model qualitatively captures the observed behaviors in real-world experiments.

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

$\begin{array}{cc}\left(\left(\lambda +\mu +\left({\eta }_p+{\eta }_s\right)\frac{d}{\mathrm{dt}}·\right)▽\left(▽·\overrightarrow{u}\left(\overrightarrow{x},t\right)\right)\right)+\left(\left(\mu +{\eta }_s\frac{d}{\mathrm{dt}}·\right){▽}^2\overrightarrow{u}\left(\overrightarrow{x},t\right)\right)-\rho \frac{d^2}{{\mathrm{dt}}^2}\overrightarrow{u}\left(\overrightarrow{x},t\right)=\overrightarrow{F}\left(\overrightarrow{x},t\right)& \left(16\right)\end{array}$

In our initial experiments we use a cylindrical test volume with radiation force applied to a smaller radius cylindrical volume with the same central axis as the test volume. This geometry ensures that there will be no variation in parameters with angle. We also assume an infinite length, which further simplifies the problem by eliminating variation in range. We apply these simplifications to the components of (16):

$\begin{array}{cc}▽\left(▽·\overrightarrow{u}\left(\overrightarrow{x},t\right)\right)=▽\left(\frac{1}{r}\frac{d}{\mathrm{dr}}\left({\mathrm{ru}}_r\right)+\frac{1}{r}\frac{d}{d\varphi }{u}_{\varphi }+\frac{d}{\mathrm{dz}}{u}_z\right)& \left(17\right)\end{array}$

Because of cylindrical symmetry and an infinite length we can assume that all derivatives with respect to ϕ and z are equal to zero. Thus (17) becomes:

$\begin{array}{cc}▽\left(▽·\overrightarrow{u}\left(\overrightarrow{x},t\right)\right)=▽\left(\frac{1}{r}\frac{d}{\mathrm{dr}}\left({\mathrm{ru}}_r\right)\right)& \left(18\right)\end{array}$

We now take the gradient to yield:

$\begin{array}{cc}▽\left(▽·\overrightarrow{u}\left(\overrightarrow{x},t\right)\right)=\left[\frac{d}{\mathrm{dr}}\left(\frac{1}{r}\frac{d}{\mathrm{dr}}\left({\mathrm{ru}}_r\right)\right)+\frac{d}{\mathrm{dr}}\left(\frac{1}{r}\frac{d}{d\varphi }{u}_{\varphi }\right)+\frac{d^2}{\mathrm{drdz}}\left(u_z\right)\right]\hat{r}+& \left(19\right)\end{array}$ $\left[\frac{1}{r}\frac{d^2}{d\varphi \mathrm{dr}}\left({\mathrm{ru}}_r\right)+\frac{1}{r}\frac{d^2}{d{\varphi }^2}{u}_{\varphi }+\frac{1}{r}\frac{d^2}{d\varphi \mathrm{dz}}\left(u_z\right)\right]\hat{\varphi }+$ $\left[\frac{d}{\mathrm{dz}}\left(\frac{1}{r}\frac{d}{\mathrm{dr}}\left({\mathrm{ru}}_r\right)\right)+\frac{d}{\mathrm{dz}}\left(\frac{1}{r}\frac{d}{d\varphi }{u}_{\varphi }\right)+\frac{d^2}{\mathrm{drdz}}\left(u_z\right)\right]\hat{z}$

We again appreciate that all derivatives with respect to ϕ and z are equal to zero. Thus (19) can be simplified to:

$\begin{array}{cc}▽\left(▽·\overrightarrow{u}\left(\overrightarrow{x},t\right)\right)=\left[\frac{d}{\mathrm{dr}}\left(\frac{1}{r}\frac{d}{\mathrm{dr}}\left({\mathrm{ru}}_r\right)\right)\right]\hat{r}+0\hat{\varphi }+0\hat{z}& \left(20\right)\end{array}$

We consider another component of (16):

$\begin{array}{cc}{▽}^2\overrightarrow{u}\left(\overrightarrow{x},t\right)=& \left(21\right)\end{array}$ $\left[{▽}^2{u}_r-\frac{1}{r^2}{u}_r-\frac{2}{r^2}\frac{d}{d\varphi }{u}_{\varphi }\right]\hat{r}+\left[{▽}^2{u}_{\varphi }-\frac{1}{{\varphi }^2}{u}_{\varphi }-\frac{2}{r^2}\frac{d}{d\varphi }{u}_r\right]\hat{\varphi }+\left[{▽}^2{u}_z\right]\hat{z}$

Again, we recognize that all derivatives with respect to ϕ and z are equal to zero. Furthermore, the displacement component with respect to ϕ is also equal to zero. Thus (21) simplifies to:

$\begin{array}{cc}{▽}^2\overrightarrow{u}\left(\overrightarrow{x},t\right)=\left[{▽}^2{u}_r-\frac{1}{r^2}{u}_r\right]\hat{r}+0\hat{\varphi }+\left[{▽}^2{u}_z\right]\hat{z}& \left(22\right)\end{array}$

We expand the Laplacians to yield:

$\begin{array}{cc}{▽}^2\overrightarrow{u}\left(\overrightarrow{x},t\right)=\left[\frac{1}{r}\frac{d}{\mathrm{dr}}\left(r\frac{d}{\mathrm{dr}}{u}_r\right)+\frac{1}{r^2}\frac{d^2}{d{\varphi }^2}{u}_r+\frac{d^2}{{\mathrm{dz}}^2}{u}_r+\frac{1}{r^2}{u}_r\right]\hat{r}+0\hat{\varphi }+\left[\frac{1}{r}\frac{d}{\mathrm{dr}}\left(r\frac{d}{\mathrm{dr}}{u}_z\right)+\frac{1}{r^2}\frac{d^2}{d{\varphi }^2}{u}_z+\frac{d^2}{{\mathrm{dz}}^2}{u}_z\right]\hat{z}& \left(23\right)\end{array}$

Again, we recognize that the derivatives with respect to ϕ and z are equal to zero. Thus (23) becomes:

$\begin{array}{cc}{▽}^2\overrightarrow{u}\left(\overrightarrow{x},t\right)=\left[\frac{1}{r}\frac{d}{\mathrm{dr}}\left(r\frac{d}{\mathrm{dr}}{u}_r\right)+\frac{1}{r^2}{u}_r\right]\hat{r}+0\hat{\varphi }+\left[\frac{1}{r}\frac{d}{\mathrm{dr}}\left(r\frac{d}{\mathrm{dr}}{u}_z\right)\right]\hat{z}& \left(24\right)\end{array}$

The partial differential equation of (16) can be split into three different equations; one for each direction. We begin by considering the result in the z direction: Note that the displacement in the z direction is a function only of r and t as the ϕ and z dependencies are eliminated by radial symmetry and the infinite length of the model:

$\begin{array}{cc}\left(\left(\lambda +\mu +\left({\eta }_p+{\eta }_s\right)\frac{d}{\mathrm{dt}}·\right)▽\left(▽·u_z\left(r,t\right)\right)\right)+\left(\left(\mu +{\eta }_s\frac{d}{\mathrm{dt}}·\right){▽}^2{u}_z\left(r,t\right)\right)-\rho \frac{d^2}{{\mathrm{dt}}^2}{u}_z\left(r,t\right)=F_z\left(t\right)& \left(25\right)\end{array}$

Substituting (20) and (24) back into (25) yields:

$\begin{array}{cc}\left(\left(\lambda +\mu +\left({\eta }_p+{\eta }_s\right)\frac{d}{\mathrm{dt}}·\right)0\right)+\left(\left(\mu +{\eta }_s\frac{d}{\mathrm{dt}}·\right)\frac{1}{r}\frac{d}{\mathrm{dr}}\left(r\frac{d}{\mathrm{dr}}{u}_z\left(r,t\right)\right)\right)-\rho \frac{d^2}{{\mathrm{dt}}^2}{u}_z\left(r,t\right)=F_z\left(t\right)& \left(26\right)\end{array}$

Expanding (26) yields:

$\begin{array}{cc}\left(\left(\mu +{\eta }_s\frac{d}{\mathrm{dt}}·\right)\frac{1}{r}\left(r\frac{d^2}{{\mathrm{dr}}^2}{u}_z\left(r,t\right)+\frac{d}{\mathrm{dr}}{u}_z\left(r,t\right)\right)\right)-\rho \frac{d^2}{{\mathrm{dt}}^2}{u}_z\left(r,t\right)=F_z\left(t\right)& \left(27\right)\end{array}$

Simplifying further yield:

$\begin{array}{cc}\left(\left(\mu +{\eta }_s\frac{d}{\mathrm{dt}}·\right)\left(\frac{d^2}{{\mathrm{dr}}^2}{u}_z\left(r,t\right)+\frac{1}{r}\frac{d}{\mathrm{dr}}{u}_z\left(r,t\right)\right)\right)-\rho \frac{d^2}{{\mathrm{dt}}^2}{u}_z\left(r,t\right)=F_z\left(t\right)& \left(28\right)\end{array}$

Note that (28) includes no dependence upon the radial displacement. Thus the potential system of partial differential equations becomes a single partial differential equation.

Oscillatory Force:

We solve equation 28 by considering a solution of the form below, in which the displacement is oscillatory.u_z(r,t)=S(r)e^jωt  (29)

We further assume that the applied force is also oscillatory.F_z(r,t)=F(r)e^jωt  (30)

Substituting (29) into (28) yields:

$\begin{array}{cc}\left(\left(\mu +{\eta }_s\frac{d}{\mathrm{dt}}·\right)\left(\frac{d^2}{{\mathrm{dr}}^2}S\left(r\right)e^{j\omega t}+\frac{1}{r}\frac{d}{\mathrm{dr}}S\left(r\right)e^{j\omega t}\right)\right)-\rho \frac{d^2}{{\mathrm{dt}}^2}S\left(r\right)e^{j\omega t}=F\left(r\right)e^{j\omega t}& \left(31\right)\end{array}$

Evaluating the derivatives with respect to time yields:

$\begin{array}{cc}\left(\left(\mu +j{\mathrm{\omega \eta }}_s\right)\left(\frac{d^2}{{\mathrm{dr}}^2}S\left(r\right)e^{j\omega t}+\frac{1}{r}\frac{d}{\mathrm{dr}}S\left(r\right)e^{j\omega t}\right)\right)+{\mathrm{\rho \omega }}^{2}S\left(r\right)e^{j\omega t}=F\left(r\right)e^{j\omega t}& \left(32\right)\end{array}$

We set aside the e^jωt term, which occurs throughout (32) to yield:

$\begin{array}{cc}\left(\left(\mu +j{\mathrm{\omega \eta }}_s\right)\left(\frac{d^2}{{\mathrm{dr}}^2}S\left(r\right)+\frac{1}{r}\frac{d}{\mathrm{dr}}S\left(r\right)\right)\right)+{\mathrm{\rho \omega }}^{2}S\left(r\right)=F\left(r\right)& \left(33\right)\end{array}$

To simplify notation we replace (μ+jωη_s) with a general complex shear modulus of G so that (33) becomes:

$\begin{array}{cc}G\left(\frac{d^2}{{\mathrm{dr}}^2}S\left(r\right)+\frac{1}{r}\frac{d}{\mathrm{dr}}S\left(r\right)\right)+{\mathrm{\rho \omega }}^{2}S\left(r\right)=F\left(r\right)& \left(34\right)\end{array}$

The solution of this equation depends upon the exact form of the forcing function. One simple form assumes a constant force F within some transmit beam radius r_tx, and a force of 0 outside that radius. For this simple forcing function the solution has two domains, one inside the radius r_tx and one outside that radius. The solution is:

$\begin{array}{cc}\begin{array}{c}{S}_1\left(r\right)=c_1{J}_0\left(\frac{\sqrt{{\mathrm{\rho \omega }}^2}}{\sqrt{G}}r\right)+c_2{Y}_0\left(\frac{\sqrt{{\mathrm{\rho \omega }}^2}}{\sqrt{G}}r\right)+\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\mathrm{for}r\le r_{\mathrm{tx}}\\ S_2\left(r\right)=c_3{J}_0\left(\frac{\sqrt{{\mathrm{\rho \omega }}^2}}{\sqrt{G}}r\right)+c_4{Y}_0\left(\frac{\sqrt{{\mathrm{\rho \omega }}^2}}{\sqrt{G}}r\right)\mathrm{for}r>r_{\mathrm{tx}}\end{array}& \left(35\right)\end{array}$

Where J₀ is the zero^th order Bessel Function of the first kind and Y₀ is the zero^th order Bessel Function of the second kind. We can further simplify this solution by considering our boundary conditions. The derivative of the solution at r=0 must be 0. Since Y₀ does not have a finite derivative at 0 we know that c₂=0. This leaves our solution as:

$\begin{array}{cc}\begin{array}{c}{S}_1\left(r\right)=c_1{J}_0\left(\mathrm{Ar}\right)+\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\mathrm{for}r\le r_{\mathrm{tx}}\\ S_2\left(r\right)=c_3{J}_0\left(\mathrm{Ar}\right)+c_4{Y}_0\left(\mathrm{Ar}\right)\mathrm{for}r>r_{\mathrm{tx}}\end{array}& \left(36\right)\end{array}$

Note that to simplify notation we have replaced

$\frac{\sqrt{{\mathrm{\rho \omega }}^2}}{\sqrt{G}}$ with A. We now solve for the arbitrary constants by considering the boundary conditions and continuity conditions of the problem. We assume that the clot is rigidly adhered to the chamber wall and therefore S₂(R)=0, where R is the radius of the test chamber. Furthermore, the two solutions must be continuous at their junction so that S₁(r_b)=S₂(r_b). Finally, the two solutions must have continuous derivatives at their junction so that S₁′(r_b)=S₂′(r_b). We can express these three conditions as shown below:

$\begin{array}{cc}{c}_3{J}_0\left(\mathrm{AR}\right)+c_4{Y}_0\left(\mathrm{AR}\right)=0& \left(37\right)\end{array}$ $c_1{J}_0\left({\mathrm{Ar}}_b\right)+\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}-c_3{J}_0\left({\mathrm{Ar}}_b\right)-c_4{Y}_0\left({\mathrm{Ar}}_b\right)=0$ $c_1{J}_0^{\prime }\left({\mathrm{Ar}}_b\right)-c_3{J}_0^{\prime }\left({\mathrm{Ar}}_b\right)-c_4{Y}_0^{\prime }\left({\mathrm{Ar}}_b\right)=0$

Recognizing that J₀′(Ar)=−AJ₁(Ar) and Y₀′(Ar)=−AY₁(Ar), and reformulating the above expressions into a single system of linear equations yields:

$\begin{array}{cc}\left[\begin{array}{ccc}0& J_0\left(\mathrm{AR}\right)& Y_0\left(\mathrm{AR}\right)\\ J_0\left({\mathrm{Ar}}_b\right)& -J_0\left({\mathrm{Ar}}_b\right)& -Y_0\left({\mathrm{Ar}}_b\right)\\ -{\mathrm{AJ}}_1\left({\mathrm{Ar}}_b\right)& {\mathrm{AJ}}_1\left({\mathrm{Ar}}_b\right)& {\mathrm{AY}}_1\left({\mathrm{Ar}}_b\right)\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_3\\ c_4\end{array}\right]=\left[\begin{array}{c}0\\ -\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\\ 0\end{array}\right]& \left(38\right)\end{array}$

This system of equations can be solved using Gaussian Elimination. We are primarily concerned within the region of force application, where r≤r_tx. Thus the coefficient of greatest interest is c₁. Solving for c₁ yields:

$\begin{array}{cc}{c}_1=\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}{J}_0\left(\mathrm{AR}\right)}\left(\frac{J_0\left(\mathrm{AR}\right)Y_1\left({\mathrm{Ar}}_b\right)-Y_0\left(\mathrm{AR}\right)J_1\left({\mathrm{Ar}}_b\right)}{Y_0\left({\mathrm{Ar}}_b\right)J_1\left({\mathrm{Ar}}_b\right)-J_0\left({\mathrm{Ar}}_b\right)Y_1\left({\mathrm{Ar}}_b\right)}\right)& \left(39\right)\end{array}$ $\begin{array}{cc}{c}_3=-\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}{J}_0\left(\mathrm{AR}\right)}\left(\frac{Y_0\left(\mathrm{AR}\right)J_1\left({\mathrm{Ar}}_b\right)}{Y_0\left({\mathrm{Ar}}_b\right)J_1\left({\mathrm{Ar}}_b\right)-J_0\left({\mathrm{Ar}}_b\right)Y_1\left({\mathrm{Ar}}_b\right)}\right)& \left(40\right)\end{array}$ $\begin{array}{cc}{c}_4=\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\left(\frac{J_1\left({\mathrm{Ar}}_b\right)}{Y_0\left({\mathrm{Ar}}_b\right)J_1\left({\mathrm{Ar}}_b\right)-J_0\left({\mathrm{Ar}}_b\right)Y_1\left({\mathrm{Ar}}_b\right)}\right)& \left(41\right)\end{array}$

The above expressions offer a rigorous solution to the analytical formulation presented here. Unfortunately these rigorous solutions are prone to numerical instability when the operands of the Bessel Functions become large. In these cases numerical errors in the evaluation of the Bessel Functions are reinforced by subtraction so that the evaluation of c₁, c₃, and c₄ may be effectively useless.

The numerical instabilities of (39), (40), and (41) can be mitigated by using simpler expressions for the Bessel Functions with large operands. We use the following simplifications for large operands.

$\begin{array}{cc}{J}_0\left(x\right)\approx \sqrt{\frac{2}{\pi x}}\cos \left(x-\frac{\pi }{4}\right)\mathrm{for}\mathrm{large}x& \left(42\right)\end{array}$ $\begin{array}{cc}{J}_1\left(x\right)\approx \sqrt{\frac{2}{\pi x}}\cos \left(x-\frac{3\pi }{4}\right)\mathrm{for}\mathrm{large}x& \left(43\right)\end{array}$ $\begin{array}{cc}{Y}_0\left(x\right)\approx \sqrt{\frac{2}{\pi x}}\sin \left(x-\frac{\pi }{4}\right)\mathrm{for}\mathrm{large}x& \left(44\right)\end{array}$ $\begin{array}{cc}{Y}_1\left(x\right)\approx \sqrt{\frac{2}{\pi x}}\sin \left(x-\frac{3\pi }{4}\right)\mathrm{for}\mathrm{large}x& \left(45\right)\end{array}$

By applying these expressions to equations 39-41, and then employing trigonometric identities, we are able to reformulate c₁, c₃, and c₄ so that they are numerically stable. Note that this reformulation is only valid for large operands, so an empirically determined transition between the two expressions must be employed.

The received signal is a weighted sum of the signals received from each of the annuli within the receive beam. Formally the signals may be considered complex exponentials and the resultant sum a complex exponential, which is then analyzed to determine the displacement. For small displacements such that the first terms of the Taylor Series are reasonable approximations to the complex exponential (cos(x)≈1 and (sin(x)≈x), the complex exponential can be ignored so that the estimated displacement will approximate the sum of displacements across the annuli.

Thus the estimated displacement is weighted sum of the actual displacements over an axisymmetric region. If the effective receive beam radius is less than the transmit beam radius then the estimated displacement is:

$\begin{array}{cc}{d}_1=2\pi {\int }_0^{r_{\mathrm{rx}}}r·s_1\left(r\right)\mathrm{dr}=2\pi {\int }_0^{r_{\mathrm{rx}}}r·\left(c_1{J}_0\left(\mathrm{Ar}\right)+\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\right)\mathrm{dr}& \left(46\right)\end{array}$

The integral of (46) is readily computed with the knowledge that

${\int }_0^{x}r\left(J_0\left(\mathrm{Ar}\right)\right)\mathrm{dr}=\frac{x}{A}{J}_1\left(\mathrm{ax}\right).$

$\begin{array}{cc}{d}_1=c_1\frac{2\pi r_{\mathrm{rx}}}{A}{J}_1\left({\mathrm{Ar}}_{\mathrm{rx}}\right)+\frac{2\pi r_{\mathrm{rx}}^{2}F}{2\sqrt{{\mathrm{\rho \omega }}^2}}& \left(47\right)\end{array}$

If however the receive beam is larger than the transmit beam, then the solution takes on the form:

$\begin{array}{cc}{d}_2=2\pi \left({\int }_0^{r_{\mathrm{tx}}}r·s_1\left(r\right)\mathrm{dr}+{\int }_{r_{\mathrm{tx}}}^{r_{\mathrm{rx}}}r·s_2\left(r\right)\mathrm{dr}\right)& \left(48\right)\end{array}$ $\begin{array}{cc}{d}_2=2\pi \left({\int }_0^{r_{\mathrm{tx}}}r·\left(c_1{J}_0\left(\mathrm{Ar}\right)+\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\right)\mathrm{dr}+{\int }_{r_{\mathrm{tx}}}^{r_{\mathrm{rx}}}r·\left(c_3{J}_0\left(\mathrm{Ar}\right)+c_4{Y}_0\left(\mathrm{Ar}\right)\right)\mathrm{dr}\right)& \left(49\right)\end{array}$ $\begin{array}{cc}{d}_2=2\pi \left(c_1\frac{r_{\mathrm{tx}}}{A}{J}_1\left({\mathrm{Ar}}_{\mathrm{tx}}\right)+\frac{{\mathrm{Fr}}_{\mathrm{tx}}^2}{2\sqrt{{\mathrm{\rho \omega }}^2}}+\frac{c_3}{A}\left(r_{\mathrm{rx}}{J}_1\left({\mathrm{Ar}}_{\mathrm{rx}}\right)-r_{\mathrm{tx}}{J}_1\left({\mathrm{Ar}}_{\mathrm{tx}}\right)\right)+\frac{c_4}{A}\left(r_{\mathrm{rx}}{Y}_1\left({\mathrm{Ar}}_{\mathrm{rx}}\right)-r_{\mathrm{tx}}{Y}_1\left({\mathrm{Ar}}_{\mathrm{tx}}\right)\right)\right)& \left(50\right)\end{array}$

Displacement due to Oscillatory Force from a small beam:

In some cases it will be helpful to consider the oscillatory displacement in the limit of a very small ultrasound beam. We will consider equation 36 in the limit as r_tx approaches zero.

$\begin{array}{cc}\underset{r_b\to 0}{\lim }{S}_1\left(r\right)=\underset{r_b\to 0}{\lim }\left(c_1{J}_0\left(\mathrm{Ar}\right)+\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\right)\mathrm{for}r\le r_{\mathrm{tx}}& \left(51\right)\end{array}$ $\begin{array}{cc}{c}_1=\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}{J}_0\left(\mathrm{AR}\right)}\left(\frac{J_0\left(AR\right)Y_1\left(A{r}_b\right)-Y_0\left(\mathrm{AR}\right)J_2\left({\mathrm{Ar}}_b\right)}{Y_0\left({\mathrm{Ar}}_b\right)J_1\left({\mathrm{Ar}}_b\right)-J_0\left({\mathrm{Ar}}_b\right)Y_1\left({\mathrm{Ar}}_b\right)}\right)& \left(52\right)\end{array}$

We recognize that the only term dependent upon r_b is the constant, c₁. Thus we take the limit of c₁ as r_b goes to zero.

$\begin{array}{cc}\underset{r_b\to 0}{\lim }{c}_1=\underset{r_b\to 0}{\lim }\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}{J}_0\left(\mathrm{AR}\right)}\left(\frac{J_0\left(\mathrm{AR}\right)Y_1\left({\mathrm{Ar}}_b\right)-Y_0\left(\mathrm{AR}\right)J_1\left({\mathrm{Ar}}_b\right)}{Y_0\left({\mathrm{Ar}}_b\right)J_1\left({\mathrm{Ar}}_b\right)-J_0\left({\mathrm{Ar}}_b\right)Y_1\left({\mathrm{Ar}}_b\right)}\right)& \left(53\right)\end{array}$ $\begin{array}{cc}\underset{r_b\to 0}{\lim }{c}_1=\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\left(-1+\frac{Y_0\left(\mathrm{AR}\right)}{J_0\left(\mathrm{AR}\right)}\right)& \left(54\right)\end{array}$

Thus the displacement is equal to:

$\begin{array}{cc}\underset{r_b\to 0}{\lim }{S}_1\left(r\right)=\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\left(-1+\frac{Y_0\left(\mathrm{AR}\right)}{J_0\left(\mathrm{AR}\right)}\right)J_0\left(\mathrm{Ar}\right)+\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\mathrm{for}r\le r_{\mathrm{tx}}& \left(55\right)\end{array}$

We simplify further by considering the displacement at the center, i.e. the beam location:

$\begin{array}{cc}\underset{r_b\to 0}{\lim }{S}_1\left(0\right)=\frac{F}{\sqrt{{\mathrm{\rho \omega }}^2}}\frac{Y_0\left(\mathrm{AR}\right)}{J_0\left(\mathrm{AR}\right)}& \left(56\right)\end{array}$ $\mathrm{Where}=\frac{\sqrt{{\mathrm{\rho \omega }}^2}}{\sqrt{G}}$

By examining this expression we recognize that the displacement will be maximum, i.e. the system will be in resonance, when J₀(AR)=0. The first zero occurs at J₀(2.4048). Thus the resonant frequency can be solved as follows:

$\begin{array}{cc}\frac{\sqrt{{\mathrm{\rho \omega }}^2}}{\sqrt{G}}R=2.4048& \left(57\right)\end{array}$ $\begin{array}{cc}f=\frac{2.4048}{2\pi R}\frac{\sqrt{G}}{\sqrt{\rho }}& \left(58\right)\end{array}$

Constant Force:

The above analysis is valid only for oscillatory forces. In the static case, where the force is constant over time, we must perform a separate analysis. We begin by reconsidering equation (28), copied below for clarity

$\begin{array}{cc}\left(\left(\mu +{\eta }_s\frac{d}{\mathrm{dt}}·\right)\left(\frac{d^2}{{\mathrm{dr}}^2}{u}_z\left(r,t\right)+\frac{1}{r}\frac{d}{\mathrm{dr}}{u}_z\left(r,t\right)\right)\right)-\rho \frac{d^2}{{\mathrm{dt}}^2}{u}_z\left(r,t\right)=F_z\left(t\right)& \left(59\right)\end{array}$

Since we are specifically interested in the static, or DC problem, we eliminate all the derivatives with respect to t, as these must be equal to zero. We also change notation to indicate that we are considering the static shear modulus,

$\begin{array}{cc}{G}_s\left(\frac{d^2}{{\mathrm{dr}}^2}{u}_z\left(r\right)+\frac{1}{r}\frac{d}{\mathrm{dr}}{u}_z\left(r\right)\right)=F_z& \left(60\right)\end{array}$

As with the oscillatory problem, the solution of this equation depends upon the exact form of the forcing function. One simple form assumes a constant force F within some transmit beam radius r_tx, and a force of 0 outside that radius. For this simple forcing function the solution has two domains, one inside the radius r_tx, and one outside that radius. The solution is:

$\begin{array}{cc}\begin{array}{cc}{S}_1\left(r\right)=c_1+\frac{{\mathrm{Fr}}^2}{4G_s}& \mathrm{for}r\le r_{\mathrm{tx}}\\ S_2\left(r\right)=c_3+c_4\log \left(r\right)& \mathrm{for}r>r_{\mathrm{tx}}\end{array}& \left(61\right)\end{array}$

We now solve for the arbitrary constants by considering the boundary conditions and continuity conditions of the problem. We assume that the clot is rigidly adhered to the chamber wall and therefore S₂(R)=0, where R is the radius of the test chamber. Furthermore, the two solutions must be continuous at their junction so that S₁(r_b)=S₂(r_b). Finally, the two solutions must have continuous derivatives at their junction so that S₁′(r_b)=S₂′(r_b). Rather than repeat the analysis used for the oscillatory solution, we instead present the solution below:

$\begin{array}{cc}{c}_1=\frac{{\mathrm{Fr}}_{\mathrm{tx}}^2}{2G_s}\left(\log \left(r_{\mathrm{tx}}\right)-\log \left(R\right)-\frac{r_{\mathrm{tx}}}{2}\right)& \left(62\right)\end{array}$ $\begin{array}{cc}{c}_3=-\frac{{\mathrm{Fr}}_{tx}^2}{2G_s}\log \left(R\right)& \left(63\right)\end{array}$ $\begin{array}{cc}{c}_4=\frac{F{r}_{\mathrm{tx}}^2}{2G_s}& \left(64\right)\end{array}$

The estimated displacement is the average of the actual displacement over the receive beam. If the effective receive beam radius is less than the transmit beam radius then the estimated displacement is:

$\begin{array}{cc}{d}_1=2\pi {\int }_0^{r_{\mathrm{rx}}}r·s_1\left(r\right)\mathrm{dr}=2\pi {\int }_0^{r_{\mathrm{rx}}}r·\left(c_1+\frac{{\mathrm{Fr}}^2}{4G_s}\right)\mathrm{dr}& \left(65\right)\end{array}$

The integral of (47) is readily computed.

$\begin{array}{cc}{d}_1=\pi r_{\mathrm{rx}}^2{c}_1+\frac{\pi {\mathrm{Fr}}_{\mathrm{rx}}^4}{8G_s}& \left(66\right)\end{array}$

If however the receive beam is larger than the transmit beam, then the solution takes on the form:

$\begin{array}{cc}{d}_2=2\pi \left({\int }_0^{r_{\mathrm{tx}}}r·s_1\left(r\right)\mathrm{dr}+{\int }_{r_{\mathrm{tx}}}^{r_{\mathrm{rx}}}r·s_2\left(r\right)\mathrm{dr}\right)& \left(67\right)\end{array}$ $\begin{array}{cc}{d}_2=2\pi \left(\frac{c_1{r}_{\mathrm{tx}}^2}{2}+\frac{{\mathrm{Fr}}_{\mathrm{tx}}^4}{16G_s}+\left(\frac{c_3{r}^2}{2}+\frac{c_4{r}^2}{2}\log \left(r\right)-\frac{c_4{r}^2}{4}\right){❘}_{r_{\mathrm{tx}}}^{r_{\mathrm{rx}}}\right)& \left(68\right)\end{array}$ $\begin{array}{cc}{d}_2=2\pi \left(\frac{c_1{r}_{\mathrm{tx}}^2}{2}+\frac{{\mathrm{Fr}}_{\mathrm{tx}}^4}{16G_s}+\frac{c_3}{2}\left(r_{\mathrm{rx}}^2-r_{\mathrm{tx}}^2\right)+\frac{c_4}{2}\left(r_{\mathrm{rx}}^2\log \left(r_{\mathrm{rx}}\right)-r_{\mathrm{tx}}^2\log \left(r_{\mathrm{tx}}\right)\right)-\frac{c_4}{4}\left(r_{\mathrm{rx}}^2-r_{\mathrm{tx}}^2\right)\right)& \left(69\right)\end{array}$

Viscoelastic Tissue Models:

Successful application of the presented model requires the selection of an appropriate viscoelastic model, as implemented through the complex, frequency dependent shear modulus G. We have explored both the Kelvin-Voigt model (Illustration 1) and the Jeffrey model (Illustration 2). The Kelvin-Voigt model is generally well suited to modeling clot behavior however it fails to capture certain poroelastic behaviors that may be observed in blood clots. In this sense the Jeffrey model is sometimes superior.

Kelvin-Voigt Model G=μ+jωη  Illustration 1

Jeffrey Model

$G=-{\mathrm{\omega \eta }}_1\frac{{\mathrm{\omega \eta }}_2-i\mu }{\mu +i\omega \left({\eta }_1+{\eta }_2\right)}$ Illustration 2

While specific embodiments have been described in detail in the foregoing detailed description and illustrated in the accompanying drawings, it will be appreciated by those skilled in the art that various modifications and alternatives to those details could be developed in light of the overall teachings of the disclosure and the broad inventive concepts thereof. It is understood, therefore, that the scope of the present disclosure is not limited to the particular examples and implementations disclosed herein, but is intended to cover modifications within the spirit and scope thereof as defined by the appended claims and any and all equivalents thereof.

Claims

1. A method for characterizing a mechanical property of a test sample, the method comprising: transmitting a plurality of waveforms into a test chamber, including at least one forcing waveform and at least two sensing waveforms;receiving, via at least one transducer, at least two sensing waveforms from the test chamber corresponding to the at least two transmitted sensing waveforms;estimating, by a processor, a resonance of the test sample from the at least two received sensing waveforms; anddetermining, by the processor, a modulus value associated with stiffness of a clot formed of the test sample based on the estimated resonance; andoutputting, by the processor, a hemostatic parameter derived from the determined modulus value, wherein the outputted hemostatic parameter is used to evaluate and/or treat a hemostatic dysfunction.

2. The method of claim 1, wherein the at least one forcing waveform comprises an ultrasonic waveform that carries sufficient energy so that its absorption and reflection within the test sample generates a radiation force.

3. The method of claim 2, wherein the generated radiation force is induced along a direction of propagation of the at least one forcing waveform.

4. The method of claim 2, wherein the generated radiation force induces a shear wave that travels within the test chamber and reflects from one or more walls of the test chamber.

5. The method of claim 4, wherein a single reflected shear wave is detected by the transducer and its arrival time is used as a basis for estimating the modulus of the clot formed of the test sample in the test chamber.

6. The method of claim 1, wherein each of the at two transmitted sensing waveforms has a magnitude sufficient to return the at least two received sensing waveforms but not sufficient to induce shear waves.

7. The method of claim 1, wherein the estimated modulus value comprises a shear modulus value for the modulus of the clot formed of the test sample in the test chamber.

8. The method of claim 1, wherein the estimated modulus value comprises a Young's modulus value or a Lamé constant value for the modulus of the clot formed of the test sample in the test chamber.

9. The method of claim 1, wherein the modulus value associated with the stiffness of the clot formed of the test sample is determined by comparing the estimated resonance to a reference model derived from parameters associated with the at least one forcing waveform and parameters associated with the test chamber.

10. The method of claim 9, wherein the reference model comprises a plurality of curves each having a time-displacement profile associated with a given modulus value.

11. The method of claim 1, wherein the modulus value associated with the stiffness of the clot formed of the test sample is determined by comparing the estimated resonance to an experimentally-derived model associated with the test sample.

12. The method of claim 1, wherein the at least one transmitted forcing waveform is transmitted as one or more acoustic pulses.

13. The method of claim 1, wherein the at least one transmitted forcing waveform is transmitted as an ultrasonic beam.

14. The method of claim 1, wherein the at least two transmitted sensing waveforms are transmitted as two or more acoustic pulses.

15. The method of claim 1, wherein the at least two transmitted sensing waveforms are transmitted as two or more ultrasonic beams.

16. The method of claim 1, further comprising: repeating the transmitting, receiving, estimating, and determining steps to generate a shear-modulus curve, wherein parameters associated with the shear-modulus curve are extracted as the hemostatic parameter.

17. The method of claim 1, further comprising: transmitting a second plurality of waveforms into a second test chamber of a cartridge, including at least one second forcing waveform and at least two second sensing waveforms, wherein the cartridge comprises the test chamber and the second test chamber;receiving, via at least one second transducer, at least two second sensing waveforms from the second test chamber corresponding to the at least two transmitted second sensing waveforms;estimating, by the processor, a resonance of the second test sample from the at least two received second sensing waveforms; anddetermining, by the processor, a second modulus value associated with stiffness of a second clot formed of the second test sample based on the estimated resonance of the second test sample; andoutputting, by the processor, a second hemostatic parameter derived from the determined second modulus value, wherein the outputted second hemostatic parameter is used along with the first hemostatic parameter to evaluate and/or treat the hemostatic dysfunction.

18. The method of claim 17, wherein four shear-modulus curves, including the shear-modulus curve, are generated from characterizations of at least four test chambers of the cartridge, including the test chamber.

19. The method of claim 1, wherein the outputted hemostatic parameter is used to evaluate for platelet damage, factor and fibrinogen consumption, and/or presence of residual anti-coagulant in the test sample.

20. The method of claim 1, wherein the method is used for a point-of-care hemostasis testing.