A General Framework for Modeling Sub- and Ultraharmonics of Ultrasound Contrast Agent Signals with MISO Volterra Series

Sub- and ultraharmonics generation by ultrasound contrast agents makes possible sub- and ultraharmonics imaging to enhance the contrast of ultrasound images and overcome the limitations of harmonic imaging. In order to separate different frequency components of ultrasound contrast agents signals, nonlinear models like single-input single-output (SISO) Volterra model are used. One important limitation of this model is its incapacity to model sub- and ultraharmonic components. Many attempts are made to model sub- and ultraharmonics using Volterra model. It led to the design of mutiple-input singe-output (MISO) Volterra model instead of SISO Volterra model. The key idea of MISO modeling was to decompose the input signal of the nonlinear system into periodic subsignals at the subharmonic frequency. In this paper, sub- and ultraharmonics modeling with MISO Volterra model is presented in a general framework that details and explains the required conditions to optimally model sub- and ultraharmonics. A new decomposition of the input signal in periodic orthogonal basis functions is presented. Results of application of different MISO Volterra methods to model simulated ultrasound contrast agents signals show its efficiency in sub- and ultraharmonics imaging.


Introduction
Medical diagnostic using ultrasound imaging was greatly improved with the introduction of ultrasound contrast agents. In ultrasound imaging, contrast agents are microbubbles [1]. Historically, the important difference between the acoustic impedance of the tissue and the gas encapsulated within the microbubbles was the first step to improve the contrast of echographic images. However, the contrast was still improved by taking into account the nonlinear behavior of microbubbles. In fact, when microbubbles were insonified by a sinusoidal excitation, they respond by generating harmonic components [2]. For example, second harmonic imaging (SHI) [3] consists in transmitting a signal at frequency 0 and receiving echoes at twice the transmitted frequency 2 0 . However, harmonic generation during the propagation of ultrasound in the nonperfused tissue limits the contrast [4].
Many years ago, experimental studies have shown the existence of subharmonics at 0 /2 [5] and ultraharmonics at ((3/2) 0 , (5/2) 0 , . . .) [6] in the microbubble response under specific conditions of frequency and pressure. The absence of these components in the backscattered signal by the tissue has enabled the introduction of sub-and ultraharmonics as an alternative of the harmonic imaging in order to enhance the contrast. Sub-and ultraharmonic imaging consists of transmitting a signal of frequency 0 and extracting components at 0 /2, (3/2) 0 , (5/2) 0 , . . ..
Many models have been developed to understand the dynamics of the microbubble [2]. Microbubble oscillation can be accurately described using models such as Rayleigh-Plesset modified equation [7][8][9]. However, to enable optimal separation of harmonic components, other nonlinear models like single-input single-output (SISO) Volterra model have been preferred [10]. A well known limitation of SISO Volterra model is its capacity to model exclusively harmonic components sub-and ultraharmonics are not modeled [11].

Computational and Mathematical Methods in Medicine
To overcome this difficulty, Boaghe and Billings [12] have proposed a multiple-input single-output (MISO) Volterrabased method. Input signals are specified by having subharmonic component at frequency 0 / . This approach has been applied in ultrasound medical imaging [13].
However, neither Boaghe and Billings [12] nor Samakee and Phukpattaranont [13] have clearly justified the required conditions to design a MISO Volterra decomposition able to model sub-and ultraharmonics.
To answer this untreated point, we propose a more general framework which firstly gives a clear justification regarding the choice of the model and secondly can offer interesting alternatives.
This paper is organized as follows: after recalling Volterra model and presenting the general framework of MISO Volterra methods, simulations of contrast ultrasound medical imaging followed by results are presented. Finally, a discussion completed by a conclusion closes the paper.
The vector of kernels h is calculated to minimize the mean square error (MSE) between ( ) and̂( ) according to the equation where E is the symbol of the mathematical expectation. Vector h is calculated using the least squares method if (X X) is invertible. Otherwise, regularization techniques can be used. Nevertheless, as reported in [12], it is not possible to model sub-and ultraharmonics with SISO Volterra model under this formulation. This is due to the fact mentioned in [12] that SISO Volterra model can only model frequencies at integer multiples of the input frequency.
To overcome this limitation, Boaghe and Billings [12] proposed a MISO Volterra-based solution and not any more a SISO Volterra. This point is discussed in Section 3.

General Framework of MISO Volterra Model
According to Boaghe and Billings' claims [12], it is possible to model sub-and ultraharmonic components of the signal ( ) if the excitation signal to Volterra model has the subharmonic component at 0 / . The solution proposed by Boaghe and Billings [12] to show up the sub-harmonic component at frequency 0 / is to decompose the input signal ( ) into multiple input signals ( ), each signal having frequency components at 0 and 0 / . From our point of view, Boaghe and Billings' approach [12] claimed two conditions that are intrinsically coupled by the choice of the decomposition method as follows: (i) the input signal to Volterra model has sub-harmonic frequency component at 0 / ; (ii) Volterra system is a MISO system described by The block diagram of MISO Volterra model is presented in Figure 2.
A third condition that is not really explained in [12], however, it is a crucial condition to carry out this modeling procedure. It is the orthogonality condition between each multiple input of MISO Volterra model. Taking into account this third condition makes it possible to generalize Boaghe and Billings' approach presented in [12] as follows: where are coefficients to be adjust and Ψ ( , 0 , ) is the periodic orthogonal basis functions having a spectral component at 0 / . Different bases could be proposed. In this study, two bases are presented as follows.
(1) In [12] a first periodic basis of orthogonal functions is proposed as follows: where is the sampling period, * represents the convolution product, and Rect 1/ 0 ( ) is the rectangular function equal to 1 when −1/2 0 < < −1/2 0 and equal to zero otherwise. Note that this approach is MISO1. (2) In the present work, a new periodic basis of orthogonal functions is presented as follows: wherẽ( ) = H( ( )) is the Hilbert transform of ( ) and 0 = 2 0 . Note that this second is MISO2 approach. For our application in contrast medical imaging, the subharmonic frequency is 0 /2 [5][6][7], so = 2.

Computational and Mathematical Methods in Medicine
The two signals 1 ( ) and 2 ( ) for the two previous bases are represented in Figure 3.
It is obvious to show that for the two bases, the signals 1 ( ) and 2 ( ) are orthogonal because ∑ 1 ( ) 2 ( ) = 0 (From a statistical point of view, the two signals 1 ( ) and Finally, if the components ( ) are orthogonal to each other, then this also means that the output of Volterra model ( ) can be decomposed as follows: where the componentŝ( ) are also orthogonal to each other. A proof of this propriety is given in Appendix A.
The consequence of this statement is that MISO Volterra model can be considered as parallel SISO Volterra models as depicted in Figure 4.

Simulations
To validate the different proposed bases and to quantify its performances for application in contrast ultrasound medical imaging, realistic simulations are proposed. To carry out the simulations, the free simulation program bubblesim developed by Hoff [7] was used to calculate the oscillations and scattered echoes for a specified contrast agent and excitation pulse. A modified version of Rayleigh-Plesset equation was chosen. The model presented by Church [15] and then modified by Hoff [7] is based on the theoretical description of microbubbles as air-filled particles with surface layers of elastic solids. In order to simulate the mean behavior of a microbubble cloud, we hypothesized that the response of a cloud of microbubbles was times the response of a single microbubble with the mean properties.
The incident burst to the microbubble is a sinusoidal wave of frequency 0 = 4 MHz (The resonance frequency of a microbubble of 1.5 m is about 2.25 MHz. Therefore, the emission frequency at 4 MHz is nearly the double of the resonance frequency.) To ensure the presence of sub-and ultraharmonics with moderate destruction of microbubbles, Forsberg

Results
In this research, the performances of different modeling methods are evaluated qualitatively and quantitatively.

Qualitative Evaluation.
To evaluate qualitatively the two MISO methods, MISO1 (with the basis proposed in [12]) and MISO2 (with the new basis proposed in the present work) with respect to SISO Volterra method, temporal representations of ( ) and̂( ), and spectral representations | ( )| 2 and |̂( )| 2 of the nonlinear system backscattered by the contrast agent in nonlinear mode are presented in Figure 5. Results presented in Figure 5 are obtained for a signal to tissue ratio SNR = ∞ and using Volterra model of order = 3 and memory = 19.
To better distinguish the different harmonic components of ultrasound signal, six cycles of 0.05 s are presented in Figure 5(a), and a bandwidth of 13 MHz covering the 3 harmonics potentially accessible in ultrasound imaging is presented in Figure 5(b). For both types of representations, the fundamental frequency, harmonics, sub-and ultraharmonics are well apparent. In Figure 5(a) (top), only harmonic components at 0 , 2 0 , and 3 0 are modeled by SISO Volterra. This result confirms that SISO Volterra system is unable to correctly model sub-and ultraharmonics at frequencies

Quantitative Evaluation.
To determine accurately the performances of the two methods and to know which Volterra approach provides the best performances a quantitative study is necessary. The relative mean square error (RMSE) defined as follows: is evaluated for different noise levels at the system output. The noise level, adjusted as a function of SNR, is Gaussian and white. Ten realizations are made to evaluate the fluctuations of RMSE. RMSE for SNR = ∞, 20, 15, and 10 dB is reported in Figure 6. A zoom in Figure 6(d) shows the fluctuations of the EQMR around a mean value. The main result of these simulations shows that regardless the SNR values, MISO Volterra methods provide a much better RMSE than SISO Volterra method. In fact, a gap between SISO Volterra method and the two methods MISO1 and MISO2 going from 5 to 16 dB can be obtained depending on the SNR conditions. These results confirm that SISO Volterra method is not suitable for sub-and ultraharmonic modeling. A zoom on Figure 6 a small advantage in favor of MISO2 method with respect to MISO1 method for memory values smaller than 6 is noted.
Finally, the more the memory increases, the more the RMSE decreases, indicating that the different methods tend asymptotically toward the optimal solution.

Discussions and Conclusions
In the present research, we proposed a general framework that describes harmonic, sub-, and ultraharmonics modeling using Volterra decomposition. This framework allowed us to highlight three essential criteria instead of two, to accurately model sub-and ultraharmonics: (i) as suggested in [12], the basis should be periodic of period 0 / ; (ii) as suggested in [12], Volterra system should be a MISO system; (iii) as suggested in this work, the decomposition of the input signal to Volterra model ( ) must be done with an orthogonal basis.
This general framework has also justified the different steps of the decomposition thus allowing to propose new periodic orthogonal bases more efficient. It is the same for the choice of the order of Volterra model, which was limited to three. In fact, for more or less severe constraints on the ultrasound transducers bandwidth, the order can be reduced or increased.
This more general formulation provides a methodological basis for optimal sub-and ultraharmonics contrast imaging and opens a new research axis for more efficient periodic orthogonal bases of MISO Volterra systems and also for new MISO systems based on Hammerstein models or Wiener models.

A. Decomposition of MISO Volterra Model of 2 Inputs to 2 SISO Volterra Models
A MISO Volterra model with inputs is equivalent to SISO Volterra models if and only if the mean square error to be minimized between ( ) and̂( ) is the same in both cases. We will determine the conditions that must be satisfied by the inputs 1 ( ) and 2 ( ) of MISO Volterra when = 2, to have this equivalence.
For the 2 SISO Volterra models of inputs 1 ( ) and 2 ( ) and outputs 1 ( ) and 2 ( ), respectively, the error to be minimized is where ⊥ means orthogonal. Elsewhere,̂1( ) and̂2( ) are calculated according to (1). That implies that One possible solution is that each term of the equation is equal to zero. For the first term, we get The last equation implies that 1 ( ) ⊥ 2 ( ). For the other terms in (A.7), we obtain the same conclusion 1 ( ) ⊥ 2 ( ). Therefore,̂1( ) and̂2( ) are orthogonal if 1 ( ) and 2 ( ) are also orthogonal.
Therefore, a MISO Volterra model with two inputs could be treated as two SISO Volterra models if the two inputs are orthogonal. This result could be generalized for MISO Volterra model with inputs.