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 clot, 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.
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.
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.
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
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.
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
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 axi symmetric, 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 odulus) 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
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
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. In 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
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. Namekawaa, 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. Pat. 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.
The modulus estimation panel 206 of
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:
This expression can be rearranged so that the modulus can be estimated directly from
measured period:
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
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 ½ 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 tn1, the time of the second peak as tp2, the time of the second trough as tn2, and the time of the third peak as as tp3. In this case the period can be related to each of the peaks and troughs as follows:
t
n1
=T/2+e1
t
p2
=T+e
2
t
n2=3T/2+e3
t
p3=2T+e4
where e1, e2, e3, and e4 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)}=(2tn1+tp2+⅔tn2+½tp3)/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)}=(tn2−tn1+tp3−tp2)/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
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:
The results shown in
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
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
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).
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
In an example embodiment, shown in
In the example embodiment of
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.
In the embodiment of
In the embodiment of
In the embodiment of
In the embodiment of
Another alternative embodiment of the present disclosure is shown in
Another alternative embodiment of the present disclosure is shown in
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.
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.
We expand equation 1 by expressing it in cylindrical coordinates and expanding the vector velocity into its constituent components.
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θ=Fr=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:
We follow a similar strategy to expand equation 2, to yield:
Collecting equations 6-9 we have:
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:
The mathematical formulation described above can be implemented numerically using a staggered grid approach as shown in
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.
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):
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:
We now take the gradient to yield:
We again appreciate that all derivatives with respect to ϕ and z are equal to zero. Thus (19) can be simplified to:
We consider another component of (16):
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:
We expand the Laplacians to yield:
Again, we recognize that the derivatives with respect to ϕ and z are equal to zero. Thus (23) becomes:
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:
Substituting (20) and (24) back into (25) yields:
Expanding (26) yields:
Simplifying further yields:
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)ejωt (29)
We further assume that the applied force is also oscillatory.
F
z(r,t)=F(r)ejωt (30)
Substituting (29) into (28) yields:
Evaluating the derivatives with respect to time yields:
We set aside the ejωt term, which occurs throughout (32) to yield:
To simplify notation we replace (μ+jωηs) with a general complex shear modulus of G so that (33) becomes:
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 rtx, and a force of 0 outside that radius. For this simple forcing function the solution has two domains, one inside the radius rtx and one outside that radius. The solution is:
Where J0 is the zeroth order Bessel Function of the first kind and Y0 is the zeroth 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 Y0 does not have a finite derivative at 0 we know that c2=0. This leaves our solution as:
Note that to simplify notation we have replaced with
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 S2(R)=0, where R is the radius of the test chamber. Furthermore, the two solutions must be continuous at their junction so that S1(rb)=S2(rb). Finally, the two solutions must have continuous derivatives at their junction so that S1′(rb)=S2′(rb). We can express these three conditions as shown below:
Recognizing that J′0(Ar)=−AJ1(Ar) and Y′0(Ar)=−AY1(Ar), and reformulating the above expressions into a single system of linear equations yields:
This system of equations can be solved using Gaussian Elimination. We are primarily concerned within the region of force application, where r≤rtx. Thus the coefficient of greatest interest is c1. Solving for c1 yields:
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 c1, c3, and c4 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.
By applying these expressions to equations 39-41, and then employing trigonometric identities, we are able to reformulate c1, c3, and c4 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:
The integral of (46) is readily computed with the knowledge that
If however the receive beam is larger than the transmit beam, then the solution takes on the form:
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 rtx approaches zero.
We recognize that the only term dependent upon rb is the constant, c1. Thus we take the limit of c1 as rb goes to zero.
Thus the displacement is equal to:
We simplify further by considering the displacement at the center, i.e. the beam location:
By examining this expression we recognize that the displacement will be maximum, i.e. the system will be in resonance, when J0(AR)=0. The first zero occurs at J0(2.4048). Thus the resonant frequency can be solved as follows:
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
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,
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 rtx, and a force of 0 outside that radius. For this simple forcing function the solution has two domains, one inside the radius rtx and one outside that radius. The solution is:
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 S2(R)=0, where R is the radius of the test chamber. Furthermore, the two solutions must be continuous at their junction so that S1(rb)=S2(rb). Finally, the two solutions must have continuous derivatives at their junction so that S1′(rb)=S2′(rb). Rather than repeat the analysis used for the oscillatory solution, we instead present the solution below:
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:
The integral of (47) is readily computed.
If however the receive beam is larger than the transmit beam, then the solution takes on the form:
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.
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.
Number | Date | Country | |
---|---|---|---|
Parent | 16656492 | Oct 2019 | US |
Child | 17185516 | US | |
Parent | 15641803 | Jul 2017 | US |
Child | 16656492 | US | |
Parent | 14660700 | Mar 2015 | US |
Child | 15641803 | US |