Technical Field

Embodiments of the present invention described herein relate to methods and systems relating to vibrational shear wave elastography, in particular, methods and systems for determining the direction of a travelling shear wave within tissue, analysing the quality of a detected shear wave and acquiring image data of a steady state shear wave field within tissue.

Background to the Invention and Prior Art

Imaging is an attractive tool for monitoring disease, e.g., tumour response to cancer therapy, as it can map physical, functional or molecular properties of the entire tumour, as opposed to biopsy which only samples small volumes of tissue. Clinically, it can allow non-invasive assessment of tumour response, at an early stage, improving patient management. Preclinically, xenograft tumour models are a vital tool allowing researchers to test novel cancer therapies and imaging can be used to non-invasively assess the treatment at multiple time points, providing greater insight into the progression (or regression) of disease, whilst improving animal welfare and reducing animal numbers. The identification of imaging biomarkers of tumour response that may be used both clinically or preclinically, requires in vivo validation.

It is now understood that elastic properties of the tumour play an important role in the development of cancer. For example, increased elastic modulus is characteristic of many solid tumours, and a driver of cancer progression and as such the elastic modulus of tissue is established as a diagnostic and prognostic biomarker of cancer, and there is increasing evidence to suggest these properties may be useful as biomarkers of response to cancer therapy.

Ultrasound shear wave elastography is currently used in the clinic and in the laboratory to image stiffness properties of tissue. In preclinical cancer research shear wave elastography can be used to image and monitor tissue changes in xenograft tumour models in response to treatment. The implementation of preclinical shear wave elastography techniques in xenograft models is challenging due to the small size (<lcm) of the tumour requiring imaging with a high degree of spatial resolution. Using an external source of vibration, the delivery of an appropriate mechanical stimulus is also challenging due to the limited access to the tumour. In recent years there has been substantial interest and development in shear wave elastography methods that can be used in the clinic for improved diagnosis and treatment monitoring.

Summary of the Disclosure

The present disclosure addresses the above, by providing methods and systems for determining the direction of a travelling shear wave within tissue suitable for high frequency travelling shear waves.

In view of the above, from a first aspect, the present disclosure relates to a method for determining the direction of a travelling shear wave within tissue, the method comprising: obtaining images of the tissue containing at least part of the travelling shear wave, the images obtained using a medical imaging modality; processing the images of the tissue to obtain a 2D or 3D amplitude and/or phase representation of the travelling shear wave; calculating a spatial shear wave autocorrelation function in dependence on the 2D or 3D amplitude and/or phase representation; and determining the direction of the travelling shear wave using the autocorrelation function.

Several advantages are obtained from embodiments according to the above-described aspect. For example, this method enables higher frequency (e.g., 500-1000 Hz) travelling shear waves to be used which achieve greater resolution than reverberant waves and do not necessarily require multiple vibration sources. This allows measurement of shear wave speed at higher frequencies with wavelengths similar to or below the dimension of the tumour. Determining the direction of the travelling shear wave is advantageous because this allows the shear wave speed to be accurately estimated.

In some embodiments, imaging the tissue comprises acquiring imaging data using a synchronised line-by-line scanning sequence.

In some embodiments, processing the images of the tissue to obtain a 2D or 3D amplitude and/or phase representation of the travelling shear wave comprises computing frame-to- frame localised displacement induced by a vibrational source used to generate the travelling shear wave.

In some embodiments, processing the images of the tissue to obtain a 2D or 3D amplitude and/or phase representation of the travelling shear wave further comprises computing temporal oscillation data for each pixel based on the frame-to-frame localised displacement. In some embodiments, processing the images of the tissue to obtain a 2D or 3D amplitude and/or phase representation of the travelling shear wave further comprises applying spectral filtering to the images by applying a fast Fourier transform in the temporal domain to the temporal oscillation data to obtain amplitude and/or phase data for each pixel.

In some embodiments, after processing the images of the tissue to obtain a 2D or 3D amplitude and/or phase representation of the travelling shear wave, the method further comprises applying spatial filtering to the images using a kernel median filter applied to the 2D or 3D amplitude and/or phase representation of the travelling shear wave.

In some embodiments, the kernel median filter has a size smaller than a wavelength of the travelling shear wave. This is advantageous as this reduces the risk of "blurring out" the shear wave phase data.

In some embodiments, the method further comprises: using the determined direction, fitting the autocorrelation function to a cosine function to determine the wavenumber; and determining the speed of the travelling shear wave using the determined wavenumber and the frequency of a vibrational source used to generate the travelling shear wave. This is advantageous as the wave speed through tissue correlates with the stiffness of the tissue because mechanical waves travel faster through stiffer tissue than through softer tissue. Therefore, by determining wave speed, the stiffness of the tissue can be inferred. The stiffness of the tissue gives diagnostic information about the status of disease, e.g., cancerous tumours will be harder than the surrounding tissue.

In some embodiments, the frequency of the travelling shear wave (and thus the vibrational source used to generate the travelling shear wave) is between 500 Hz and 1000 Hz. This is advantageous as high frequencies (e.g., 500 to 1000 Hz) achieve greater resolution than lower frequencies. The frequency range 500 to 1000 Hz is good for small preclinical tumours. Lower ranges may be useful for clinical tumours.

In some embodiments, the medical imaging modality comprises ultrasound or MRI.

In some embodiments, prior to the step of determining the direction of the travelling shear wave, the amplitude of the travelling shear wave is set to unity. This is advantageous because this renders the determination of the direction of the travelling shear wave dependent on phase, not the amplitude of the wave. In some embodiments, the method further comprises: using the determined direction, fitting the autocorrelation function to a cosine function; calculating a goodness of fit parameter which measures the fit of the autocorrelation function to the cosine function.

In some embodiments, the obtaining images step comprises generating the travelling shear wave within tissue using a vibrational source vibrating at a frequency and imaging the tissue using the medical imaging modality to obtain the images of the tissue.

From a second aspect, the present disclosure relates to a method for analysing the quality of a detected shear wave within tissue, the method comprising: generating a travelling shear wave within tissue using a vibrational source vibrating at a first frequency; detecting first vibrational energy measured at the first frequency; detecting total vibrational energy measured at all frequencies; and evaluating the quality of the travelling shear wave based on comparing the first vibrational energy to the total vibrational energy.

From a third aspect, the present disclosure relates to a method of acquiring image data of a steady state shear wave field within tissue, the method comprising: calculating a time interval based on at least a frequency of the steady state shear wave field and a desired number of samples; calculating a synchronisation delay required to maintain coherence between lines based on at least the desired number of samples, the time interval and the frequency; and acquiring a plurality of samples per image line in accordance with the calculated time interval and the calculated synchronisation delay.

This is known as a synchronised line-by-line scanning method. This method takes advantage of the continuous nature of the high frequency harmonic vibration which generates a steady state (non-transient) shear wave field. Unlike sequential A-line acquisition used in standard B-mode imaging, here we repeatedly scan individual A lines enabling the use of high pulse repetition frequency to avoid aliasing the high frequency shear waves. Subsequent A-lines transmissions are then synchronised in such a way that phase coherence is maintained at the shear wave frequency, eliminating the need for phase-unwrapping. The line-by-line approach works because the vibration source is continuously on, providing a "steady-state" shear wave field (as opposed to transient shear wave elastography) meaning we can reconstruct the shear wave field slowly one line at a time (relative to planar wave imaging) because the shear wave field should not be changing.

In some embodiments, the acquiring of the plurality of samples per image line comprises: performing n samples of a first image line with the time interval; and after a time period equal to the synchronisation delay, performing n samples of a second image line with the time interval.

From a fourth aspect, the present disclosure relates to a system for determining the direction of a travelling shear wave within tissue, the system comprising: a vibrational source configured to vibrate at a frequency to generate a travelling shear wave within tissue; a medical imaging modality configured to image the tissue to obtain images of the tissue; a processor; and a memory including computer program code; the memory and the computer program code configured to, with the processor, cause the processor to: process the images of the tissue to obtain a 2D amplitude and/or phase representation of the travelling shear wave; calculate a spatial shear wave autocorrelation function in dependence on the 2D amplitude and/or phase representation; and determine the direction of the travelling shear wave using the autocorrelation function.

The different embodiments described in relation to the first aspect apply equally to the fourth aspect.

From a fifth aspect, the present disclosure relates to a system for analysing the quality of a travelling shear wave within tissue, the system comprising: a vibrational source configured to vibrate at a first frequency to generate a travelling shear wave within tissue; a detector configured to detect first vibrational energy measured at the first frequency; a detector configured to detect total vibrational energy measured at all frequencies; a processor; and a memory including computer program code; the memory and the computer program code configured to, with the processor, cause the processor to: evaluate the quality of the travelling shear wave based on comparing the first vibrational energy to the total vibrational energy.

From a sixth aspect, the present disclosure relates to a system for acquiring image data of a steady state shear wave field within tissue, the system comprising: a processor; and a memory including computer program code; the memory and the computer program code configured to, with the processor, cause the processor to: calculate a time interval based on at least a frequency of the steady state shear wave field and a desired number of samples; calculate a synchronisation delay required to maintain coherence between lines based on at least the desired number of samples, the time interval and the frequency; and acquire a plurality of samples per image line in accordance with the calculated time interval and the calculated synchronisation delay. In some embodiments, the memory and the computer program code are further configured to, with the processor, cause the processor to automatically update the time interval and the synchronisation delay in response to changes in the frequency and/or the desired number of samples.

Brief Description of the

Embodiments of the invention will now be further described by way of example only and with reference to the accompanying drawings. Figures 2 to 30 are accompanied by the description in the appended papers. In particular, Figures 2 to 12 are explained further in the "The Use of Travelling Shear Waves in Preclinical Vibrational Shear Wave Elastography" paper. Figures 13 to 23 are explained further in the "High Frequency Ultrasound Vibrational Shear Wave Elastography for Preclinical Research" paper. Figures 24 to 30 are explained further in the "Preclinical three-dimensional vibrational shear wave elastography for mapping of tumour biomechanical properties in vivo" paper.

Figure 1 is a block diagram of a computer system according to an embodiment of the present invention.

Figure 2 is a schematic showing data processing steps involved in some direction detection embodiments. A series of 2D frames in the Z-X plane are used to compute frame-to-frame localised displacement (a) induced by a 1000 Hz vibration source. The temporal oscillation at a given pixel (b) is used to compute the spectrum (c), from which oscillation amplitude (d) and phase (e) can be extracted. Square, 3x3 mm kernels are used to compute the 2D autocorrelation function with spatial lags Az and Ax, shown in (e).

Figure 3 is an example of a 2D autocorrelation function in the Z-X plane with the polar coordinate distance from the centre position (Ax=Az=0), r, and angle 0 measured below the positive X axis (a). Example of interpolation to polar coordinates using 88 branches at uniform Ad and Ar intervals (black dots in (b)), isocontour level of 0.7 (white line shape marked A), with the minimum width 0 (solid lines marked B) and opposing angle (dashed lines marked C) (b-c). Interpolated AC values along 0 _min are shown in (d), with a cosine fit applied (dashed line).

Figure 4 is a B-mode image showing the stiff inclusion (arrow) in the breast phantom (a), phase map using a 700 Hz shear wave vibration (b), shear wave speed image computed using the directional detector with 5mm kernel (c), and dmin image showing the detected shear wave direction in degrees (d). The circle indicates the location of the stiff inclusion in (b), (c) and (d). Figure 5 shows the effect of varying AC interpolation angular interval Ad (90°, 45°, 20° and 10° from left to right respectively) on shear wave speed estimate (m/s, top row) and detected angle dmin (degrees, bottom row) using 5mm kernel and maximum spatial lag 1mm. The shear wave speed at three selected positions are shown: A in the stiff inclusion, B in the background, and C in a region with reverberation.

Figure 6 shows the effect of varying AC interpolation maximum spatial lag (0.5mm, 1.0mm, 1.5mm and 2.0mm respectively) on shear wave speed estimate (m/s, top row), and AC cosine fit mean normalised residual (middle row) using 36 angular steps (interval Ad = 10°) and 5mm kernel size. In the bottom row the linear AC profiles (solid lines) and their corresponding fits (dashed lines) used to determine the shear wave speed from three positions (A, B, C) in Figure 4 are shown: stiff inclusion (A), background (B) and the reverberant region (C).

Figure 7 shows the effect of varying AC kernel size (2mm, 3mm, 4mm, 5mm, 6mm, 8mm and 10mm respectively from left to right, and from top to bottom) on shear wave speed (m/s) images, AC interpolation using 36 angles (interval Ad = 10°) and 1mm maximum spatial lag, A circular ROI centred within the stiff inclusion (inner circle A), and an annulus in the background region (B) are shown. The mean (bold line) and standard deviation (dashed line) of shear wave speeds (centre bottom plot) from within the stiff inclusion (square), and background annulus region (circle) are used to estimate contrast to noise ratio (bottom right plot).

Figure 8 shows the effect of median filtering kernel size on shear wave phase (row a), shear wave speed (row b), cosine fit normalised residual (row c) for no median filtering, 0.7mm, 1.3mm and 2.5 mm kernel size respectively from left to right. Example autocorrelation profiles for three positions in the phantom are shown (row d) : stiff inclusion (B), background (C) and noisy region (A). The shear wave speed and normalised residuals for all three positions are plotted in figures e and f respectively.

Figure 9 shows AC profile (solid line) and cosine fit (crosses) applied to the data in the stiff inclusion (left column), and two further locations in the background (middle and right hand side columns from) in the CIRS phantom for kernel size of 3, 5, 8, 11 and 15 mm respectively (rows, from top to bottom). The effect on the best fit of varying r _max is indicated by the marker shape (o - 1mm; A - 2mm; v - 3mm; * - 4mm; x - 5mm). Figure 9 demonstrates that good shear wave speed estimates are obtained with small kernel sizes when using travelling shear waves. Figure 10 shows phase maps in a 500 Hz reverberant shear wave field in a homogeneous region of the CIRS phantom (a-e) showing kernel sizes of 3, 5, 8, 11 and 15 mm respectively. The real part of the autocorrelation maps are shown (f-j) for the corresponding kernel sizes, together with a dashed line indicating the result of the direction detector when applied to the reverberant field data.

Figure 11 shows a cosine fit (left column) applied to the data shown in Figure 10 for kernel size of 3, 5, 8, 11 and 15 mm respectively, from top to bottom. Plots in the middle and right hand columns show the reverberant field fits [Zvietcovich et al] applied to data in the axial (Z) and lateral (X) data from the autocorrelation data shown in Figure 9, for increasing kernel size. The effect of varying search length is indicated by the marker shape (o - 1mm; A - 2mm; v - 3mm; * - 4mm; x - 5mm). Figure 11 demonstrates larger kernel sizes are required to obtain accurate shear wave speed estimates under reverberant shear wave conditions.

Figure 12 shows a directional detector applied to vibrational shear wave elastography data acquired from a xenograft mouse tumour ex-vivo. The B-mode image (a), phase map (b), shear wave speed image (c), and 6 _m!n image (d) are shown. Vibration frequency was 1000 Hz, kernel size and r _max were set to 2 and 0.75 mm respectively.

Figure 13 is a diagram illustrating the experimental setup for vibrational shear wave elastography in the xenograft tumours.

Figure 14 is a schematic indicating the delays between transmissions events for the line- by-line imaging sequence for a 1000 Hz shear wave signal. The bold line (reference numeral 1402) represents temporal segments of the harmonic shear signal that is measured by the system. The thinner line (reference numeral 1404) which is beneath the bold line represents the underlying harmonic shear signal we wish to detect but is only visible in the figure between 8.5 ms and 9 ms. The gap, where the thinner line is visible, is an interval during which no sampling takes place where we are switching acquisition to the next image A line. The thinner line indicates the harmonic shear wave signal at fixed location. A regular short delay, 6t, is used to sample a section (bold line) of the shear wave oscillations at regular intervals (6x shear wave frequency). A longer time delay, At, is used between the last transmission of the current A-line, and the first transmission event of the next A-line, so phase coherence is maintained between A-lines. The value of At is automatically calculated based on the number of samples to be acquired, the value of 6t, and the shear wave frequency. Figure 15 shows B-mode images (a-c) and signal to noise ratio images [dB] (d-f) of a uniform section of the CIRS breast elastography phantom imaged with the L22-14vX corresponding to planar wave imaging with no transmit apodisation (a, d), planar wave imaging with Hanning transmit apodisation (b, e), and focused beam line-by-line imaging (c, f).

Figure 16 shows a detected 1000 Hz shear wave oscillation phase [rad] (a, b), Conformance [%] (c, d), and estimated shear wave speed [m/s] (e, f) for planar wave imaging with Hanning transmit apodisation (a, c, e), and focused beam imaging (b, d, f) in the CIRS phantom.

Figure 17 shows detected 600 Hz shear wave oscillations: time domain (a) and spectrum (b) for location with signal Conformance greater than 99%; time domain (c) and spectrum (d) for location with signal Conformance of 60%.

Figure 18 shows shear wave detection in an ex-vivo tumour with focused beam imaging. Example B-mode image with the contactor indicated by the single-headed arrow (a) and the direction of motion of the contactor (double-headed arrow), phase [rad] (b), and Conformance [%] (c) with a shear wave frequency of 500 Hz. Phase images (d-f), and Conformance images (g-i) are shown for 600, 800 and 1000 Hz shear wave frequencies respectively.

Figure 19 shows estimated shear wave speed [m/s] in ex-vivo tumour for frequencies of 500, 600, 800 and 1000 Hz (columns, left to right), using a Ixlmm, 2x2mm and 3x3mm kernel (rows, from top to bottom).

Figure 20 shows mean shear wave speed in ex-vivo tumour vs vibration frequency estimated with Ixlmm, 2x2mm and 3x3mm kernel, error bars indicate standard deviation (a). Shear wave speed histograms using 1mm (b), 2mm (c) and 3mm (d) square kernels for shear wave speed estimation, horizontal bars represent the mean ± one standard deviation for each vibration frequency.

Figure 21 shows B-mode outline of three tumours (rows) imaged in vivo using VSWE. Outer outline includes the whole tumour including hyper-echoic rim, inner outline indicates body of the tumour. Columns represent images obtained at 500, 700 and 1000 Hz vibration frequency respectively. The contactor was positioned towards the top right of the tumour in all images. The effect of tissue motion can be observed by the apparent undulation of the skin surface, examples of which are indicated by the arrows. Figure 22 shows VSWE Conformance images for in-vivo tumours corresponding B-mode images shown in Figure 21. Conformance images represent different tumours (rows) and vibration frequencies (columns) of 500, 700 and 1000 Hz respectively. Outlines are those drawn using the B-mode images, as indicated in Figure 21.

Figure 23 shows VSWE shear wave speed images for in-vivo tumours corresponding to the B-mode and Conformance images shown in Figures 21 and 22 respectively. Shear wave speed imaged represent different tumours (rows) and vibration frequencies (columns) of 500, 700 and 1000 Hz respectively. Outlines are those drawn using the B-mode images, as indicated in Figure 21.

Figure 24 shows the design of our 3D-vibrational shear wave elastography (3D-VSWE) platform, (a) schematic showing the apparatus used to acquire vibrational shear wave elastography data. Two shakers were used as vibrational sources and were coupled to the tumours using carbon fibre rods and contactors. Shear wave fields were imaged using line-by-line focused imaging acquired with a L22-14vX linear array ultrasound probe (Philips) and a Vantage system (Verasonics). The probe was scanned superiorly in a step and shoot motion to collect 3D data, (b) Representative B-mode image data and phase maps acquired using vibrational frequencies of 700, 1000 and 1200 Hz. Images are the central slices of 3D datasets in the transverse, sagittal and coronal planes.

Figure 25 shows the performance of our 3D-vibrational shear wave elastography (3D- VSWE) platform, (a) shows orthogonal slices through 3D maps of c _s , C and Goodness of Fit (GoF) across the tumour. The Conformance is a measure of the percentage of the vibrational frequency detected at the vibrational source frequency and can be mapped to allow the user to assess the penetration of the shear waves during positioning of the contactors. Similarly, the average GoF across the tumour allows real-time assessment of the quality of the data that can be obtained, (b) and (c) show the percentage change in spatial-median c _s and slope, respectively, between repeat measurements acquired 24 hours apart. Data shown were generated using 2 mm kernels with 0.5 mm spacing and one vibrational source.

Figure 26 illustrates how 3D-vibrational shear wave elastography is sensitive to ZD6126- induced tissue damage, (a and b) show MDA-MB-231 human triple negative breast cancer, (d and e) show U-87 MG human glioblastoma tumours. Representative 3D-VSWE derived parametric maps of c _s shows a marked and global decrease in c _s , 24h following treatment with ZD6126 in both models, corroborated by the presence of extensive regions of haemorrhagic necrosis (outlined in black) visible on H&E-stained slide in stark contrast to the presence of large area of viable tumour with varied number of necrotic foci in the vehicle-control cohort, (b and e) show the quantitative analysis of the relative changes in spatial-median c _s values acquired at different frequencies, with either 200mg/kg ZD6126 or vehicle, in both models (box plot show mean ± SD; Mann-Witney test with 5% significant level), c and f show the quantitative analysis of the extent of necrosis quantified on H&E-stained section.

Figure 27 shows comparison of 3D-VSWE derived parametric maps of shear wave speed c _s with aligned haematoxylin and eosin stained histological sections, (a) shows MDA-MB- 231 and (b) shows U-87 MG vehicle tumours. Necrotic regions in the tissue micrographs have been coloured black; in the c _s maps the corresponding necrotic regions tend to have low c _s represented by dark colours and black. The spatial distribution of necrotic and viable tissue in vehicle tumours reveals that VSWE provides a good representation of the heterogeneous viable versus necrotic tissue within the tumour. Note that c _s maps are the central plane of the 3D data which represents a 2 mm thick slice of tissue whereas the histological sections were about 8 pm thick.

Figure 28 illustrates how choice of a kernel dimension of 2 mm was based on the analysis of the shear wave speed (c _s ), Goodness of Fit and repeatability, for kernel sizes of 0.5, 1, 2 and 3 mm. (a) shows shear wave speed was similar for all kernel sizes, (b) On average GoF increased vibrational frequency and decreased with kernel size, (c) Here the repeatability was assessed using the percentage change in c _s between repeat measurements acquired 24 hours apart, (d) The conformance, which was the percentage of vibrational energy detected at the shear wave drive frequency given varied little with kernel size as expected because it is assessed at the centre of the kernel. A kernel dimension of 2 mm offered a good repeatability and small GoF error. Data shown are for one shaker. Trends were similar for data acquired with two shakers.

Figure 29 shows box plots and whisker plots of c _s , C and GoF in MDA-MB-231 tumours as measured on three consecutive days. We compared median c _s , mean C and mean GoF for one and two vibrational sources (shakers) from all tumours (n = 5). There was no significant difference in C and GoF for two vibrational sources compared to one. Median c _s was greater for two shakers but not significantly so.

Figure 30 shows representative histology sections stained with H8<.E and CD31 antibody. MDA-MB-231 and U-87 MG tumours treated with ZD6126 and vehicle are shown. Dashed black line shows necrotic region. Inset has been magnified to show necrotic and viable cell appearances. CD31 staining has been shown to demonstrate the functionality of ZD6126 on tumour.

Figure 30 is a flow diagram in accordance with embodiments of the present invention.

Figure 31 is a flow diagram in accordance with embodiments of the present invention.

Figure 32 is a flow diagram in accordance with embodiments of the present invention.

Figure 33 is a flow diagram in accordance with embodiments of the present invention.

Figure 34 is a flow diagram in accordance with embodiments of the present invention.

Figure 35 is a flow diagram in accordance with embodiments of the present invention.

Figure 36 is a flow diagram in accordance with embodiments of the present invention.

Figure 37 is a flow diagram in accordance with embodiments of the present invention.

Figure 38 is a flow diagram in accordance with embodiments of the present invention.

Figure 39 is a flow diagram in accordance with embodiments of the present invention. ion of the Embodiments

Overview

Embodiments of the present invention relate to methods of vibrational shear wave elastography. In particular, some embodiments relate to determining the direction of a travelling shear wave within tissue, e.g., tumours. Some embodiments relate to determining the direction of a travelling shear wave within tissue using only phase information (not the amplitude of the wave) such that the determination of the direction of the travelling wave is dependent on phase only. Some embodiments relate to analysing the quality of the shear wave within the tissue.

Some embodiments of the present invention are suitable for use with monochromatic shear waves at frequencies higher than those conventionally used in the clinic to generate travelling waves in tumours. Monochromatic shear wave data may include monochromatic shear wave data derived via Fourier transform filtering of polychromatic shear wave datasets. Elastography maps elastic properties and stiffness of soft tissue. The hardness of the tissue gives diagnostic information about the status of disease, e.g., cancerous tumours will be harder than the surrounding tissue. Elastography can be used to guide or even replace biopsies. Elastography is non-invasive. Elastography works by inducing a distortion in the tissue, observing and processing the tissue response to the distortion to infer mechanical properties of the tissue, and then displaying results to the operator, usually in the form of an image. Usually, the observation of the tissue response is displayed to the operator along side a conventional image of the tissue, so that measurements of the tissue response can be mapped to a location within the tissue. Stiffness of the tissue can be calculated based on the principle that mechanical waves travel faster through stiffer tissue than through softer tissue. Being able to detect the wave speed through tissue is therefore important.

Prior art techniques use lower frequency (e.g., <300 Hz) reverberant waves to ensure waves can effectively reverberate around the tumour. However, these do not provide high enough resolution to be used with small (e.g., <10 mm) tumours as the wavelength (>10 mm) associated with these lower frequency waves is typically larger than the size of the tumour, and hence it is not possible to measure shear wave speed on the scale of a few millimetres. To overcome this problem higher vibration frequencies may be used, however higher shear wave attenuation impedes effective reverberant conditions, resulting in attenuated travelling shear waves. The use of multiple vibration sources to overcome the effects of attenuation is both costly and practically challenging due to the limited access to the tumour.

To solve these problems, some embodiments of the present invention use high frequency (e.g., 500-1000 Hz) travelling shear waves which achieve greater resolution than reverberant waves and do not necessarily require multiple vibration sources. This allows measurement of shear wave speed at higher frequencies with wavelengths similar to or below the dimension of the tumour.

The response of the shear wave speed estimator in a phantom may be characterised in terms of kernel size, maximum spatial lag, angular interval and/or spatial median filtering kernel size. Reliable shear wave speed estimates can be obtained using travelling waves with small kernel sizes, as opposed to a prior art reverberant shear wave method which is more readily implemented using lower vibration frequencies and benefits from kernel sizes larger than the shear wavelength. The angular interval is preferably less than 20° in order to get better precision in the directional detector, thus reducing error in shear wave speed estimates. Reducing the interval further provides diminishing improvements in the precision of shear wave speed estimates.

The maximum spatial lag preferably should be chosen such that it covers the central peak of the autocorrelation function down to 0. The exact value that achieves this will depend on the wavelength we expect to measure. Preferably the maximum lag to be used is set adaptively during processing.

In one example, vibrational shear wave elastography in accordance with some embodiments of the present invention may be used in preclinical tumours by using travelling shear waves at frequencies up to and potentially above 1000 Hz, for example, frequencies between 500 and 1000 Hz. The required shear waves can be generated using a single vibration source. A shear wave speed estimator is based on experimental evaluation of the localised shear wave field spatial autocorrelation function. The direction of the travelling shear wave, a property that cannot be assumed a-priori with arbitrary positioning of the vibration source, can be determined from the autocorrelation function. By using the spatial autocorrelation function as the basis of the measurement, alternative estimates of shear wave speed may be possible if simple travelling wave conditions are not met. An example of an alternative estimate of shear wave speed is the reverberant shear wave method, where under truly reverberant conditions the autocorrelation function would be expected to follow a theoretically derived profile for a shear wave imaging system sensitive to displacement in the axial direction only. In this approach the theoretical profile can be fitted to the experimental autocorrelation function in 2D, or perhaps in ID (axial, lateral, or both) to determine an estimate of wavelength (see Figures 10 and 11). Incidentally, there is no need for a direction detector in this reverberant approach because by definition a reverberant field comprises multiple waves randomly and evenly distributed in all directions.

The direction detection approach in accordance with some embodiments assumes non- reverberant, locally-planar shear wave propagation, and the kernel size can be smaller than the shear wavelength, allowing for localised estimation of shear wave speed, i.e., better resolution. 2D autocorrelation data may be used to determine the local direction of the shear wave in the region defined by the kernel. Subsequently, a ID cosine curve may be fitted only to the real part of the autocorrelation data aligned with the direction of the shear wave. The relatively small size of the kernel therefore inherently allows for any variation in shear wave propagation through the imaging plane due to diffraction. The cosine fitting process yields a shear wave wavelength estimate, A The wavelength, and vibration frequency f are then used to calculate the shear wave speed, c: c = A f .

Autocorrelation is a known mathematical tool for finding repeating patterns. The autocorrelation function is used to determine the direction of the travelling shear wave before estimating shear wave speed. This is necessary because the source of the vibration and the interaction of the shear wave with the tumour tissue means no prior assumption can be made regarding the direction of the shear wave, unlike for example in Acoustic Radiation Force Impulse (ARFI) based transient shear wave elastography methods where the axial 'push' beam is expected to generate shear waves propagating laterally. The use of the autocorrelation function is explained in detail in the subsection titled "Shear Wave Speed Estimation" of the appended paper "The Use of Travelling Shear Waves in Preclinical Vibrational Shear Wave Elastography". The autocorrelation function allows the coordinate of 0 which minimised the corresponding r coordinate (referred to as 0 _mzn ) to be determined. dmin is the direction of the travelling shear wave.

Once the direction of the travelling shear wave has been determined, the wavenumber k _] can be determined. Frequency f is known from the vibration frequency. This allows the shear wave speed c to be obtained using

Ci.i = Zitf/kt

Various aspects and details of principal components of embodiments of the present invention will be described below with reference to Figures 1 to 37.

An example of a computer system used to perform embodiments of the present invention is shown in Figure 1.

Figure 1 is a block diagram illustrating an arrangement of a system according to an embodiment of the present invention. Some embodiments of the present invention are designed to run on general purpose desktop or laptop computers. Therefore, according to an embodiment, a computing apparatus 100 is provided having a central processing unit (CPU) 106, and random access memory (RAM) 104 into which data, program instructions, and the like can be stored and accessed by the CPU. The apparatus 100 is provided with a display screen 120, and input peripherals in the form of a keyboard 122, and mouse 124. Keyboard 122, and mouse 124 communicate with the apparatus 100 via a peripheral input interface 108. Similarly, a display controller 105 is provided to control display 120, so as to cause it to display images under the control of CPU 106. Image data 102 from a medical imaging modality 150 and/or frequency data 103 from the vibrational source 152, can be input into the apparatus and stored via data input 110. In this respect, apparatus 100 comprises a computer readable storage medium 112, such as a hard disk drive, writable CD or DVD drive, zip drive, solid state drive, USB drive or the like, upon which image data 126, frequency data 128, images of tissue 142, shear wave data 130, which may comprise displacement data and/or phase data, and wave direction and/or wave speed data 132 can be stored. Alternatively, data 126, 128, 130, 132, 142 could be stored on a web-based platform, e.g. a database, and accessed via an appropriate network. Computer readable storage medium 112 also stores various programs, which when executed by the CPU 106 cause the apparatus 100 to operate in accordance with some embodiments of the present invention.

In particular, a control interface program 116 is provided, which when executed by the CPU 106 provides overall control of the computing apparatus, and in particular provides a graphical interface on the display 120 and accepts user inputs using the keyboard 122 and mouse 124 by the peripheral interface 108. The control interface program 116 also calls, when necessary, other programs to perform specific processing actions when required. For example, an image acquisition program, which may be a line-by-line scanning program 136 or a plane wave imaging program, may be provided which is able to operate on image data 126 (which may be ultrasound echo data, referred to here as "image data" because it is data which is used to create the images) indicated by the control interface program 116, so as to output images of tissue 142. An image processing program 140 may be provided which is able to operate on the images of tissue 142 indicated by the control interface program 116, so as to output shear wave data 130. A shear wave data processing program 134 may be provided which is able to operate on the shear wave data 130 and frequency data 128 indicated by the control interface program 116, so as to output wave direction and/or wave speed data 131. A quality assessment program 138 may be provided which is able to operate on the shear wave data 130 and/or wave direction and/or wave speed data 132 and/or frequency data 128 to output quality metric data. The operations of the line-by-line scanning program 136, shear wave data processing program 134, image processing program 140 and quality assessment program 138 are described in more detail below.

The detailed operation of the computing apparatus 100 will now be described. Firstly, the user launches the control interface program 116. The control interface program 116 is loaded into RAM 104 and is executed by the CPU 106. The user then launches a program 114, which may be comprised of the line-by-line scanning program 136, shear wave data processing program 134, image processing program 140 and quality assessment program 138. The programs act (directly or indirectly) on the input data (image data 102, 126, and frequency data 103, 128) as described below.

The Direction Detector

Some embodiments of the present invention relate to a direction detector which detects the direction of travel of the travelling shear wave which enables the speed of the travelling shear wave within tissue to be determined. The direction detector may comprise one or more of the image acquisition program, e.g., line-by-line scanning program 136, the shear wave data processing program 134, the image processing program 140 and the quality assessment program 138.

With reference to the Figure 31, embodiments of the present invention (referred to as the direction detector) may begin with obtaining images of the tissue containing at least part of the travelling shear wave, the images obtained using a medical imaging modality. This obtaining step may comprise generating a travelling shear wave within tissue using a vibrational source 152 vibrating at a frequency (step 3102). One example of an experimental set up 1300 is shown in Figure 13 which shows a mini-shaker 1302 as the vibrational source 152 in contact with tissue 1306 (in this example, mouse tissue) via a contactor 1304. Frequency data 103 (e.g., the frequency of vibration) from the vibration source 152 (e.g., 1302) is input into the computer apparatus 100 and stored as described above.

The tissue is then imaged using a medical imaging modality 150 to obtain images of the tissue (step 3104). The medical imaging modality 150 may comprise ultrasound. In the example shown in Figure 13, the L22-14vX imaging probe 1308 is an ultrasound probe. Examples of ultrasound images 202 are shown in Figure 2a. Image data 102 from the medical imaging modality 150 is input into the computer apparatus 100 and stored as image data 126 as described above. The image data 126 may be used to obtain images of tissue 142 via a planar wave imaging sequence with or without angular compounding, a synchronised line-by-line scanning sequence, or any other suitable image acquisition sequence. Embodiments of the invention are applicable to any ultrasound beamforming/image acquisition methods. The planar wave imaging sequence or synchronised line-by-line scanning sequence may be performed by a program, for example, the line-by-line scanning program 136 described above. This program (e.g., line- by-line scanning program 136) takes image data 126 as inputs and outputs images of tissue 142 which can be stored as described above.

The images of the tissue 142 are then processed by an image processing program 140 to obtain a 2D or 3D amplitude and/or phase representation of the travelling shear wave (step 3106), also referred to as "shear wave data 130". An example of a 2D amplitude representation of the travelling shear wave is shown in Figure 2d and an example of a 2D phase representation of the travelling shear wave is shown in Figure 2e.

In some embodiments, the 2D or 3D amplitude and/or phase representation of the travelling shear wave (shear wave data 130) is obtained by the image processing program 140 as follows (see Figure 32). The frame-to-frame localised displacement induced by the vibrational source 152 is computed (step 3202), see Figure 2a which shows different frames and how the shear wave moves across the frames (images 202). The temporal oscillation data (shown in Figure 2b) for each pixel is computed based on the frame-to- frame localised displacement (step 3204). Each pixel (e.g., the circled pixel in Figure 2a) has its own temporal oscillation data which can be displayed as in Figure 2b, effectively illustrating how the travelling shear wave moves through each pixel in time. Spectral filtering (step 3206) is then applied to the images by applying a fast Fourier transform in the temporal domain to the temporal oscillation data (e.g., the data of Figure 2b). This results in the spectrum (containing amplitude and/or phase data) shown in Figure 2c for each pixel. In Figure 2c it can be seen there is a peak around the vibration frequency (1000 Hz). This peak has an amplitude (approximately 3.8 mm/s) for this example pixel. This amplitude value may then be plotted at the location of that pixel. This may be done for every pixel to obtain a 2D or 3D amplitude representation of the travelling shear wave (an example of a 2D amplitude representation is shown in Figure 2d). A 2D or 3D phase representation of the travelling shear wave (an example of a 2D phase representation is shown in Figure 2e) may be extracted from the spectrum shown in Figure 2c. In this way, the amplitude and/or phase data for each pixel is processed to obtain a 2D or 3D amplitude and/or phase representation of the travelling shear wave (step 3208).

A spatial shear wave autocorrelation function is then calculated by a shear wave data processing program 134 in dependence on the 2D or 3D amplitude and/or phase representation (step 3108). The spatial autocorrelation provides a measure of similarity between a localised region of the shear wave field, defined by a kernel with known size, and an equally sized region displaced in 2D or 3D space by a known amount (lag). An example of a 2D autocorrelation function is shown in Figure 3a. In the case of a travelling wave, the result of the spatial autocorrelation is such that when the lag follows a direction orthogonal/normal to the wave direction, a high degree of correlation is maintained because there is little to no variation in the underlying shear wave phase. When the lag is in the direction parallel to wave propagation, the autocorrelation varies rapidly as a result of the gradient in shear wave phase (lines B and C in Figures 2b and 2c). The kernel defines the dimensions of the localised region of the shear wave field to be used for computing the autocorrelation. A small kernel size thus provides better spatial resolution, and a larger kernel size provides more spatial averaging which may be beneficial if the detected shear wave field is noisy. The direction of the travelling shear wave is then determined by the shear wave data processing program 134 by determining the direction which minimises the width of the autocorrelation function (step 3110) starting from the origin (zero lag). To facilitate direction detection the autocorrelation data which is initially sampled in a 2D or 3D cartesian coordinate system Figure 2b), is converted to a polar coordinate system (Figure 2c), an example of determining the direction is shown in Figure 3. The mathematics of determining the direction are described in the subsection titled "Shear Wave Speed Estimation" of the appended paper "The Use of Travelling Shear Waves in Preclinical Vibrational Shear Wave Elastography". The direction data may be stored as wave direction data 132.

In 3D embodiments, the autocorrelation data in 3D coordinates is transformed to spherical polar coordinates (radius r, and two angles: 0 and < >). 3D direction detection is performed in the same way as 2D by determining the 0 and < > angles giving the least "width" in the autocorrelation function, and then fitting the cosine function in this direction only (identical to 2D).

The method of Figure 33 is the same as the method described in relation to Figure 31, with the addition of the step of spatial filtering. After processing the images of the tissue to obtain a 2D or 3D amplitude and/or phase representation of the travelling shear wave, the method may further comprise applying spatial filtering to the images using a kernel median filter applied to the 2D or 3D amplitude and/or phase representation of the travelling shear wave (see step 3308). The 2D or 3D median filtering is applied separately to the real and imaginary components of the detected shear wave image data, which define the amplitude and phase of the shear wave field. At a given pixel position the median filtered (post-processed) real and imaginary components are determined by computing the median values from the detected (pre-processed) values at the pixel location and its neighbours, as defined by the median filtering kernel size and shape. Applying median filtering to the detected shear wave is advantageous as it removes high frequency noise prior to computation of the autocorrelation function, thus increasing the accuracy of the direction detection. Without median filtering the autocorrelation function that is computed can be considered as resulting from the combination of the shear wave field data and any noise source. If the latter dominates, for example if no shear wave field is present, or the ultrasound echo signal is very noisy, the autocorrelation function will become highly dependent on the form of the noise signal, manifesting as a sharp spike centred around zero spatial lag if the nature of the noise signal is uncorrelated between adjacent image pixels. Median filtering is a simple and effective way to reduce the effects of noisy shear wave data, aiding recovery of autocorrelation function data that provide an improved fit to the cosine function, as shown in Figure 8. The effect of varying the median filtering kernel size is shown in Figure 8 for a number of examples in terms of the phase, shear wave speed, and normalised residual images. Alternative method to reduce noise in shear wave data may include spatial finite impulse response filters, k-space or spatial frequency domain filtering, or other data acquisition strategies aimed at improving shear wave signal to noise ratio by increasing the duration of the acquired shear wave signal, or for example by employing suitably angular compounding methods.

The determined direction of the wave may be used to determine the speed of the wave. The method of Figure 34 is the same as the method described in relation to Figure 31, with the further steps which allow the speed of the wave to be determined. In step 3412, the determined direction (wave direction data 132) is used to extract a ID profile from the autocorrelation function and fit to a cosine function to determine the wavenumber. In step 3414, the speed of the travelling shear wave is determined using the wavenumber and the frequency of the vibrational source 152 (frequency data 128). The speed of the travelling shear wave may be stored as wave speed data 132. The mathematics of these steps are described in the subsection titled "Shear Wave Speed Estimation" of the appended paper "The Use of Travelling Shear Waves in Preclinical Vibrational Shear Wave Elastography".

The method of Figure 35 is the same as the method described in relation to Figure 31, except for the additional step of setting the amplitude of the travelling shear wave to unity (step 3510). In some embodiments, prior to evaluation of the autocorrelation function, the amplitude of the wave is set homogeneously to unity, rendering shear wave speed estimation dependent on the phase data only. This is advantageous over using phase and amplitude of the shear wave because amplitude is more susceptible to noise than phase is.

Quality measures When using shear wave elastography to measure elastic properties of tissue, it is important for the user to know in real-time whether they are applying the correct frequency for the tissue and whether the vibrational source 152 is contacting the tissue in the correct place. The quality measures described below provide real-time feedback on whether the positioning and frequency of the vibrational source 152 is correct. Without these quality measures, the user would not know whether the elastography measurements had been successful until analysis was complete. The quality measures allow the user to modify the frequency and/or position of the vibrational source 152 during the elastography measurement to ensure a high-quality measurement is obtained without having to repeat the process later on after analysis has been performed.

The calculation of one or more of the quality measures described below may be performed by a quality assessment program 138 which operates as described in relation to computer apparatus 100 above.

Some embodiments of the present invention relate to real-time metrics of the shear wave field quality within a tumour to allow users to rapidly optimise the arrangement of vibrational sources and their frequencies and amplitude.

Two quality metrics (Figure 25a) have been devised which can be mapped and displayed to the user in real-time to aid with placement of the contactors and the choice of vibrational frequency:

1. The Goodness of Fit

The Goodness of Fit (GoF), a measure of the mean residual error of the cosinusoidal fit to the autocorrelation function. GoF is calculates as the root mean square of the difference between the measured autocorrelation profile and the best fit:

Where y _m and yt are the measured and fitted autocorrelation values, and N is the number of samples used to determine the fit.

GoF error indicates how well the shear wave phase and magnitude (over a cuboid kernel) conform to a travelling shear wave and provided a measure of the quality of the estimation of shear wavelength. Figures 6, 8c and 8f illustrate uses of the GoF quality measure. The GoF measure can be used to assess the effect of varying AC interpolation maximum spatial lag on sheer wave estimate (see Figure 6).

The GoF measure can also be used to assess the effect of median filtering kernel size on the fit of AC profiles to cosine fits (see Figure 8). In Figure 8, autocorrelation profiles and corresponding cosine fits are included for locations of interest. These profiles indicate the poor fit obtained when no median filtering is applied, with relatively large normalised residuals. An improved fit was obtained with a median filtering kernel size >0.7mm for the stiff inclusion and background regions. For location C (right-most location in Figure 8) corresponding to a position more distant from the contactor, and showing evidence of shear wave reflection, the choice of kernel size had a large influence on the resulting autocorrelation profile, for example at the largest kernel size tested of 2.5mm the shear wave speed was almost equal to the estimates in the stiff inclusion region. These results therefore indicate the choice of spatial median filtering kernel size may impose a degree of bias in shear wave speed estimates depending on the signal to noise ratio in shear wave detection, or relative conformance to the assumption of plane progressive shear waves.

The method for evaluating the GoF is described in relation to Figure 36. The method of Figure 36 is the same as the method described in relation to Figure 31, except for the additional steps of fitting the autocorrelation function to a cosine function using the determined direction (step 3612) and calculating a goodness of fit parameter which measures the fit of the autocorrelation function to the cosine function (step 3614). One or more of these steps may be performed by the quality assessment program 138.

2. The Conformance

The Conformance (C), enables evaluation of the quality of the detected shear wave field relative to vibrations at other frequencies and noise. C is calculated as the percentage of the total detected vibrational energy that is measured at the shear wave drive frequency: where |S/| represents spectral amplitude components, and |Sf| is the amplitude component at the shear wave frequency. The detection of shear waves depends on the amplitude of the oscillation being large enough to be detected, but also on the signal to noise ratio of the ultrasound echo data. The Conformance measure therefore represents a simple and effective method to visualise detection of shear waves that combines the factors above.

The method for evaluating the Conformance is described in relation to Figure 37. At step 3702, a travelling shear wave is generated within tissue using a vibrational source 152 vibrating at a first frequency. At step 3704, first vibrational energy measured at the first frequency (i.e., energy at the frequency of vibration of the vibrational source, said frequency referred to as "the shear wave drive frequency") is detected. At step 3705, the total vibrational energy measured at all frequencies is detected (e.g., total vibrational energy). At step 3708, the quality of the travelling shear wave is evaluated based on comparing the first vibrational energy to the total vibrational energy. One or more of these steps may be performed by the quality assessment program 138.

An example of detected shear wave oscillations is given in Figure 17, comparing the time domain signal and spectra from regions with high (>99%) and low (60%) Conformance.

Examples of phase and Conformance images obtained in an ex-vivo tumour are shown in Figure 18, for vibration frequencies from 500 Hz to 1000 Hz. With the contactor placed on the tumour skin, as indicated by the single headed arrow in Figure 18a, it is possible to visualise the penetration of the shear wave into the tumour. Higher amplifier gain settings were required at the higher frequencies. Amplifier gain was increased within a safe operating limit to improve detection of the shear waves, as determined by inspection of the real-time Conformance image. At a frequency of 500 Hz shear wave oscillations are detected throughout the tumour. As the frequency is increased the Conformance images show a progressive reduction in the ability to detect shear wave oscillations due to increased shear wave attenuation. At 600 Hz Conformance is reduced towards the distal edge of the tumour with respect to the position of the contactor. At 800 Hz the shear waves do not propagate further than 8-10 mm from the contactor, and at 1000 Hz it is difficult to detect any shear waves beyond 10 mm from the contactor. The phase images in the tumour show a progressive shear wave similar to the ones observed in the CIRS phantom at the higher frequencies of 800 and 1000 Hz. At the lower frequencies there is more evidence of possible interactions between the shear wave propagating into the tumour reflections from the tumour surface.

The B-mode images of three tumours at three vibration frequencies are shown (Figure 21), with their respective Conformance (Figure 22). Outlines of the tumour including a hyper-echoic rim around the edge (darker line) and the body of the tumour (white line) are included. The in vivo Conformance images replicate the ex-vivo findings, namely that at 500 Hz it is possible to generate and detect shear vibrations throughout the tumour, and increasing the frequency reduces the effective penetration of the shear waves within the body of the tumour.

Synchronised line-by-line processing

Vibrational shear wave elastography may be implemented with a high frequency (e.g., MHz frequencies) probe running a focused line-by-line ultrasound imaging sequence. This is advantageous as it has been found to offer improved detection of 1000 Hz shear waves over an ultrafast planar wave imaging sequence in a homogenous tissue-mimicking phantom.

To reconstruct an image and detect shear waves, embodiments of the invention use a synchronised line-by-line scanning method which can be performed by the line-by-line scanning program 136. This method differs from the prior art as it takes advantage of the continuous nature of the high frequency harmonic vibration which generates a steady state (non-transient) shear wave field. Unlike sequential A-line acquisition used in standard B- mode imaging, here we repeatedly scan individual A lines enabling the use of high pulse repetition frequency to avoid aliasing the high frequency shear waves. Subsequent A-lines transmissions are then synchronised in such a way that phase coherence is maintained at the shear wave frequency, eliminating the need for phase-unwrapping.

With reference to Figure 39, a method of acquiring image data of a steady state shear wave field within tissue is described. At step 3902, a time interval is calculated based on at least a frequency of the steady state shear wave field and a desired number of samples. At step 3904, a synchronisation delay required to maintain coherence between lines is calculated based on at least the desired number of samples, the time interval and the frequency. At step 3906, a plurality of samples is acquired per image line in accordance with the calculated time interval and the calculated synchronisation delay. Step 3906 may comprise performing n samples of a first image line with the time interval and, after a time period equal to the synchronisation delay, performing n samples of a second image line with the time interval.

Figure 14 is a schematic indicating the delays between transmissions events for the line- by-line imaging sequence for a 1000 Hz shear wave signal. The continuous thinner line 1404 indicates the harmonic shear wave signal at fixed location. A regular short delay, 6t, is used to sample a section (bold line 1402) of the shear wave oscillations at regular intervals (6x shear wave frequency). A longer time delay, At, is used between the last transmission of the current A-line, and the first transmission event of the next A-line, so phase coherence is maintained between A-lines. The value of At is automatically calculated based on the number of samples to be acquired, the value of 6t, and the shear wave frequency.

A benefit of imaging preclinical tumours using the L22-14vX probe is improved image quality. The line-by-line approach works because the vibration source is continuously on, providing a "steady-state" shear wave field (as opposed to transient shear wave elastography) meaning we can reconstruct the shear wave field slowly one line at a time (relative to planar wave imaging) because the shear wave field should not be changing.

A method of synchronised line-by-line processing 3800 is described in relation to Figure 38. At step 3802, the number of samples to acquire per line and the number of samples to be acquired per vibration cycle is selected. At step 3804, the regular short delay is computed using vibration frequency and the number of samples to be acquired per vibration cycle. At step 3806, the number of cycles to be sampled is determined using the number of samples to acquire per line and the number of samples to be acquired per vibration cycle. At step 3808, the synchronisation delay required to maintain coherence between lines is computed using the number of cycles to be sampled, the vibration frequency, the number of samples to acquire per line and the regular short delay.

An example illustrating how the line-by-line approach works is described below with reference to Figure 14 and method 3800:

1) Select number of samples to acquire per line (n = 16), and number of samples to be acquired per vibration cycle (s = 6) (step 3802).

2) With vibration frequency f (1000 Hz), compute (step 3804): i l l

St = — - = - = - = 0.0001667 seconds f x s 1000 x 6 6000

3) Determine value of n / s (16 / 6 = 2.67) (step 3806), this represents the number of cycles to be sampled (bold line 1402), this will usually be a decimal value, and is rounded up to the next integer value, we call this N (2.67 rounds up to 3). This determines the time point when the first sample of the next line should occur, as can be seen in Figure 14 the first triangle appears exactly 3 cycles after the first sample of the current line (first circle). 4) Compute At (step 3808), this is the synchronisation delay required to maintain coherence between lines. It is calculated by computing the time required for N cycles, and subtracting from this the time point of the last sample of the current line: 0.0005 seconds

5) the Verasonics hardware is programmed with values of n, 6t and At so that the transmission and receive events for each line are evenly sampled at 6t intervals, and the typically longer delay At is programmed to allow the start of the next line to be in phase with the previous line. The end result is the hardware performs n repeats of line 1 with interval 6t, followed by a longer delay At before it repeats n repeats of line 2 (at interval 6t), followed by At delay before n repeats for line 3 etc. This process repeats for 128 lines. All the data is saved and processed.

The value of n may be determined when starting up the software and may not be changed thereafter. The value of f and s may be changed interactively on the software whilst imaging. The software automatically updates the required 6t and At values. This allows easy switching of the frequency and correct real time visualisation of the shear wave phase maps. The value of s can be changed to acquire data more/less quickly.

In 3D imaging embodiments, synchronisation with the vibration source signal is required. In 3D synchronisation is required between the multiple planes. This is achieved by using the "sync" signal from the signal generator driving the vibration. This is not required for simple 2D imaging (described above). Examples of 3D data are presented in the appended paper titled "Preclinical three-dimensional vibrational shear wave elastography for mapping of tumour biomechanical properties in vivo".

Various modifications whether by way of addition, deletion, or substitution of features may be made to above described embodiment to provide further embodiments, any and all of which are intended to be encompassed by the appended claims.

Three unpublished papers written by the inventors are included below and form part of the description of the embodiments, thus forming part of the patent application. These papers describe experiments which use and explain embodiments of the present invention in detail. The Use of Travelling Shear Waves in Preclinical Vibrational Shear Wave Elastography

J. Civale, JC Bamber, EJ Harris

Ultrasound shear wave elastography is currently widely used in the clinic and in the laboratory to image stiffness properties of tissue. In preclinical cancer research shear wave elastography can be used to image and monitor tissue changes in xenograft tumour models in response to treatment. The Implementation of preclinical shear wave elastography techniques in xenograft models is challenging due to the small size (<lcm) of the tumour requiring imaging with a high degree of spatial resolution. Using an external source of vibration the delivery of an appropriate mechanical stimulus is also challenging due to the limited access to the tumour. In this paper we present vibrational shear wave elastography, a method which makes use of monochromatic shear waves at frequencies higher than those conventionally used in the clinic to generate travelling waves in tumour. We present and fully describe an ultrasound shear wave imaging system and shear wave speed estimator based on computation of the spatial shear wave autocorrelation function, the latter enabling measurement of the local directionality of the travelling shear wave. We characterise the response of the shear wave speed estimator in a phantom in terms of kernel size, maximum spatial lag, angular interval and spatial median filtering kernel size. We demonstrate the ability of the shear wave speed estimator in to distinguish a stiffer inclusion (3.4 m/s) compared to the background (2.2 m/s). We further demonstrate reliable shear wave speed estimates can be obtained using travelling waves with small kernel sizes, as opposed to a reverberant shear wave method which is more readily implemented using lower vibration frequencies and benefits from kernel sizes larger than the shear wavelength. Finally we demonstrate the feasibility of our vibrational shear wave elastography system in an ex-vivo tumour xenograft, demonstrating the ability to detect travelling shear waves at 1000Hz.

Introduction

Preclinical xenograft tumour models provide a vital tool in the field of oncology allowing researchers to test and validate new cancer treatments. Preclinical models also offer the ability to test imaging techniques capable of detecting tissue biomarkers to assess tumour response to treatment. Ultrasound (US) imaging offers the ability to monitor a number of tissue biomarkers, including for example mechanical properties where a measure of shear wave phase velocity may be used to infer the Young's modulus. In recent years there has been substantial interest and development in shear wave elastography methods that can be used in the clinic for improved diagnosis and treatment monitoring. A number of commercial devices are available for use in the clinic [Evans et al, Cosgrove et al, Chang et al], state of the art systems make use of relative long duration 'pushing' beams that take advantage of the acoustic radiation force impulse (ARFI) to generate localised displacement in the tissue. The mechanical impulse generates transient shear waves [Nightingale] which travel at relatively slow speeds (<10 m/s) and can be imaged using fast ultrasound imaging sequences. Determination of the shear wave speed enables a quantitative assessment of local tissue stiffness, however the presence of vertical banding artefacts in ARFI based shear wave speed image has been previously described [Harris, Berg]. More recently there has been renewed interest in developing shear wave methods where the mechanical impulse is provided by external sources. Examples include reverberant methods [Parker, Ormachea] where the aim is to generate a shear wave field that can be considered as the superposition of random, isotropically distributed plane waves. The reverberant field can be generated in principle by having multiple vibration sources, where reflections of shear waves from boundaries are not considered problematic because they contribute towards generating the reverberant field. Estimation of the shear wave speed is performed by measuring the spatial shear wave autocorrelation function over a finite region of interest. Whilst early clinical results in the clinic are encouraging, the reverberant method relies on a large kernel size, typically larger than the shear wavelength in order to provide meaningful estimates of shear wave speed. A further example of shear wave elastography using multiple vibration sources is a timereversal based method where a cross-correlation of shear wave signals is used to derive the necessary delays required to focus the shear wave field at a region of interest [Benech et al, Brum et al]. This method requires a learning stage where shear wave signals are detected from multiple vibration sources which are activated in turn and individually. The time reversal is then applied in software, and shear wave speed can be estimated based on a number of methods [Brum et al 2021].

Whilst it may be possible to generate a reverberant shear wave field in a xenograft tumour using a single vibration source, realistically this can only be achieved using relatively low frequencies. The shear wave speeds in tumour tissue of 2.5-3.5 m/s [Li, Li Cancer Res, Riegler] is sufficiently high that the shear wavelength at low frequencies (<200 Hz) is typically larger than the dimension of the tumour, rendering the reverberant technique unsuitable for mapping tissue stiffness heterogeneity. Moving to higher frequencies (500 Hz and above) the generation of a reverberant shear wave field becomes difficult due to shear wave attenuation which rapidly eliminates shear wave reflections which contribute towards the reverberant field. The use of multiple vibration sources could potentially overcome these limitations, however the small size of the tumours, and the need to position the imaging probe in close contact with it, effectively reduces viable access positions rendering the reverberant and time-reversal methods practically unfeasible. For these reasons we propose the use of vibrational shear wave elastography in preclinical tumours using travelling shear waves at frequencies up to and potentially above 1000 Hz. The required shear waves can be generated using a single vibration source. We introduce a shear wave speed estimator based on experimental evaluation of the localised shear wave field spatial autocorrelation function. The direction of the travelling shear wave, a property that cannot be assumed a-priori with arbitrary positioning of the vibration source, can be determined from the autocorrelation function. We fully characterise the shear wave speed estimator in a phantom in terms of kernel size, maximum spatial lag, angular interval and spatial median prefiltering kernel size. Application of vibrational shear wave elastography with travelling shear waves is then demonstrated in xenograft tumour ex- vivo. More generally, the choice to use the spatial autocorrelation function as the basis of the measurement means alternative estimates of shear wave speed may be possible if simple travelling wave conditions are not met.

Materials and Methods

Shear wave generation

Vibrations were generated by the use of a mini shaker (model 4810, Bruel & Kjaer, Denmark) coupled to the test samples via a carbon fibre rod and Delrin® contactor applied to the skin. The shaker was driven by a narrowband continuous signal (variable frequency between 200Hz to 1200 Hz) generated by a signal generator (Agilent 33120A) amplified by a variable gain 500W audio amplifier (Intimidation VLV-1000, UK).

US imaging system and shear wave detection

Ultrasound imaging was performed using a Vantage 256 imaging system (Verasonics Inc., Kirkland, WA, USA). For breast elastography phantom validation measurements an L7-4 imaging probe with 5.0 MHz centre frequency was used, for preclinical use imaging small (<1.0cm) tumours in a mouse model a L22-14vX higher frequency probe (18 MHz centre frequency) was used. Imaging data with the L7-4 probe was acquired using a planar wave imaging sequence, whereas with the L22-14vX probe a more conventional synchronised sequential line-by-line scanning sequence was used to image the ex-vivo tumour [Civale et al 2022 under review]. Shear waves were detected from image In-phase and Quadrature (IQ) data reconstructed by the Verasonics system. IQ data consisted of a 3D array with dimensions of 256 axial samples acquired at the equivalent of A/2 spacing, 128 samples in the lateral direction with separation equivalent to the probe elements spacing. For the phantom measurements with the L7-4 probe, 200 temporal samples were acquired at 7000 Hz pulse repetition frequency, for the tumour data obtained with the L22-14vX probe 52 temporal samples were obtained at a 6000 Hz pulse repetition frequency. Shear waves were detected using a 2D phase method [Loupas et al 1995] with axial and temporal averaging parameters M and N set to 5 and 3 respectively applied directly to the IQ data. Spectral filtering was performed by quantifying the vibration frequency component from a fast Fourier transform in the temporal domain, yielding a 2D amplitude and phase representation of the shear wave field. Spatial filtering was performed prior to shear wave speed estimation using a rectangular kernel median filter applied to the real and imaginary components of the shear wave field data (Figure 2).

Shear Wave Speed Estimation

The complex autocorrelation (AC) function A, with spatial lags Ax and Az may be computed from a complex shear wave field p(x, z, t) where x and z represent two spatial dimensions, and t represents time:

A(Ax,Az, t) = p*(x, z, )p(x + Ax,z + Az, t) (Eq. 1) where p* represents the complex conjugate of p.

We consider the case where the shear wave field is in a steady state, such as for example with a monochromatic and continuous source. Furthermore if we assume the shear wave field in a small region of interest can be reduce to being described as a travelling planar wave, the shear wave field p may be quantified as: p(x, z.t) = p _o e^ ^kxX+kzZ ~ ^wt+<p ^ (Eq. 2) where the shear wave wavevector k has components k _x and k _z in the lateral and axial directions respectively, a> is the angular frequency, represents a constant phase offset, and po is a constant shear wave amplitude. If the expression for p is substituted into equation 2, the dependence of A on t vanishes. If the real part of A is taken, we obtain:

‘ [A _xz Ax,Az')] = p cos k _x Ax + k _zA z) (Eq. 3) a result which shows how the AC function is independent of x and z for the special case of a non-lossy planar shear wave. The AC function may be reformulated following transformation of the cosine argument variables to polar coordinates r and 0: k _x = k ■ n _x = k cos(0 _fe ) (Eq . 4) k _z = k n _z = k sin(0 _fe ) (Eq . 5)

4x = r cos(0) (Eq . 6)

4z = r sin(0) (Eq . 7) where n _x and n _z represent unit vectors in the x and z directions respectively, and k is the magnitude of the wavevector k, resulting in:

9l[4 _%z (r, 0)] = PQ cos(/cr cos(0 _fe ) cos(0) + kr sin(0 _fe ) sin(0)) (Eq . 8)

In the special case where the autocorrelation function is sampled linearly where dk = d, the autocorrelation function is sampled through the origin and in the same orientation of the wavevector k, a simple cosine function with minimised width independent of 0 is obtained:

^ r)] = Po cos(fcr) (Eq. 9)

Experimental estimates of k were obtained from small regions of interest (ROI) over which planar shear wave conditions were assumed. The experimentally measured complex shear wave field p is quantified over the ultrasound imaging scan plane in a 2D grid with lateral and axial coordinates x, and z ₇ . Prior to evaluation of the AC function the magnitude of p is set homogeneously to unity [Ormachea et al 2019], rendering shear wave speed estimation dependent on the phase data only. The normalised autocorrelation function, with cartesian spatial lags expressed as integer multiples / and v of the pixel spacings Ax and Az respectively, was evaluated over a (2/V+l) x 2M+1) grid defining a ROI centred on pixel with coordinates x, and zj, up to a maximum spatial lag r _max

The real part of A was resampled in polar coordinates system with radius, r, and polar angle, 0, using bilinear interpolation [William et al] (Figure 3) with intervals Ar and Ad respectively. The interval Ar was set equal to the IQ data axial sampling. The default polar angle interval Ad was set automatically so that the 2D cartesian autocorrelation function was sampled homogeneously along the outer edges by intervals equivalent to the IQ data axial sampling interval. The value of dk was subsequently determined by evaluating isocontour at a user defined value (i.e., 0.7), and finding the coordinate of 0 which minimised the corresponding r coordinate (see Figure 3). With planar waves the isocontour plot results in two local minima separated by n over the full angular range, therefore to avoid ambiguity and maintain 9k aligned with the correct orientation, the slope of the autocorrelation phase data was inspected and where necessary the 9k value was reverted (by adding or subtracting n, as appropriate). The real part of A at polar angles 9k and 9k±n, were fitted to a cosine function, yielding the wavenumber k _] . Finally, a shear wave speed c was obtained

Test phantom validation and preclinical tumour testing

The performance of the imaging system and shear wave speed estimator using travelling shear waves was tested on a CIRS breast elastography phantom, (model 059, CIRS inc, Norfolk, VA, USA) which included regions of stiffness contrast. According to the manufacturers specification the Young's modulus of the material in the background is 10- 15 kPa, and approximately double this value in the inclusions. The L7-4 imaging probe was positioned so that one of the stiff inclusions was visible in the B-mode image. The shaker providing the vibration was placed in contact with the surface of the phantom with the Delrin® contactor (1mm radius spherical shaped tip) at a position aligned with the imaging plane. A 700 Hz vibration frequency was used to generate a travelling shear wave generated at the surface of the phantom. The effect of varying the autocorrelation kernel size (3 - 15 mm), the maximum spatial lag r _max (1 - 5 mm), polar angle interval A9 (10° - 90°), and median filtering kernel size (0 - 2.5 mm) on the shear wave speed estimates were quantified. Goodness of fit, measured as the mean normalised residual was quantified for the cosine fit used to estimate the shear wave speed. Subsequently a more complex shear wave field was generated using 4 contactors placed around the CIRS phantom with a 500 Hz vibration frequency, with the aim of generating a reverberant shear wave field to compare the performance of the shear wave speed estimator with the method previously published using reverberant methods [Ormachea et al].

The performance of the shear wave speed estimator using high frequency vibrational shear wave elastography was tested in an ex-vivo xenograft tumour model grown from C33A cell line in a nude athymic mouse. All animal experiments were approved by The Institute of Cancer Research Animal Welfare and Ethical Review Body, and performed in accordance with the UK Home Office Animals (Scientific Procedures) Act 1986, the United Kingdom National Cancer Research Institute guidelines for the welfare of animals in cancer research and reported according to the ARRIVE (animal research: reporting in vivo experiments) guidelines. Results

Figure 4 illustrates use of the directional detector in computing shear wave speed images indicating the stiff inclusion in the CIRS phantom. The inclusion is visible under B-mode imaging (Figure 4a), and its effect on the travelling shear wave wavelength is visible in the phase image (Figure 4b) with an apparent elongation of the shear wavelength, consistent with an increases shear wave speed (Figure 4c) associated with the stiff inclusion. The detected shear wave direction (0 _mzn ) at each measurement location is shown in Figure 4d, indicating the shear waves propagating away from the contactor position for the majority of the field of view. The effect of varying the polar angle interval is shown in Figure 5, where three positions of interest in the phantom were defined: stiff inclusion (A), background (B), and a region with evidence of shear wave interference or reverberation (C). A large interval polar angle will produce shear wave speed over-estimates due to poor alignment with the true wave direction. As A0 is reduced, more precise dk images are obtained thus reducing the extent of over-estimation in shear wave speed. In the 700 Hz example the shear wave speed estimates were found to be stable with an interval set to 20° or lower. The effect of varying the maximum spatial lag distance r _max , is shown in Figure 6. A small r _max covers a relatively small portion of the cosine curve to estimate the appropriate wavenumber, whereas a larger r _max value allows a much larger portion of the cosine function to be fitted, as shown by the curves in the bottom row of Figure 6. In our example the smallest r _max value produced smaller shear wave speed estimates compared to the larger r _max values for the locations of interest in the phantom. The normalised residual plots (Figure 6, middle row) increased as r _max was raised, particularly for the position with shear wave interference (C) where a visibly poor fit is obtained using r _max of 2.0 mm leading to an increased estimate in the shear wave speed. These results therefore indicate how the choice of r _max can affect shear wave speed estimates in locations where detection of the shear wave is poor, or where the nature of the shear wave field does not conform to the assumption of a plane travelling wave.

The choice of kernel size is important in determining the trade-off between spatial resolution and shear wave speed image contrast. In Figure 7 the effect of varying the kernel size from 2 to 10 mm on the shear wave speed images is shown. Contrast to noise ratio was calculated based on a circular region of interest centred in the stiff inclusion, and an annulus covering the background as indicated in Figure 7. A kernel size of 4 mm was found to be produce the best overall contrast to noise ratio (Figure 7, bottom right), with a reduction in contrast attributable to a reduction in the shear wave speed difference between the inclusion and the background as the kernel size is increased. The mean shear wave speed values using a 4mm kernel in the inclusion (3.3±0.3 m/s) and background (2.3±0.2 m/s) show excellent agreement with the upper end of the range of values expected assuming a density of 1000kg/m ³ and a perfectly elastic, isotropic phantom material (3.2, and 2.2 m/s respectively).

The median filtering that is applied to the detected shear wave field plays an important role of removing high frequency noise prior to computation of the autocorrelation function which would otherwise be compromised. The effect of varying the median filtering kernel size is shown in Figure 8 for a number of examples in terms of the phase, shear wave speed, and normalised residual images. Autocorrelation profiles and corresponding cosine fits are included for the locations of interest. These profiles indicate the poor fit obtained when no median filtering is applied, with relatively large normalised residuals. An improved fit was obtained with a median filtering kernel size >0.7mm for the stiff inclusion and background regions. For location C (right-most location in Figure 8) corresponding to a position more distant from the contactor, and showing evidence of shear wave reflection, the choice of kernel size had a large influence on the resulting autocorrelation profile, for example at the largest kernel size tested of 2.5mm the shear wave speed was almost equal to the estimates in the stiff inclusion region. These results therefore indicate the choice of spatial median filtering kernel size may impose a degree of bias in shear wave speed estimates depending on the signal to noise ratio in shear wave detection, or relative conformance to the assumption of plane progressive shear waves.

Having established the dependence of shear wave speed estimates on the individual parameters, further insight into the relative accuracy of the cosine fit can be obtained by inspecting the linear extracted AC data at 6 _m!n with the corresponding fit for different kernel and r _max sizes, as shown in Figure 9. The key feature of these plots is a better agreement of the AC data with the fitted cosine fit using smaller kernel sizes for all three locations in the phantom. Furthermore, less sensitivity to r _max is also shown at smaller kernel sizes, indicating the travelling shear wave is suitable to shear wave speed measurements with smaller kernel sizes and improved spatial resolution. In the case of a reverberant field, generated in the phantom by using 4 vibration sources randomly placed around the surface, and a reduction in the vibration frequency to 500 Hz, a location in the centre of the phantom away from any apparent inclusions was chosen and analysed. The resulting phase and 2D autocorrelation images for a number of kernel sizes are shown in Figure 10. With a small kernel size of 3 mm, approximately similar in size to the shear wave wavelength, the 2D autocorrelation function has some features similar to those that would be expected from those of a plane shear wave, namely a banded appearance similar to the one shown in Figure 2a. As the kernel size is increased these features are lost and the 2D autocorrelation image shows a profile similar to the one expected in a reverberant field [Parker et al, Zvietcovich et al]. Figure 11 shows the AC profiles along S^ and the axial and lateral profiles from the AC data in Figure 10, with the corresponding cosine fits (left column) and theoretically derived autocorrelation function that would be expected under fully reverberant conditions [Zvietcovich] (middle and right columns). The effect of varying the r _max value is included by plotting different series. It can be observed that shear wave speed estimates obtained using the cosine fit (traveling shear wave) is both larger than the previous measurements (3.4-4.5 m/s vs ~2.1 m/s), and highly variable dependent on the choice of kernel and r _max size, where it is clear the cosine function is a poor fit to these data. The autocorrelation data in the axial and lateral directions do however provide better fits to the expected functions for a reverberant field, and in particular for the lateral direction a shear wave speed is obtained (2.0 - 2.1 m/s) that is in good agreement with the values measured in the background region for the 700 Hz travelling shear wave data show previously (Figure 5, 6, 8) using a single vibration source aligned to the imaging plane. Furthermore improved visual fits to the theoretical reverberant profiles are obtained with the use of larger kernels, in contrast to the previous example at 700 Hz with a single vibration source (Figure 9).

Applying vibrational shear wave elastography with a single shaker at a frequency of 1000 Hz to an ex-vivo xenograft tumour in a mouse model produced a shear wave field with wavelength in the order of approximately 5mm (Figure 12b). Despite the large shear wave attenuation it was possible to map the shear wave field through the majority of the tumour section appearing in the scan plane. An example of the shear wave speed image and shear wave direction (0 _mzn ) are given in Figures 12c and 12d respectively using a kernel size of 2mm in this instance with the preclinical L22-14vX probe.

Discussion

We have presented an ultrasound vibrational shear wave imaging system and a shear wave speed estimator based on the propagation of travelling shear waves. Whilst the approach to quantify the shear wave field and compute the AC function over a finite sized ROI is similar of the work of Parker et al [Parker et al], our method differs fundamentally in that we use travelling shear waves as opposed to generating a reverberant shear wave field. The travelling shear wave approach is more favourable in circumstances where generating a suitable reverberant field is difficult. For example in small xenograft tumours a viable way to generate a reverberant field is to use a sufficiently low vibration frequency (<300 Hz) using a single source to ensure shear waves can effectively reverberate around the tumour, however the wavelength (>10 mm) associated with the these frequencies is typically larger than the size of the tumour, hence it is not possible to measure shear wave speed on the scale of a few millimetres. To overcome this problem higher vibration frequencies may be used, however higher shear wave attenuation impedes effective reverberant conditions, resulting in attenuated travelling shear waves. The use of multiple vibration sources to overcome the effects of attenuation is both costly and practically challenging due to the limited access to the tumour. For these reasons using a single vibration source to generate a travelling shear wave offers a simpler setup, and the opportunity to measure shear wave speed at higher frequencies with wavelengths similar to or below the dimension of the tumour. Our shear wave speed estimator can determine the direction of shear wave propagation directly from the AC function, this is necessary because the source of the vibration and the interaction of the shear wave with the tumour tissue means no prior assumption can be made regarding the direction of the shear wave, unlike for example in AR.FI based transient shear wave elastography methods where the axial 'push' beam is expected to generate shear waves propagating laterally.

Our method bears some similarity with other studies previously published which make use of external vibration sources without necessarily generating reverberant fields, for example clinical studies using time-harmonic shear elastography [Tzschatzsch et al 2014, Tzschatzsch et al 2016 UMB, Zhao et al]. One difference between our work and these previously published methods, aside from the lower vibration frequencies used in the clinic, is our estimator is designed to explicitly determine the direction of a travelling shear wave before estimating shear wave speed, instead of directional filtering the shear wave field [Manduca et al 2003] prior to computation of a weighted mean of shear wave speed estimates [Tzschatzsch et al 2016 Medical image analysis, Zhao et al].

Here we have validated the performance of a vibrational shear wave imaging system and shear wave speed estimator based on high frequency travelling shear waves in a CIRS breast elastography phantom. Provided the underlying assumptions of local plane wave conditions aligned with the imaging plane are maintained, the technique is reliable and able to detect changes in shear wave speed as demonstrated in the phantom. An important step in the data processing sequence was found to be the size of the 2D spatial median filtering kernel. Without median filtering the autocorrelation function that is computed can be considered as resulting from the combination of the shear wave field data and any noise source. If the latter dominates, for example if no shear wave field is present, or the ultrasound echo signal is very noisy, the autocorrelation function will become highly dependent on the form of the noise signal, manifesting as a sharp spike centred around zero spatial lag if the nature of the noise signal is uncorrelated between adjacent image pixels. Median filtering is a simple and effective way to reduce the effects of noisy shear wave data, aiding recovery of autocorrelation function data that provide an improved fit to the cosine function, as shown in Figure 8. Alternative method to reduces noise in shear wave data may include spatial finite impulse response filters [Ormachea et al], k-space or spatial frequency domain filtering, or other data acquisition strategies aimed at improving shear wave signal to noise ratio by increasing the duration of the acquired shear wave signal, or for example by employing suitably angular compounding methods [Denarie et al]. The spatial median filter investigated in this study were maintained below the size of shear wave wavelength, where larger dimensions than this would result in loss of the shear wave phase data necessary to correctly estimate shear wave speed.

The choice of r _max , the maximum spatial lag over which to fit the cosine function was found to potentially introduce bias to shear wave speed estimates. A short r _max biases the result towards any remaining noise following the median spatial filtering step. A relatively large r _max can introduce errors because the autocorrelation can eventually depart from the cosinusoidal fit where strict plane shear wave conditions are not met. Therefore it is likely there is an optimum choice of r _max in between the two conditions above, covering the main peak of the cosine function down to the first zero-crossing. A refinement in estimating shear wave speed could then be applied at the cosine fitting stage by limiting the fit to AC values above a pre-defined threshold. Under reverberant shear wave conditions, where the assumption of a travelling shear wave is no longer tenable, the choice of r _max can significantly impact shear wave speed estimate as shown in the left column of Figure 11.

The choice of optimal kernel size be used is likely to depend on the trade-off between shear wave speed image contrast-to-noise ratio and spatial resolution. This optimum choice in this trade-off is likely to be dependent on the signal to noise ratio of the imaged shear wave field data, where more spatial averaging may be required to reduce uncertainty in shear wave speed estimates. In the phantom study we have demonstrated increasing contrast to noise ratio with increasing kernel size, up to the diameter of the stiff inclusion. An interesting finding from a careful analysis of the linear AC profiles (black series in Figure 9) was that an improved fit with the cosine function was obtained with the smaller kernel sizes. This result can be understood in terms of the size of the kernel relative to the shear wavelength, when the former is small plane wave conditions are largely maintained over the locality defined by the kernel size. When the kernel is relatively large it will encompass a greater section of the shear wave field, which will inevitably display variations in the direction of shear waves due to diffraction effects with the net result of the AC function not conforming to the strict cosinusoidal profile at larger spatial lags. The behaviour of experimentally determined AC profiles under reverberant field conditions was different: despite little variability in shear wave speed estimates with different kernel sizes, improved conformance to the axial and lateral AC profiles was obtained with larger kernel sizes (Figure 11). This finding can be understood in terms of the random, spatially variant nature of a reverberant field. A larger kernel size will naturally cover more of this variability (Figure 10) and hence the AC function form which it derives will more closely resemble the properties expected from a theoretically fully reverberant field.

An important limitation of the study is the requirement of a travelling shear wave with local plane wave conditions. In practice it may not always be possible to eliminate unwanted shear wave reflections, or the propagation of surface shear waves travelling through the imaging plane causing bias in shear wave speed estimates, an example of the latter is possibly visible near the top of the phantom (phase image in Figures 4) resulting in apparent shear wave speed over-estimates (Figure 6). Possible solutions to these problems could be in the form of directional filters aimed at detecting and isolating the principle shear wave components prior to shear wave speed estimation, and the implementation of imaging in orthogonal planes or in 3D. An advantage of estimating the shear wave speed via autocorrelation is that quantitative properties of the autocorrelation function could be defined to assess whether the local shear wave field can be considered consistent with the properties of a travelling wave, or whether another shear wave speed estimation method is more appropriate, for example with reverberant conditions. This type of approach could potentially enable an automated differential shear wave speed estimation, where the estimation method to be used is determined by assessing quantitative parameters associated with the autocorrelation function. The details of how such an approach could be implemented in practice are beyond the scope of this paper and are left to a future study.

Conclusion

The use of a shear wave speed estimator based on high frequency travelling shear waves has been validated in a phantom study, and its feasibility demonstrated in an ex-vivo preclinical tumour. The technique is well suited for situations such as in preclinical studies involving small tumours where simple travelling shear waves can be generated using higher vibration frequencies, offering the potential to map spatial heterogeneity in tissue stiffness properties. References

Berg WA, Cosgrove DO, Dore CJ, Schafer FK, Svensson WE, Hooley RJ, Ohlinger R, Mendelson EB, Balu-Maestro C, Locatelli M, Tourasse C. Shear-wave elastography improves the specificity of breast US: the BE1 multinational study of 939 masses, Radiology. 2012 Feb;262(2) :435-49.

Brum Javier, Benech Nicolas, Gallot Thomas, Negreira Carlos. Shear Wave Elastography Based on Noise Correlation and Time Reversal. Frontiers in Physics, 9, 2021, p50. https://doi.org/10.3389/fphy.2021.617445

Brum J, Catheline S, Benech N, Negreira C. Shear elasticity estimation from surface wave: the time reversal approach. J Acoust Soc Am. 2008 Dec;124(6) :3377-80. doi: 10.1121/1.2998769. PMID: 19206764.

Catheline S, Benech N, Brum J, Negreira C. Time reversal of elastic waves in soft solids. Phys Rev Lett. 2008 Feb 15; 100(6) : 064301. doi: 10.1103/PhysRevLett.100.064301. Epub 2008 Feb 15. PMID: 18352476.

Chang JM, Moon WK, Cho N, Yi A, Koo HR, Han W, Noh DY, Moon HG, Kim SJ. Clinical application of shear wave elastography (SWE) in the diagnosis of benign and malignant breast diseases. Breast cancer research and treatment. 2011 Aug;129(l) :89-97.

Cosgrove DO, Berg WA, Dore CJ, Skyba DM, Henry JP, Gay J, Cohen-Bacrie C. Shear wave elastography for breast masses is highly reproducible. European radiology. 2012 May;22(5) : 1023-32.

Denarie B, Tangen TA, Ekroll IK, Rolim N, Torp H, Bjasta T, Lovstakken L. "Coherent Plane Wave Compounding for Very High Frame Rate Ultrasonography of Rapidly Moving Targets" in IEEE Transactions on Medical Imaging (2013) 32 (7), pp. 1265-1276

Evans A, Whelehan P, Thomson K, McLean D, Brauer K, Purdie C, Jordan L, Baker L, Thompson A. Quantitative shear wave ultrasound elastography: initial experience in solid breast masses. Breast cancer research, 2010 Dec; 12(6) : l-1.

Harris E, Sinnatamby R, O'Flynn E, Kirby AM, Bamber JC. A Cross-Machine Comparison of Shear-Wave Speed Measurements Using 2D Shear-Wave Elastography in the Normal Female Breast. Applied Sciences. 2021; l l(20) :9391. https://doi.org/10.3390/appll209391

Loupas, Thanasis, Powers, and Robert W. Gill. "An axial velocity estimator for ultrasound blood flow imaging, based on a full evaluation of the Doppler equation by means of a two- dimensional autocorrelation approach." IEEE transactions on ultrasonics, ferroelectrics, and frequency control 42.4 (1995) : 672-688.

Manduca A, Lake DS, Kruse SA, Ehman EL. Spatio-temporal directional filtering for improved inversion of MR elastography images. Medical Image Analysis, Volume 7, Issue 4, 2003, Pages 465-473, https://doi.org/10.1016/S1361-8415(03)00038-0.

Nightingale KR, McAleavey SA, Trahey GE. Shear wave generation using acoustic radiation force: in vivo and ex vivo results. Ultrasound Med Biol (2003); 29(2) : 1715-1723

Ormachea J, Castaneda B, Parker KJ. "Shear wave speed estimation using reverberant shear wave fields: implementation and feasibility studies." Ultrasound in medicine & biology 44.5 (2018) : 963-977. Parker KJ, et al. "Reverberant shear wave fields and estimation of tissue properties. Physics in Medicine & Biology 62.3 (2017) : 1046.

Tzschatzsch H, et al 2014 Phys. Med. Biol. 59 1641

Tzschatzsch H, Ipek-Ugay S, Nguyen Trong M, Guo J, Eggers J, Gentz E, Fischer T, Schultz M, Braun J, Sack I. Multifrequency Time-Harmonic Elastography for the Measurement of Liver Viscoelasticity in Large Tissue Windows. Ultrasound in Medicine & Biology, Volume 41, Issue 3, 2015, Pages 724-733. https :/7doi.orq/'10.1016Zi.ultrasmedbio.2014. i l .009.

Tzschatzsch H, Nguyen Trong M, Scheuermann T, Ipek-Ugay S, Fischer T, Schultz M, Braun J, Sack I. Two-Dimensional Time-Harmonic Elastography of the Human Liver and Spleen. Ultrasound in Medicine & Biology. Volume 42, Issue 11, 2016, Pages 2562-2571. https://doi.org/10.1016Zj.ultrasmedbio.2016.07.004.

Tzschatzsch H, Guo J, Dittmann F, Hirsch S, Barnhill E, Jbhrens K, Braun J, Sack I. Tomoelastography by multifrequency wave number recovery from time-harmonic propagating shear waves. Medical Image Analysis, Volume 30, 2016, Pages 1-10 https://doi.Org/10.1016/j.media.2016.01.001.

William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (1992). Numerical recipes in C: the art of scientific computing (2nd ed.). New York, NY, USA: Cambridge University Press, pp. 123-128

Zhao H, Song P, Meixner DD, et al. External vibration multi-directional ultrasound shearwave elastography (EVMUSE): application in liver fibrosis staging. IEEE Trans Med Imaging. 2014;33(ll):2140-2148. doi: 10.1109/TMI.2014.2332542

Zvietcovich F, Pongchalee P, Meemon P et al. Reverberant 3D optical coherence elastography maps the elasticity of individual corneal layers. Nat Commun 10, 4895 (2019). https://doi.orq/10.1038/s41467-019-12803-4

High Frequency Ultrasound Vibrational Shear Wave Elastography for Preclinical Research

J Civale, V Parasaram, JC Bamber, EJ Harris

Division of Radiotherapy and Imaging The Institute of Cancer Research, Sutton, UK

Preclinical evaluation of novel therapies using models of cancer is an important tool in cancer research, where imaging can provide non-invasive tools to characterise the internal structure and function of tumours. The short propagation paths when imaging tumours and organs in small animals allow the use of high frequencies for both ultrasound and shear waves, providing the opportunity for high-resolution shear wave elastography and hence its use for studying the heterogeneity of tissue elasticity, where heterogeneity may be a predictor of tissue response. Here we demonstrate vibrational shear wave elastography using a mechanical actuator to produce high frequency (up to 1000 Hz) shear waves in preclinical tumours, an alternative to the majority of preclinical ultrasound SWE studies where an acoustic radiation force impulse is required to create a relatively low- frequency broad-band shear-wave pulse. We implement vibrational shear wave elastography with a high frequency (18.5 MHz) probe running a focused line-by-line ultrasound imaging sequence which as expected was found to offer improved detection of 1000 Hz shear waves over an ultrafast planar wave imaging sequence in a homogenous tissue-mimicking phantom. We test the vibrational shear wave elastography in an ex-vivo tumour xenograft, demonstrating the ability to detect shear waves up to 10 mm from the contactor position at 1000Hz. By reducing the kernel size used for shear wave speed estimation to 1mm we are able to produce shear wave speed images with spatial resolution of this order. Finally, we present VSWE data from xenograft tumours in vivo, demonstrating the feasibility of the technique in mice under isoflurane sedation. Mean shear wave speeds in the tumours are in good agreements with those reported by previous authors. Characterising the frequency dependence of shear wave speed demonstrates the potential to quantify the viscoelastic properties of tumours in-vivo.

Introduction

Imaging is an attractive tool for monitoring tumour response to cancer therapy as it can map physical, functional or molecular properties of the entire tumour, as opposed to biopsy which only samples small volumes of tissue. Clinically, it can allow non-invasive assessment of tumour response, at an early stage, improving patient management ¹ . Preclinically, xenograft tumour models are a vital tool allowing researchers to test novel cancer therapies and imaging can be used to non-invasively assess the treatment at multiple time points, providing greater insight into the progression (or regression) of disease, whilst improving animal welfare and reducing animal numbers ² . The identification of imaging biomarkers of tumour response that may be used both clinically or preclinically, requires in vivo validation ³ . Here preclinical tumour models, offer the opportunity for direct correlation of imaging signals with histological correlates of tumour biology.

It is now understood that elastic properties or the tumour play an important role in the development of cancer. For example, increased elastic modulus is characteristic of many solid tumours, and a driver of cancer progression ⁴ and as such the elastic modulus of tissue is established as a diagnostic and prognostic biomarker of cancer ⁵ , and there is increasing evidence to suggest these properties may be useful as biomarkers of response to cancer therapy ⁶ .

Ultrasound shear wave elastography (SWE) uses measurements of the characteristics of shear waves propagating through tissue to determine shear wave speed. A characterisation of shear wave speed dispersion may then be used to characterise tissue viscoelastic properties, for example by quantifying the complex shear modulus. Both the frequency and amplitude of the shear waves used are important to consider. Shear waves must have sufficient amplitude to generate tissue displacements that are detectable using ultrasound. The minimum detectable displacement depends on the system imaging characteristics. Higher frequency shear waves (with shorter wavelengths) can improve the spatial resolution of SWE ⁷ ' ⁹ , however, the attenuation of shear wave amplitude increases with frequency limiting the distance from the shear wave source over which oscillations can be detected.

Common methods used to generate shear waves in tissue include the use of an acoustic radiation force impulse (ARFI) or an external vibrational source. ARFI typically generates broad-band shear waves ¹⁰ with properties which depend on the acoustic and geometrical properties of the impulse, but also on the properties of the tissue with which the ARFI beam interacts. The upper limit of ARFI shear wave spectrum in tissue is usually limited to less than 500 Hz ^11,12 . Preclinical applications of ARFI-based elastography include single tracking location shear wave elastography imaging (STL-SWEI) ¹³ in pancreatic tumour liver metastasis, and harmonic motion imaging methods using the acoustic radiation force to generate localised vibrations at a frequency of 50 Hz in liver metastasis ¹⁴ . Shear wave elastography in preclinical studies has been implemented with a high frequency probe (15 MHz) to characterise breast tumour stiffness in a xenograft model ¹⁵ , and subsequently to study the response to anti-angiogenic therapy ¹⁶ . Elastography techniques using external vibration sources include sonoelastography ¹⁷ , where a reduction in vibration amplitude was correlated to increases tissue stiffness in excised prostate tissue ¹⁸ , and shear wave elastography in mammary tissue ¹⁹ , with vibration frequencies restricted to less than 500 Hz.

The use of external vibration sources enables more control of the excitation, with narrowband single or multiple frequencies, of relatively high amplitudes. In the context of preclinical imaging, in which tumours are small (typically < 10 mm) and superficial, there is opportunity to use high frequency shear waves that have been generated by an external source. High spatial resolution SWE is desirable to be able to discern spatial heterogeneity in the elastic properties of tumours as it is well established that tumours are not homogenous, and heterogenous response to therapy drives treatment failure ²⁰ . Despite the small dimensions of preclinical tumours, attenuation of high frequency shear waves presents challenges.

The aim of the work described in this paper was to assess the feasibility of vibrational shear wave elastography (VSWE) for preclinical imaging systems. To image the finer detail of tumours a sufficiently high ultrasound imaging frequency is needed to provide the required spatial resolution. A high vibration frequency is also potentially advantageous in detecting spatial inhomogeneity in tissue mechanical properties. For these reasons we present a VSWE system using a high centre frequency ultrasound imaging probe, with vibration frequencies in excess of 500 Hz using shear waves generated by a single external vibration source. In this paper we begin by comparing VSWE implementations that do not require coherent compounding, namely an unsteered ultrafast planar wave against a custom designed scanned focused beam imaging sequence compatible with a high frequency continuous monochromatic shear wave signal. To assess the penetration depths at which high frequency shear waves can be detected, we quantify and map shear wave signal Conformance, a metric that quantifies the percentage of the detected oscillation which is accounted for by the shear wave (and vibrational source) frequency. Subsequently we use the focused beam sequence to present shear wave data obtained in xenograft tumour models, both ex-vivo and in-vivo, demonstrating the feasibility of VSWE in estimating shear wave speed using high vibrational frequencies of between 500 and 1000 Hz.

Materials and Methods

Test phantom and tumour model VSWE imaging sequences were tested in a breast elastography biopsy phantom (model 059, CIRS, VA, USA). A homogeneous background region of the phantom was identified and imaged using simple planar wave, and focused beam imaging. According to the manufacturer's specification the Young's modulus in the background region of the phantom was 20±5 kPa, equivalent to a shear wave speed of 2.6±0.4 m/s assuming a mass density of 1000 kg/m ³ and an isotropic incompressible perfectly elastic material.

The feasibility of VSWE was tested in xenograft tumour models. Xenograft tumours were grown following subcutaneous injection of cancer cells into the right flank of athymic nude mice (26g weight). An initial ex-vivo feasibility study was performed using an ex-vivo tumour grown from C33A cell line. When the tumour reached a maximum diameter of 10 mm the mouse was terminated and the tumour was immediately imaged using VSWE at four frequencies (500, 600, 800 and 1000Hz). Ex-vivo imaging of the tumour allowed an initial assessment of the potential of the focused beam shear wave imaging sequence without the influence of breathing motion. Subsequently, three separate tumours grown in separate mice injected with MDA-MB-231 cells line were imaged using VSWE at three different frequencies (500, 700 and 1000 Hz) in-vivo to assess performance of VSWE in vivo. Mice were anaesthetised with 2% isoflurane inhalation, and placed on a heated pad to maintain a temperature of 38°C during imaging. All animal experiments were approved by The Institute of Cancer Research Animal Welfare and Ethical Review Body, and performed in accordance with the UK Home Office Animals (Scientific Procedures) Act 1986, the United Kingdom National Cancer Research Institute guidelines for the welfare of animals in cancer research and reported according to the ARRIVE (animal research: reporting in vivo experiments) guidelines.

Shear wave generation

Shear waves were generated in phantom and the tumour by coupling vibrations from a mini shaker (model 4810, Bruel & Kjaer, Denmark) via a 4mm thick cylindrical carbon fibre rod placed in physical contact with the skin overlying the tumour or to the outer surface of the phantom or tissue being imaged (Figure 13). The shaker was supported by a manual precision stage (Thorlabs, UK) allowing precise positioning (<20pm) of the contactor. A minimal amount of pressure was applied at the contactor to maintain physical contact at the surface of the phantom or tumour skin. The carbon fibre rod and contactor were oriented horizontally parallel to and aligned with the ultrasound imaging plane. Due to the relatively large size of the phantom the contactor remained a small distance away (15-20 mm) from the lateral edge of the ultrasound field of view. The nature of the subcutaneous tumours is such that a large enough portion of it protrudes above the skin of the surrounding tissues, allowing positioning of the contactor on the lateral side of the tumour in the gap between the imaging probe and the skin of the tissue surrounding the tumour. Due to the small size of the tumour, it is possible to visualise the contactor in the ultrasound imagine field of view. The shaker was driven by a narrowband continuous signal (up to 1000 Hz) generated by a signal generator (Agilent 33120A, UK) connected to a 500 W audio amplifier (Intimidation VLV-1000, UK) with independent gain control.

US imaging system and shear wave detection

Ultrasound imaging was performed using a Vantage 256 imaging system (Verasonics Inc., Kirkland, WA, USA) in conjunction with a high frequency (18.5 MHz nominal centre frequency) L22-14vX imaging probe (Verasonics Inc.). The performance of the imaging system was initially tested in a CIRS breast elastography phantom, (model 059, CIRS inc, Norfolk, VA, USA). The Vantage system's image reconstruction software was used to compute image In-phase/Quadrature (IQ) data allowing visualisation of B-mode images and computation of shear waves oscillations. Ultrafast imaging with unsteered planar transmitted waves was compared to a conventional scanned focused beam sequence modified for tracking of shear wave oscillations at high frequency. The rationale for this comparison was to determine whether ultrafast imaging with planar wave unfocused transmission was capable of adequately detecting shear wave oscillations, when compared to a focused beam approach. Synthetic focusing methods, for example coherent compounding using multiple steered planar transmissions, were not considered at this time. For the planar wave sequence the imaging depth was set to 12 mm with a pulse repetition frequency set to an exact multiple (lOx) of the shear wave vibration frequency for a total of 200 frames of IQ data. The modified focused beam line-by-line sequence consisted of 52 repeat transmission events for each A-line to sample shear wave oscillations, before moving to the next A-line in the sequence. A total of 128 separate sets of A line transmissions and line reconstructions were performed for an imaging depth of 12 mm. The transmit focus was set to 7mm depth, with a transmit aperture size of 1mm. The focused beam sequence shares some similarity to stroboscopic techniques used to track shear waves ^10,21 . The pulse-echo repetition frequency was set to 6 times the vibration frequency. A longer delay in transmission events was implemented between the last transmission event of an A-line, and the first transmission event of the next A-line to maintain shear wave phase coherence between acquisition of successive A-lines (Figure 14). The duration of the longer delay necessary between A-lines was calculated and set automatically according to the sampling parameters and vibration frequency to maintain phase coherence with the shear wave vibrations at source. A full set of 128 A-line data sets were collected this way, and transferred to the host computer for IQ data reconstruction. Phantom measurements using the different imaging sequences were performed immediately after another maintaining the same contactor position and excitation voltage. A comparison of data acquisition parameters between planar wave imaging and focused line-by-line imaging is given in Table 1. In order to maximise the number of cycles that could be acquired in the line-by-line imaging sequence, the number of samples per cycle (6x) and the Vantage's receive bandwidth were both reduced (50% vs 200%) when compared to the planar wave imaging sequence. These steps were necessary to maximise the temporal duration and hence number of cycles that could be acquired for the line-by-line imaging sequence.

The focused beam line-by-line and planar wave imaging sequences were evaluated in terms of signal to noise ratio. Rectangular and Hanning window transmit apodisation functions were tested for planar wave imaging. To compute signal to noise ratio images, the VSWE sequences were run to obtain a set of IQ data frames without any vibration for each transmission method. One set of data was obtained whilst imaging a homogeneous section of the CIRS breast elastography phantom, and another whilst imaging degassed water in the absence of any scatterers to compute the signal and noise estimates respectively. Signal to noise ratio was quantified in dB using

SNR = 10 logic ( Ps I PN ) where Ps and PN represent the power in the IQ data for the signal and noise measurements. Power in the IQ data sets was computed by multiplying the IQ data with its complex conjugate at each pixel, followed by temporal and local spatial (0.8mm square) averaging.

The reconstructed image IQ data was processed in a similar way as reported by Ormachea et al ²² et al. A 2D phase axial velocity estimator ²³ was applied to the IQ data to compute local tissue axial velocities, amplitude and phase data at the shear wave frequency were obtained from a Fast Fourier transform applied in the temporal domain to the local tissue shear velocity data. The Conformance, C, of the detected shear wave was calculated as the percentage of the total detected vibrational energy that was measured at the shear wave drive frequency: where |S/| represents spectral amplitude components, and |Sf| is the amplitude component at the shear wave frequency. The detection of shear waves depends on the amplitude of the oscillation being large enough to be detected, but also on the signal to noise ratio of the ultrasound echo data. The Conformance measure therefore represents a simple and effective method to visualise detection of shear waves that combines the factors above. A further sanity check is provided by inspecting the shear wave phase images, allowing the user to confirm the shear wave field that is being generated conforms to the expected form, for example shear waves propagating away from the contactor position.

Shear wave speed was estimated from the shear wave phase images using a square shaped kernel based 2D autocorrelation approach similar to that used by Parker et al ²⁴ . Our method differs from the method used by Parker et al ²⁴ in that it assumes non- reverberant, locally-planar shear wave propagation, and the kernel size can be smaller than the shear wavelength, allowing for localised estimation of shear wave speed. The 2D autocorrelation data is used to determine the local direction of the shear wave in the region defined by the kernel. Subsequently, a ID cosine curve is fitted only to the real part of the autocorrelation data aligned with the direction of the shear wave. The relatively small size of the kernel therefore inherently allows for any variation in shear wave propagation through the imaging plane due to diffraction. The cosine fitting process yields a shear wave wavelength estimate, A The wavelength, and vibration frequency f are then used to calculate the shear wave speed, c: c = A f.

Results

B-mode images from a homogeneous section of the CIRS elastography phantom are shown in Figure 15(a-c) for planar wave imaging with rectangular and Hanning transmit apodisations, and with the focused line-by-line sequence. Amplifier gain was maintained constant, and the contactor position was left unchanged between the different imaging scans. The images show inhomogeneous brightness across the image for both planar wave imaging sequences. The focused beam method however shows a higher and more uniform level of brightness across the image. Figure 15(d-f) show the signal to noise ratio measurements for the imaging sequences, quantifying the variations that were observed in the B-mode images. The highest signal to noise ratios were obtained using the focused beam line-by-line imaging sequence which showed a consistent improvement, typically between 10-15 dB, in signal to noise ratio with depth over the planar wave transmission with rectangular apodisation sequence. Applying Hanning windowing to the planar wave transmission improved signal to noise ratio by approximately 5dB. The echogenicity and signal to noise ratio using the focused beam imaging was also found to be more homogeneous across the lateral direction when compared to the planar transmissions.

The CIRS breast phantom was also used to evaluate the detection of 1000 Hz shear waves generated in the phantom. For these measurements the contactor was positioned aligned with, but just outside, of the imaging plane so that shear waves could be detected travelling across the imaging plane. Figure 16 shows a comparison of the phase, Conformance and shear wave speed images obtained using planar wave imaging with Hanning transmit apodization, and the focused line-by-line sequence. The phase images (Figure 16a, b) show the progressive wave nature of the shear waves. The improved ultrasound signal to noise ratio using the focused line-by-line sequence is reflected in the relative appearance of the shear wave phase images, and the greater extent of higher Conformance values (>80%) at positions further away from the contactor. These images illustrate the improved detection of the shear wave oscillations using the focused beam imaging method as an improvement in the apparent penetration depth of the shear wave. The shear wave speed maps show a mostly homogeneous distribution across the phantom with the exception of the most distal region from the probe, and for the planar wave imaging, a region close to the transducer surface with elevated shear wave speed compared to the surroundings (Figure 16e, f). The increased shear wave speed in the distal region can be understood as arising due to the poor signal to noise ratio in this region due to attenuation of the ultrasound echoes. The mean and standard deviation of the shear wave speed images were 2.24±0.21 and 2.32±0.25 m/s for the focused beam and planar wave imaging respectively. If only the middle third section of the phantom only is considered, eliminating the proximal and distal regions to the probe, the mean and standard deviation were found to be 2.15±0.06 and 2.17±0.07 m/s for the focused beam and planar wave imaging respectively. These results are close to the shear wave speed expected (2.3-2.9 m/s) given the manufacturer's specification of a Young's modulus of 20±5 kPa for the phantom material, assuming an elastic isotropic material with a mass density of 1000 kg/m ³ . The results are also in good agreement with the values reported by Ormachea et al ²² on a similar phantom.

The improved performance of the focused beam line-by-line imaging sequence over simple planar wave imaging was not surprising since the latter is also known to be more susceptible to grating lobe artefacts ²⁵ . Vibrational shear wave elastography was therefore tested in the xenograft ex-vivo tumour using the focused beam imaging sequence. An example of the detected shear wave oscillations is given in Figure 17, comparing the time domain signal and spectra from regions with high (>99%) and low (60%) Conformance. Examples of phase and Conformance images obtained in the ex-vivo tumour are shown in Figure 18, for vibration frequencies from 500 Hz to 1000 Hz. With the contactor placed on the tumour skin, as indicated in Figure 18a, it is possible to visualise the penetration of the shear wave into the tumour. Higher amplifier gain settings were required at the higher frequencies. Amplifier gain was increased within a safe operating limit to improve detection of the shear waves, as determined by inspection of the real-time Conformance image. At a frequency of 500 Hz shear wave oscillations are detected throughout the tumour. As the frequency is increased the Conformance images show a progressive reduction in the ability to detect shear wave oscillations due to increased shear wave attenuation. At 600 Hz Conformance is reduced towards the distal edge of the tumour with respect to the position of the contactor. At 800 Hz the shear waves do not propagate further than 8-10 mm from the contactor, and at 1000 Hz it is difficult to detect any shear waves beyond 10 mm from the contactor. The phase images in the tumour show a progressive shear wave similar to the ones observed in the CIRS phantom at the higher frequencies of 800 and 1000 Hz. At the lower frequencies there is more evidence of possible interactions between the shear wave propagating into the tumour reflections from the tumour surface.

The appearance of the tumour in B-mode imaging was used to segment the shear wave speed images to show data in the tumour region only. The shear wave speed images computed using the three different kernel sizes show the appearance of an apparent stiffer region in the upper part of the tumour (Figure 19). The mean shear wave speed in the tumour increases with frequency for all three kernel sizes (Figure 20a). The main effect of varying kernel size is to alter the spatial resolution of the final shear wave speed image. The smaller kernel size (1mm) produces more detailed images but is more susceptible to noise. Moving to a larger kernel size increases the area used for spatial averaging in computing the autocorrelation and is an effective way to reduce some of the variance in the shear wave speed image. A comparison of the mean and standard deviation of the shear wave speeds at different frequencies for the different kernel sizes (Figure 20, b-d) illustrates the reduction of the variance in shear wave speed as the kernel size is increased.

Finally, VSWE was tested in vivo with the focused beam sequence on xenograft tumours in mice under sedation with isoflurane. The B-mode images of three tumours at three vibration frequencies are shown (Figure 21), with their respective Conformance (Figure 22). Outlines of the tumour including a hyper-echoic rim around the edge and the body of the tumour are included. The in vivo Conformance images replicate the ex-vivo findings, namely that at 500 Hz it is possible to generate and detect shear vibrations throughout the tumour, and increasing the frequency reduces the effective penetration of the shear waves within the body of the tumour. At the higher frequency of 1000 Hz shear waves were mostly detected along the top (Figure 22c, 22f and 22i) and bottom of the tumour (Figure 22i). A visualisation of the detected shear wave phase temporal evolution for the three tumours at each frequency can be observed in the supplementary video. These data were used to compute shear wave speed using a 1mm kernel size (Figure 23) which reveals regions of greater shear wave speed consistent with the location of the upper rim, including the skin, of the tumour for mouse 1 and mouse 3 (Figure 23a-c, 23g-i). Analysis of the shear wave speed data obtained within the tumour body (white outline in Figure 23) reveals mean and standard deviations in shear wave speed of 1.9±0.8, 2.5±1.2 and 3.4±1.6 m/s in the three tumours respectively (Table 2). The relatively large variation of shear wave speed within each tumour (up to 50%) is driven by localised regions of increased shear wave speed that are broadly consistent across the three measurement frequencies.

Discussion

To characterise shear wave speed in tumours in preclinical studies using shear wave elastography, high frequencies are desirable for both the ultrasound imaging and for the shear wave vibration source in order to detect and visualise small (<~2mm) features in the tumours. In practice this requires ultrasound imaging frequencies at high frequencies (15 MHz and above), and vibration frequencies approaching 1000 Hz to obtain shear wavelengths of only a few (<~5mm) millimetres. We have demonstrated the feasibility of a vibrational shear wave elastography method implemented with a focused beam, line-by- line, scanning method. Our approach using the conventional focused beam line-by-line method to reconstruct the image and detect the shear waves has some similarity to the work published previously by others ²⁶ ' ²⁸ . Our method differs from those of previous authors in that we take advantage of the continuous nature of the high frequency harmonic vibration which generates a steady state (non-transient) shear wave field. Unlike sequential A-line acquisition used in standard B-mode imaging, here we repeatedly scan individual A lines enabling the use of high pulse repetition frequency to avoid aliasing the high frequency shear waves. Subsequent A-lines transmissions are then synchronised in such a way that phase coherence is maintained at the shear wave frequency, eliminating the need for phase-unwrapping.

The focused beam method was found to produce image data with improved signal to noise ratio (>10 dB) and more homogeneous echogenicity when compared to unsteered planar transmit waves. Line-by-line focused beam imaging was chosen here as an alternative to multi-angle, steered, spatial compounding methods because it represented a simple and efficient way to achieve transmit focusing and maximise image quality. It has been shown by Montaldo et al ²⁵ that images with similar quality to those obtained with focused beam method can be obtained using synthetic focusing by coherent spatial compounding with multiple steered planar transmissions. Indeed, an advantage of coherently adding images obtained with planar wave transmissions over multiple angles is that the overall frame rate compares favourably with the focused beam approach providing an increase of approximately an order of magnitude in overall frame rate. A higher frame rate is a desirable feature in a preclinical ultrasound elastography system where the temporal duration of anaesthesia is an important consideration, however it should not significantly compromise the ability to accurately map shear wave oscillations. A limiting factor of coherent spatial compounding methods is the degree of motion, due to shear wave oscillations or otherwise, that can result in loss of coherence ^{29, 30} . This is likely to be an important consideration when imaging high frequency oscillations (1000 Hz) where in order to adequately track the oscillations a choice has to be made in the trade-off between imaging depth, number of compounded images, and effective frame rate. A study comparing the relative merits of the focused beam imaging method versus synthetic focusing with coherent compounding methods, likely to require motion correction ²⁹ , is beyond the scope of this paper and will be the subject of future studies.

The focused beam approach places greater demands on the Vantage 256 system's memory for image data reconstruction (Table 1). Therefore, to maximise the number of temporal samples of the shear wave oscillations that could be acquired for a single frame, data reduction in the digitised RF signal was achieved by lowering the Vantage receive bandwidth setting (50% vs 200%). The reduced bandwidth did not produce any observable image degradation effects in mapping shear wave oscillations. The number of temporal samples (52) acquired per image line using this method was sufficient to sample up to 8 cycles of the shear wave oscillations with the minimum possible of 6 samples per cycle compatible with the IQ phase change detector ²³ used to track the shear oscillations. This choice proved adequate for monochromatic shear waves as described here, and is potentially sufficient to acquire shear waves signals containing several distinct harmonic frequencies.

The increase in ultrasound echo signal to noise ratio using the focused beam line-by-line, compared to planar wave imaging, is thought to be one of the principle factors in the improvement in shear wave phase images (Figure 16). This improvement was quantified in terms of higher Conformance of the detected shear wave signals, leading to an increase of approximately 5mm in 1000 Hz shear wave penetration depth in the CIRS phantom. The focused beam line-by-line imaging sequence was therefore used to characterise the shear wave fields generated in an ex-vivo xenograft tumour model. Measurements were performed over a range of vibration frequencies to investigate the nature of the shear waves fields that may be generated using a single contactor. At low frequencies (500 Hz) it was possible to detect shear wave oscillations throughout the tumour with good Conformance levels (>80%). As the vibration frequency was increased, the penetration of the shear wave was reduced due to attenuation, and the nature of the shear wave began to more clearly resemble that of an attenuated progressive wave. At 1000 Hz vibration frequency the penetration of the shear wave was limited to within 10 mm from the contactor position, rendering characterisation of the shear wave speed difficult in the regions distal to the contactor. Shear wave speed images computed from acquisitions at different frequencies indicate an overall increase in the shear wave speed with frequency, providing a potential basis to quantify the viscoelastic properties of tissue. An important advantage of using ultrasound to characterise frequency dependent shear wave speed is the relative speed of the measurement compared to other modalities such as for example magnetic resonance elastography ^31,32 . In principle data could be obtained at several frequencies over the course of a few seconds using ultrasound imaging, whereas the long scan times ^5,32 associated with magnetic resonance imaging (> 20 minutes) make this unfeasible for in-vivo applications.

In this study the choice of kernel size had an impact on the level of detail that can be observed in shear wave speed images. A small kernel size allows potential visualisations of small regions (1mm) of the tumour, conversely increasing the kernel size has the effect of smoothing out some of the variance in the shear wave image if a mean value is required. In the example illustrated here, a small region with increased shear wave speed is visible towards the top of the tumour was observed at all frequencies and kernel sizes (Figure 19). In the absence of an independent gold-standard measurement of tissue stiffness this feature roughly coincided with a hyper-echoic region in the B-mode image. For in-vivo studies the choice of kernel size to be used is largely going to be dependent on a compromise between improved spatial resolution, with a small kernel size, versus the need to reduce variance by averaging over a wider area with a larger kernel in the presence of noisy shear wave data.

Our results using the focused beam approach represent an optimisation of imaging quality, using transmit focusing, at the expense of overall scanning speed, which was not a significant factor for phantom work and ex-vivo tissue studies presented here. In vivo tissue motion may lead to artefacts in the accurate tracking of shear wave displacements, resulting in wave speed image artefacts or errors. For example, Ahmed et al ³⁰ used ARFI based STL-SWE to measure shear wave speed in tumours in mice under isoflurane sedation. Their imaging method consisted of coherent compounding of the echoes from three angled planar transmit waves, and they found artefacts in the shear wave speed images which they associated with breathing motion inhalation events. To overcome this problem these authors monitored breathing motion and applied a gating technique to estimate shear wave speed only during "quiet" times of the breathing cycle. To some extent the effect of breathing motion in our in-vivo data can be observed in terms of the modulation of the skin surface apparent in the B-mode images (Figure 21b, 21c, 21d). The amplitude of this motion was however limited to 0.5mm or less, a value smaller than the shear wavelength. With the tumour positioned on the flank of the mouse and imaged from above, the breathing motion was predominantly aligned in a vertical directional parallel to the imaging beam axis, and approximately orthogonal to the direction of shear waves. Under these conditions and with mice under isoflurane sedation, VSWE results using the focused beam sequence were not significantly compromised by tissue motion without the need to apply any breathing motion gating.

The nature of the shear wave fields generated in xenograft tumours is potentially complex. At 500 Hz where shear wave attenuation is expected to be lower there is possibility of interaction between the outgoing shear wave from the contactor and reflections arising from within the body and boundaries of the tumour. The presence of surface waves (Rayleigh waves) near the boundaries of the tumour is also possible. The wave speed associated with surface waves is not quantitatively equal to that of shear waves travelling in the bulk, where surface waves may be approximately 5% lower than the shear wave speed in the bulk of the tumour for a viscoelastic soft tissue ³³ . The penetration of surface waves within the tumour would be expected to be less than one wavelength from the boundary, and therefore will be dependent on both the exact wave speed and vibration frequency. For data presented here the penetration of surface waves could potentially be as large as 5 mm for the combination of shear wave speed and frequency of 3 m/s and 500 Hz respectively. It was not possible to clearly distinguish between surface Rayleigh waves and true shear waves travelling through the body of the tumour. Another factor to consider is the nature of the contact between the tumour and surrounding tissues, for example the skin. The Conformance images also suggest the shear vibrations, particularly at the higher frequency of 1000Hz, travel around the outside of the tumour and through the skin more readily than within the body. This is not surprising as the skin is directly in contact with the external vibration source. These locations were also consistent with increased shear wave speed with respect to the body of the tumour immediately below (Figure 23 and Table 2). Vibration in the skin may be better described as guided waves (Lamb waves), which potentially include partial standing waves or components which travel perpendicular to the imaging plane leading to wave speed overestimates. Shear wave speed estimates in the body of the tumour showed relatively large standard deviations which were driven by localised regions of higher shear wave speed that were broadly consistent across vibration frequencies in each tumour. Regions of low shear wave speed (<2m/s) also contribute towards the variance, and visually correlate with regions of low Conformance.

Despite the above limitations, and considering the range of factors potentially influencing the mechanical properties of xenograft tumours (collagen content, tumour location, extent of tissue necrosis), the shear waves speeds measured in the body of the tumours are in good agreement with values reported by other authors. Li et al ³⁴ used MR elastography to measure the complex shear modulus of xenograft tumours in mice at a vibration frequency of 1000 Hz. These authors also studied MDA-MB-231 tumours, and reported mechanical properties with equivalent shear wave speeds of 3.3±0.2 m/s (orthotopic), and 2.2±0.2 m/s (intracranial). The same authors reported complex shear modulus for other orthotopic tumours including BT-474 (2.9±0.2 m/s), PANCI (3.3±0.2 m/s), and other intracranial tumours from GTML/Trp53 (2.1±0.2 m/s) to U87 MG (2.5±0.2 m/s), and finally subcutaneous SW620 (2.3±0.2 m/s) ³⁵ . Page et al ³⁶ , reported the shear elastic modulus of SA-LIV tumours using MR elastography at 600 Hz (1.3 m/s). Ahmed et al ³⁰ used STL-SWE, an ARFI based ultrasound shear wave imaging modality to measure the shear wave speed of PDAC xenograft tumours (2.3±0.5 m/s) in mice.

There is scope for further optimisation of VSWE by improving the coupling of vibration from the contactor to the tumour tissue, particularly for frequencies of 1000 Hz where Conformance measures in the body of the tumour was low. This may potentially be achieved by optimising the contactor shape design. A measurement of shear vibrations in 3D, using a 2D array probe for example, and an improved understanding of the nature of shear vibrations in the tumour and surrounding tissue, by modelling or other means, are the next steps to improve characterisation of mechanical properties of preclinical tumours.

Conclusion

We have demonstrated the feasibility of high spatial resolution (<5mm) measurements of shear wave speed using vibrational shear wave elastography in an human tumour xenograft model at high vibration frequencies (500-1000 Hz) in vivo. Vibrational shear wave elastography has therefore potential in quantifying frequency dependent shear wave speed, and hence tissue viscoelastic properties.

Tables Table 1. Parameters for the planar wave and line-by-line imaging sequence used in the study for an imaging depth of 12mm and a shear wave frequency of 1000 Hz. Receive bandwidth indicates the Verasonics Vantage system's sampling frequency setting, where 200% is equivalent to 4 samples per wavelength, and 50% to 2 samples for every 2 wavelengths.

Table 2. Summary of mean shear wave speed (±standard deviation) measured in the body and rim of three in-vivo tumours.

References

[1] Brindle, K. New approaches for imaging tumour responses to treatment. Nature Reviews Cancer (2008), 8(2), 94-107.

[2] Workman, P., Aboagye, E. O., Balkwill, F., Balmain, A., Bruder, G., Chaplin, D. J., ... & Eccles, S. A. Guidelines for the welfare and use of animals in cancer research. British journal of cancer (2010); 102(11), 1555-1577.

[3] O'Connor, J. P., Jackson, A., Asselin, M. C., Buckley, D. L., Parker, G. J., 8<. Jayson, G. C. Quantitative imaging biomarkers in the clinical development of targeted therapeutics: current and future perspectives. The lancet oncology (2008); 9(8), 766-776

[4] Kalli M, & Stylianopoulos T. Defining the role of solid stress and matrix stiffness in cancer cell proliferation and metastasis. Frontiers in oncology (2018); 8, 55.

[5] Runel G, Lopez-Ramirez N, Chlasta J, Masse I. Biomechanical Properties of Cancer Cells. Cells (2021); 10(4) :887. https://doi.org/10.3390/cellsl0040887

[6] Li, J., Jamin, Y., Boult, J. K. R., Cummings, C., Waterton, J. C., Ulloa, J., ... 8<. Robinson, S. P. Tumour biomechanical response to the vascular disrupting agent ZD6126 in vivo assessed by magnetic resonance elastography. British journal of cancer (2014); 110(7), 1727-1732.

[7] Palmeri ML, Rouze NC, Wang MH, Ding X, Nightingale KR. Quantifying the impact of shear wavelength and kernel size on shear wave speed estimation, 2010 IEEE International Ultrasonics Symposium (2010), pp. 13-16, doi: 10.1109/ULTSYM.2010.5935798

[8] Rouze NC, Wang MH, Palmeri ML, Nightingale KR. Parameters Affecting the Resolution and Accuracy of 2-D Quantitative Shear Wave Images, IEEE Transactions in Ultrasonic, Ferroelectrics and Frequency Contol (2012); 59(8), 1729 - 1740. [9] Tzschatzsch H, Ipek-Ugay S, Trong MN, Guo J, Eggers J, Gentz E, Fischer T, Schultz M, Braun J, Sack I. Multifrequency Time-Harmonic Elastography for the Measurement of Liver Viscoelasticity in Large Tissue Windows. Ultrasound in Medicine & Biology (2015); 41, 3, 724-733, doi: 10.1016/j.ultrasmedbio.2014.11.009.

[10] Nightingale KR, McAleavey SA, Trahey GE. Shear wave generation using acoustic radiation force: in vivo and ex vivo results. Ultrasound Med Biol (2003); 29(2): 1715-1723

[11] Tanter M, Bercoff J, Athanasiou A, Deffieux T, Gennisson JL, Montaldo G, Muller M, Tardivon A, Fink M. Quantitative Assessment of Breast Lesion Viscoelasticity: Initial Clinical Results Using Supersonic Shear Imaging. Ultrasound in Medicine & Biology (2008); 34 (9): 1373-1386. https://doi.Org/10.1016/j.ultrasmedbio.2008.02.002.

[12] Mitri FG, Urban MW, Fatemi M, Greenleaf JF. Shear wave dispersion ultrasonic vibrometry for measuring prostate shear stiffness and viscosity: an in vitro pilot study. IEEE Trans Biomed Eng. (2011); 58 (2):235-42. doi: 10.1109/TBME.2010.2053928.

[13] Ahmed R, Ye J, Gerber SA, Linehan DC, Doyley MM. Preclinical Imaging Using Single Track Location Shear Wave Elastography: Monitoring the Progression of Murine Pancreatic Tumor Liver Metastasis In Vivo. IEEE Transactions on Medical Imaging (2020) Jul;39(7):2426-2439. DOI: 10.1109/tmi.2020.2971422.

[14] Payen T, Palermo CF, Sastra SA, et al. Elasticity mapping of murine abdominal organs in vivo using harmonic motion imaging (HMI). Physics in Medicine and Biology (2016) Aug;61(15):5741-5754. DOI: 10.1088/0031-9155/61/15/5741.

[15] Chamming's F, Latorre-ossa H, Le Frere-belda MA, Fitoussi V, Quibel T, Assayag F, Fournier LS. Shear wave elastography of tumour growth in a human breast cancer model with pathological correlation. European Radiology (2013) 23(8), 2079-86. doi: http://dx.doi.org/10.1007/s00330-013-2828-8

[16] Chamming's F, Le-Frere-Belda MA, Latorre-Ossa H, Fitoussi V, Redheuil A, Assayag F, Pidial L, Gennisson JL, Tanter M, Cuenod CA, Fournier LS. Supersonic Shear Wave Elastography of Response to Anti-cancer Therapy in a Xenograft Tumor Model. Ultrasound in Medicine & Biology (2016); 42 (4); 924-930. https://doi.Org/10.1016/j.ultrasmedbio.2015.12.001.

[17] Taylor LS, Porter BC, Rubens DJ, Parker KJ. Three-dimensional sonoelastography: principles and practices. Phys. Med. Biol. (2000); 45, 1477-1494.

[18] Taylor LS, Rubens DJ, Porter BC, Wu Z, Baggs RB, di Sant'Agnese PA, ... & Parker KJ. Prostate cancer: three-dimensional sonoelastography for in vitro detection. Radiology (2005), 237(3), 981-985.

[19] Wang, Y., & Insana, M. F. Viscoelastic properties of rodent mammary tumors using ultrasonic shear-wave imaging. Ultrasonic imaging (2013): 35(2), 126-145.

[20] Bedard, P., Hansen, A., Ratain, M. et al. Tumour heterogeneity in the clinic. Nature (2013): 501, 355-364. https://doi.org/10.1038/naturel2627

[21] Gennisson JL, Baldeweck T, Tanter M, Catheline S, Fink M, Sandrin L, Cornillon C, Querleux B. Assessment of elastic parameters of human skin using dynamic elastography, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control (2004) 51(8), pp. 980-989. [22] Ormachea, Juvenal, Benjamin Castaneda, and Kevin J. Parker. "Shear wave speed estimation using reverberant shear wave fields: implementation and feasibility studies." Ultrasound in medicine & biology 44.5 (2018): 963-977.

[23] Loupas, Thanasis, J. T. Powers, and Robert W. Gill. "An axial velocity estimator for ultrasound blood flow imaging, based on a full evaluation of the Doppler equation by means of a two-dimensional autocorrelation approach." IEEE transactions on ultrasonics, ferroelectrics, and frequency control 42.4 (1995): 672-688.

[24] Parker, Kevin J., et al. "Reverberant shear wave fields and estimation of tissue properties." Physics in Medicine & Biology 62.3 (2017): 1046.

[25] Montaldo G, Tanter M, Bercoff J, Benech N, Fink M. "Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography." IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control (2009) 56(3), pp. 489- 506.

[26] Tzschatzsch, Heiko, et al. "In vivo time-harmonic ultrasound elastography of the human brain detects acute cerebral stiffness changes induced by intracranial pressure variations." Scientific reports 8.1 (2018): 1-9.

[27] Rabin C, Benech N. Quantitative breast elastography from B-mode images. Med Phys (2019); 46(7):3001-3012. doi: 10.1002/mp.13537

[28] Huang C, Song P, Mellema DC, et al. Three-dimensional shear wave elastography on conventional ultrasound scanners with external vibration. Phys Med Biol (2020); 65(21):215009. doi: 10.1088/1361-6560/aba5ea

[29] Denarie B, Tangen TA, Ekroll IK, Rolim N, Torp H, Bjasta T, Lovstakken L. "Coherent Plane Wave Compounding for Very High Frame Rate Ultrasonography of Rapidly Moving Targets" in IEEE Transactions on Medical Imaging (2013) 32 (7), pp. 1265-1276

[30] Ahmed R, Ye J, Gerber SA, Linehan DC, Doyley MM. Preclinical Imaging Using Single Track Location Shear Wave Elastography: Monitoring the Progression of Murine Pancreatic Tumor Liver Metastasis In Vivo. IEEE Trans Med Imaging. (2020) Jul;39(7):2426-2439. doi: 10.1109/TMI.2020.2971422

[31] Muthupillai R, Rossman PJ, Lomas DJ, Greenleaf JF, Riederer SJ, Ehman RL. Magnetic resonance imaging of transverse acoustic strain waves. Magnetic Resonance in Medicine (1996); 36:266-274

[32] Schregel, K., Nowicki, M. O., Palotai, M., Nazari, N., Zane, R., Sinkus, R., ... & Patz, S.. Magnetic Resonance Elastography reveals effects of anti-angiogenic glioblastoma treatment on tumor stiffness and captures progression in an orthotopic mouse model. Cancer Imaging (2020); 20: 1-16.

[33] Royston TJ, Dai Z, Chaunsali R, Liu Y, Peng Y, Magin RL. Estimating material viscoelastic properties based on surface wave measurements: a comparison of techniques and modeling assumptions. J Acoust Soc Am. (2011 Dec;130(6):4126-38. doi: 10.1121/1.3655883

[34] Li J, Zormpas-Petridis K, Boult JK, Reeves EL, Heindl A, Vinci M, Lopes F, Cummings C, Springer CJ, Chesler L, Jones C, Bamber JC, Yuan Y, Sinkus R, Jamin Y, Robinson SP. Investigating the Contribution of Collagen to the Tumor Biomechanical Phenotype with Noninvasive Magnetic Resonance Elastography. Cancer Res (2019); 79 (22): 5874-5883. [35] Li J, Jamin Y, Boult JK, Cummings C, Waterton JC, Ulloa J, Sinkus R, Bamber JC, Robinson SP. Tumour biomechanical response to the vascular disrupting agent ZD6126 in vivo assessed by magnetic resonance elastography. Br J Cancer. (2014); 110(7): 1727-32

[36] Page G, Tardieu M, Besret L, Blot L, Lopes J, Sinkus R, Van Beers BE, Garteiser P. Assessing Tumor Mechanics by MR Elastography at Different Strain Levels. J Magn Reson Imaging. (2019) 50(6): 1982-1989.

Preclinical three-dimensional vibrational shear wave elastography for mapping of tumour biomechanical properties in vivo

Vaideesh Parasaram, John Civale, Jeff Bamber, Simon Robinson, Yann Jamin and Emma Harris

Simple Summary: The increased tissue stiffness associated with cancer has been found to be a significant barrier to effective treatment and indicates increased likelihood of cancer progression. Techniques to map tumour stiffness throughout the whole tumour in three- dimensions will assist with preclinical research that aims to understand the relationship between the stiffness, the underlying tumour biology, and the response of tumours to therapy. We have developed an approach that measures the tumour stiffness in murine models of cancer which are commonly used for cancer research. This approach uses a high frequency vibrational source, ultrasound imaging and a three-dimensional analysis, which has advantages in terms spatial resolution and rapid acquisition times. Here we present the first successful demonstration of the non-invasive three-dimensional measurement of tumour stiffness in two preclinical tumour models and its ability to detect a change in tumour stiffness in response to an anti-cancer drug.

Preclinical investigation of the biomechanical properties of tissues and their treatment- induced changes are essential to support drug-discovery, clinical translation of biomarkers of treatment response, and studies of mechanobiology. Here we describe the first use of preclinical 3D elastography to map the shear wave speed (c _s ), which is related to tissue stiffness, in vivo and demonstrate the ability of our a novel 3D vibrational shear wave elastography (3D-VSWE) system to detect tumour response to a therapeutic challenge. We investigate use of one or two vibrational sources at vibrational frequencies of 700, 1000 and 1200 Hz. The within-subject coefficients of variation of our system were found to be excellent, for 700 and 1000 Hz, 5.4% and 6.2%, respectively. The relative change in c _s measured with our 3D-VSWE upon treatment with an anti-vascular therapy ZD6126 in two tumour xenografts reflected changes in tumour necrosis. U-87 MG vehicle vs drug: Ac _s = -24.7 ± 2.5 % vs 7.5 ± 7.1%, (p=0.002) and MDA-MB-231 vehicle vs drug: Ac _s = - 12.3 ± 2.7 % vs 4.5 ± 4.7%, (p=0.02). Unlike previous preclinical elastography, our system enables rapid (<5 minutes) 3D mapping of quantitative tumour viscoelastic properties in vivo, allowing exploration of regional heterogeneity within tumours and speedy recovery of animals from anesthesia so that longitudinal studies (e.g., during tumour growth or following treatment) may be conducted with excellent temporal sampling. 1. Introduction

Mounting evidence supports the role of increased tissue stiffness in malignant transformation, tumour progression and metastasis [1-4]. Mechanical stress associated with rapid tissue growth, compressed vasculature and lymphatics, and extracellular matrix (ECM) structure and rigidity is the major contributor to this phenomenon. Both elevated solid stress and interstitial fluid pressure (IFP), which influence tumour viscoelastic properties, are two major obstacles to efficient drug delivery, and consequently, therapies targeting the stromal component of tumours are being investigated for therapeutic gain [5].

Innovative elastography techniques based on the use of ultrasound, magnetic resonance (MR) and optical imaging are being developed to non-invasively visualise and measure the viscoelastic properties of tissue in preclinical tumour models and cancer patients in vivo. These emerging methods can inform on the underlying tumour microstructure and treatment-induced changes to tumour integrity with direct clinical translation through the development and validation of biomarkers that will enhance decision-making in the oncology clinic including for diagnosis, treatment planning and treatment response assessment [6-16]. They can also be expected to enhance our understanding of tumour mechanobiology and identify key genomic and pathophysiological drivers of tumour stiffness [17]. To support these efforts, there is a need for a dedicated small animal elastography platform that enables quantitative and reproducible mapping of the biomechanical properties (elastic and viscous moduli) to monitor their acute and long-term evolution upon disease progression or treatment. Rapid 3D mapping of these properties will allow exploration of regional heterogeneity within tumours and time-efficient animal procedures so that longitudinal studies (e.g., during tumour growth or following treatment) may be conducted with excellent temporal sampling.

Shear wave elastography (SWE), which is becoming widely used in the clinic as a diagnostic tool, relies on the generation of shear waves and the subsequent rapid measurement of shear wave speed to enable the quantitative assessments of elastic and viscous moduli [18]. Measuring the shear wave speed at multiple shear wave frequencies allows characterisation of tissue viscous modulus through measurement of the shear wave speed dispersion [18]. The measurement of the biomechanical properties of tissue enabled by ultrasound is an attractive imaging option for a high-throughput preclinical shear wave elastography platform. Commercial clinical scanners typically generate shear waves using an acoustic radiation force impulse push delivered inside the tissue. This uses a broadband transient shear wave pulse with a frequency spectrum that varies inside the tissue and depends on factors such as the depth of the push location, the tissue properties at the push location and the distance the shear waves have propagated. How these biases are corrected for by the scanner manufacturer results in variation in shear wave speed measurements for different depths, and between ultrasound SWE systems [19-22]. To overcome these issues, an external continuous vibrational source can be used to generate narrow-band shear waves with controlled frequency and wave amplitude [23]. This is the approach adopted in the development of both clinical [24] and preclinical [25] magnetic resonance elastography and is also being explored for clinical use with clinical ultrasound systems [26].

Here, we describe the development of a three-dimensional preclinical high frequency vibrational shear wave elastography (3D-VSWE) system, employing up to two external mechanical vibrational sources to generate continuous travelling shear waves in tissue, instead of acoustic radiation force impulses. The attenuation of shear waves increases with frequency, however, since mice tumours are small, with propagation distances typically less than 20 mm, there is opportunity to use external mechanical sources of high frequency shear waves, to achieve the spatial resolution required to investigate tumour heterogeneity and to avoid shear wave reflections from tissues beyond the tumour [27]. Because higher frequency shear waves are more rapidly attenuated, a delicate balance exists between using higher frequencies for high spatial resolution and maintaining the propagation of shear waves throughout the tumour. To help achieve this we have developed and validated two real-time metrics of the shear wave field quality within a tumour to allow users to rapidly optimise the arrangement of vibrational sources and their frequencies and amplitude. We demonstrate that our 3D-VSWE system has excellent repeatability and show proof-of-concept of its ability to map regional differences in, and detect treatment-induced changes to, elastic modulus resulting from changes in tumour integrity.

2. Materials and Methods

2.1 Cell culture

Human triple negative luc-MDA-MB-231 LM2-4 breast cancer cells and luc-U87 MG human glioblastoma cells were authenticated by short tandem repeat (STR) profiling and tested negative for mycoplasma infection prior to tumour implantation. Cells were grown under aseptic conditions in Dulbecco's modified Eagle's medium (DMEM) (Invitrogen, Life Technologies) supplemented with 10% (v/v) foetal bovine serum. 2.2 Animals

Adult (5- to 6-week-old) female athymic NCr-Foxnlnu mice (n=29, Charles River, UK) were used. Mice were housed in specific pathogen-free rooms in autoclaved, aseptic microisolator cages with a maximum of 5 animals per cage and allowed access to sterile food and water ad libitum. On the day of inoculation, cells were washed, trypsinised and counted before re-suspending them in a mixture of DMEM :Matrigel (Corning, UK) (1: 1) for injection. 5x106 luc-MDA-MB-231-LM2-4 (referred to as MDA-MB-231 tumours from here on) or luc-U87-MG cells (referred to as U-87 MG tumours from here on) per mouse were injected subcutaneously in the right flank. Tumours were measured twice a week using callipers and were imaged when the tumours reached a volume of ~350 mm3. For imaging, anaesthesia was induced by an intraperitoneal 5ml/kg injection of a combination of fentanyl citrate (0.315mg/ml) plus fluanisone (lOmg/ml) (Hypnorm, Janssen Pharmaceutical, Oxford, UK) and midazolam (5mg/ml) (Roche, Welwyn Garden City, UK) and water (1: 1 :2). They were then secured sideways on a warm platform maintained at 38°C, with tumours facing up.

Repeatability and kernel size optimisation study: Initially a cohort of mice (n = 5) with MDA- MB-231 tumours were imaged for three consecutive days using 3D-VSWE. This was to both test the repeatability of the set-up and measurements, and to use the data to optimise the number of vibrational sources required (one or two shakers) and the kernel size used for shear wave speed calculation.

Therapeutic challenge: To assess the sensitivity of our approach to changes in tissue integrity, we chose to use a therapeutic challenge with ZD6126 (N-acetylcolchinol-O- phosphate), a vascular disrupting agent shown to induce central necrosis in a wide range of tumour models 24 hours after a single dose of 200 mg/kg [28]. ZD6126 was formulated in 20% of 5% sodium carbonate and 80% phosphate-buffered saline and administered intra-peritoneally. Mice with established MDA-MB-231 or U-87 MG were imaged prior to and 24h after treatment with either 200 mg/kg ZD6126 (n=6 and n=4, respectively) or vehicle alone (n=5 and n=4, respectively). One dataset at 1000 Hz using 2 shakers in the repeatability study and one dataset in the therapeutic challenge were excluded due to incomplete collection of data.

2.3 Vibrational Shear Wave Elastography (VSWE) imaging using ultrasound

In our 3D-VSWE system, presented in Figure 24a, a signal generator was used to generate continuous sinusoidal signals of frequencies 700 Hz, 1000 Hz and 1200Hz to drive one or two vibrational sources (mini shakers, model 4810, Bruel & Kjaer, Denmark). Shear waves were generated in tumours by coupling the two mini shakers via carbon fibre rods and custom-made Delrin® contactors placed in physical contact with the skin overlying the tumour. The signal generator (Agilent 33120A, UK) was connected to a 2 channel 500 W audio amplifier (VLV-1000 Intimidation, UK) generated one or two narrowband continuous signals of 700 Hz, 1000 Hz and 1200 Hz to drive the shakers.

A line-by-line focused-beam ultrasound imaging sequence was used to detect shear waves and the resulting shear wave fields were measured by moving a ID array transducer across the tumour in a step and shoot manner to build up the three-dimensional data. A Vantage 256 imaging system (Verasonics Inc., Kirkland, WA, USA) in conjunction with a high frequency (18.5 MHz nominal centre frequency) L22-14vX imaging probe (Verasonics, USA.) was used. The Vantage system's image reconstruction software was used to compute image in-phase and quadrature (IQ) sampled echo data allowing calculation and visualisation of B-mode images (see below). Shear wave oscillations were detected by repeat transmissions along individual A-lines with a pulse repetition frequency equal to 6 times the shear wave vibration frequency.

Radiofrequency data from the ultrasound scans were collected in the form of IQ data for analysis. Time points of scans, frequency of the vibrating shaker and number of shakers used during the scans formed variables for different scenarios of scanning giving rise to a total of 18 datasets per mouse in the repeatability study and 12 datasets per mouse in the therapeutic challenge study. Given the step length of 0.2 mm in Y-direction, a total of 50 to 75 images in ZX plane were acquired for tumours depending on their Y-axis dimension.

See Figure 24 and corresponding brief description.

Shear wave phase maps (Figure 24b) show the phase of the displacement of the tissue caused by propagation of shear waves through the tumour and, in addition to shear wave amplitude maps, can be visualized by the user in real-time to access the quality of the phase data. To generate shear wave phase maps, a 2D ultrasound phase-shift based axial displacement estimator was applied to the IQ data to measure the local tissue axial displacement that was produced by the shear wave during each ultrasound pulse-to-pulse interval [29]. Fast Fourier transform of the tissue displacement versus time data was used to obtained shear wave amplitude and phase at each spatial location at the shear wave drive frequency.

2.4 Shear wave speed estimation Shear wave speed was estimated by calculating the autocorrelation function of the 3D shear wave phase maps in the direction of local shear wave propagation. We employed a cube shaped kernel-based autocorrelation approach similar to that used by Hoyt et al [30]. A directional detector was used to extract ID autocorrelation data aligned with the local shear wave propagation direction. A cosinusoidal fit to the autocorrelation data was used to estimate the shear wave wavelength, A, and shear wave speed, c _s , using the relationship: c _s = A f, where f is the frequency. The dimension of the cuboid, or kernel, of data used to calculate the autocorrelation function directly influences VSWE spatial resolution and the ability to map heterogeneity of c _s .

2.5 Data quality metrics

We devised two quality metrics (Figure 25a) which could be mapped and displayed to the user in real-time to aid with placement of the contactors and the choice of vibrational frequency: i) the Conformance (C), enabled us to evaluate the quality of the detected shear wave field relative to vibrations at other frequencies and noise. C was calculated as the percentage of the total detected vibrational energy that was measured at the shear wave drive frequency: where |Si| represents spectral amplitude components, and |Sf| is the amplitude component at the shear wave frequency. ii) the Goodness of Fit (GoF), a measure of the mean residual error of the cosinusoidal fit to the autocorrelation function. GoF error indicates how well the shear wave phase and magnitude (over a cuboid kernel) conform to a travelling shear wave and provided a measure of the quality of the estimation of shear wavelength.

2.6 High frequency B-mode imaging

High spatial resolution B-mode images of tumours were obtained using a Vevo 770 and a 30 MHz transducer (Fujifilm Visualsonics Inc, Toronto, Canada) prior to the acquisition of VSWE data for the improved alignment of histological sections with imaging data. Both Verasonics and Visualsonics transducers were positioned using custom made holders such that the images from both systems could be easily co-registered.

2.7 Analysis of ultrasound datasets

IQ data obtained for a set of beam lines for a 3D ultrasound scan were first converted into B-mode images in the plane of the linear array for visualisation of tumour using simple log compression of the magnitude of the IQ data. For 3D analysis, B-mode images were imported to Microscopy Image Browser (MIB) for 3D segmentation of tumour outlines [31]. The 3D tumour outline was used to define a grid of measurement points for shear wave speed calculation inside the tumour. Shear wave speed calculation was performed for a kernel of tissue as described above. The size of the kernel and percentage of overlap between the kernels determined the number of points that were placed on the image for analysis. The calculation was repeated creating an (m x n x p) matrix of points across the 3D data set. Median value of all shear wave speeds estimated at all points was considered as the centroid of the c _s probability distribution for a tumour and has been used for data reduction to produce summary results for illustration in graphical form. Metrics that were used for quantitative analysis comprised of C, GoF, c _s and Slope (slope of c _s versus shaker frequency). For the study of repeatability, C, GoF and c _s were calculated at various kernel sizes ranging from 0.5 mm to 3 mm to assess the effect of kernel size on these parameters. The within-subject coefficient, CVWS was calculated as previously described using the equation [32] :

CVWS (%) = 100 x ([E(A/m)2]/2n)l/2, in which m is the mean of the 2 paired repeat determinations of median c _s whose absolute difference is A, and the sum is taken over the n test-retest paired duplicate measurements.

2.8 Histopathological validation

After the final imaging, animals were euthanised with an overdose of sodium pentobarbitol (Dolethal®). Tumours were excised, fixed in 10% neutral buffered formalin and embedded in paraffin. Aligned tumour sections (8pm) were subsequently cut from the axial plane to reflect the central B-mode ultrasound images in the XZ direction and stained with hematoxylin and eosin (H8<.E) and for the murine vascular endothelial marker CD31 antibody (DIA-310, Dianova, Hamburg, Germany) in combination with Rat Histofine Max PO (414311F, Nichirei Bioscience, Tokyo, Japan). Whole-slide images were digitized using a Zeiss Axioscan slide scanner (20x magnification, 0.46pm resolution, Zeiss, Germany) and analysed using QuPath software[33]. The percentage of necrosis (large area of tissue damage) was calculated from regions of interest drawn manually by an experienced observer (YJ) blinded to the study parameters. This was expressed as the ratio of area of necrotic region to total tumour area.

For the qualitative comparison of histological images and shear wave speed distribution, the digitized whole-slide H&E-stained images were visually aligned with the Visualsonics B-mode images (including rotation/reflection) using anatomical landmarks, including the shape of the tumour and attached skin. The observer was blinded to the processed 3D- VSWE-derived parametric shear wave speed maps, which were created by selecting the plane of shear wave speed measurement points at the centre of the tumour (i.e., the central m x n measurement points) using a 2 x 2 x 2 mm kernel for calculation.

2.9 Statistical analysis

Statistical analysis was performed using GraphPad Prism (GraphPad Software Inc., La Jolla, USA). Any significant differences in quantitative parameters were determined using Wilcoxon test for the repeatability study, and any significant differences between groups were determined using the Mann-Whitney U-Test with a 5% level of significance.

3. Results

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

3.1 Optimization and repeatability of preclinical VSWE

We optimised our 3D-VSWE system and tested its repeatability over three consecutive days in a cohort of five mice bearing subcutaneous tumour xenografts derived from the injection of MDA-MB-231 breast cancer cells. Overall, we observed that c _s and C for the whole tumour remained similar for kernel sizes of dimension 0.5, 1, 2, and 3 mm (Appendix A, Figure 28a and 28d). The GoF error was significantly less for 2 mm compared to 1 mm (Figure 28c.). The mean change in c _s over 24 hours was similar for all kernel sizes (Figure 28b). With respect to vibrational frequency, c _s and GoF both increased in magnitude with increasing frequency (700, 1000 and 1200 Hz), while C showed a reverse trend demonstrating that shear waves of higher frequency did not penetrate the tumour to the same depths as lower frequencies (Figure 28d) leading to greater noise in the displacement data. Subsequent results have been obtained using a 2 mm kernel dimension. There were no significant differences in median c _s , C or GoF values measured using one or two vibrational sources (Figure 29). Consistent with the viscoelastic properties of tissue [18], there was also a trend indicating increasing c _s values with increasing vibrational frequencies (Figures 28 and 29). There was no significant difference in measured whole tumour c _s over 24 hours. Example maps of c _s , C and GoF and repeatability data for a 2mm kernel size are shown in Figure 26. The within-subject coefficients of variation (CVWS) using one shaker were 5.2 %, 6.4 % and 14.1 % at 700, 1000 and 1200 Hz, respectively, and were greater when using two vibrational sources. CVWS was poor for the slope of c _s vs frequency (> 20%). Finally, kernel size had no effect on within-subject coefficient of variation.

See Figure 25 and corresponding brief description.

See Figure 26 and corresponding brief description.

3.2 3D-VSWE can detect treatment-induced changes in tumour tissue integrity

We set out to show proof-of-concept of the ability of our 3D-VSWE preclinical system to detect therapy-induced changes in tissue elastic modulus resulting from reduction in tissue integrity due to necrosis. To this end, we chose to measure changes in tumour c _s following the well-characterised pathophysiological response to the vascular disrupting agent ZD6126 in two xenograft models. We did not observe any change in shear wave speed dispersion, as indicated by the slope of the linear fit to tumour-median shear wave speed as a function of vibration frequency, from baseline in either treated or vehicle tumours in either model. In the vehicle cohort, there was no significant difference in measured whole tumour c _s over 24 hours, with measured tumour-median values of c _s at baseline of 4.0 ± 0.7 m/s (n=9) and c _s = 4.3 ± 0.3 m/s (n=4) for MDA-MB-231 and U-87 MG tumours, respectively (mean ± SD). Twenty-four hours after treatment with ZD6126, an overall decrease in shear wave speed was apparent across the 3D parametric c _s maps of ZD6126- treated tumours, contrasting with the stable pattern seen in vehicle tumours (Figure 26a). This was corroborated by the quantitative analysis at 700 Hz showing a reduction in c _s of both tumour models (Figure 26b and 26e), resulting in significant differences in the relative changes in c _s over 24 h between treated and vehicle cohort in both models (U-87 MG: Ac _s 24h = -24.7 ± 2.5 % vs 7.5 ± 7.1%, p=0.002 and MDA-MB-231 : Ac _s 24h = -12.3 ± 2.7 % vs 4.5 ± 4.7%, p=0.02). These results were corroborated at 1000 Hz only for the MDA- MB-231 tumours (Ac _s 24h = -23.6 ± 2.4 % vs -0.0 ± 3.0% respectively, p=0.007). Histopathological analysis of the tumours excised 24h after ZD6126 treatment showed an extensive central area of tissue damage (haemorrhagic necrosis) surrounded by a thin viable rim in all tumours - the characteristic pattern of response associated with successful treatment with a high dose of ZD6126 [6]. Quantitative analysis (Figure 26c and 26f) shows significantly higher mean necrosis area in ZD6126-treated cohort compared to vehicle cohort in both tumour models (U-87 MG: 86.0 ± 4.5 % vs 15.7 ± 5.0%, respectively, p=0.001 and MDA-MB-231 : 91.9 ± 3.9 % vs 29.6 ± 4.0%, respectively, p=0.02). Representative images of ZD6126 and vehicle group tumours are shown in Appendix A Figure 30.

3.3 3D-VSWE can provide map of tumour tissue integrity

Control U-87-MG and MDA-MB-231 tumours present necrosis foci varying in number, location, shape and size, which enabled us to evaluate the ability of our 3D-VSWE to map spatial heterogeneity in tissue integrity. We performed a qualitative comparison of the spatial distribution of c _s with the distribution of necrosis and viable regions in aligned haematoxylin and eosin (H8<.E) stained tumour slides. The maps of c _s showed good visual agreement with the spatial distribution of viable (high c _s ) and manually-segmented regions of necrosis (low c _s ) in both tumour models, further highlighting the sensitivity of our 3D- VSWE to tissue damage and demonstrating ability to resolve localised tissue damage (Figure 27).

See Figure 6 and corresponding brief description.

4. Discussion

Authors should discuss the results and how they can be interpreted from the perspective of previous studies and of the working hypotheses. The findings and their implications should be discussed in the broadest context possible. Future research directions may also be highlighted. Beyond their potential for early-disease detection, robust imaging methodologies enabling the non-invasive measure of tissue biomechanics and a comprehensive understanding of their topological variations are key to accelerate both the preclinical and clinical development of novel biomechanics-targeted therapeutic strategies, as well as providing more generic biomarkers of response to therapy, including extensive cell death, via the intricate relation between biomechanical properties and microstructure remodelling [34]. Tumour mechanical properties as measured by elastography are emerging as valuable imaging biomarkers for response assessment to chemotherapy and radiotherapy treatment, both preclinically and clinically [8-10,12,15,16,35-38] . Here, we have developed and implemented a dedicated small animal three-dimensional 3D-VSWE system which can reproducibly measure shear wave speed and its variation with frequency. Our development of 3D VSWE for preclinical imaging aimed to overcome the established limitations of 2D imaging and provide whole-tumour mapping of c _s and increase the repeatability of measurements of c _s , as has been shown clinically[39] . The within-subject coefficient of variation value for c _s measurement of our 3D-VSWE was found to be excellent for 700 and 1000 Hz (CVcs < 6.5%), comparing favourably with the repeatability of other elastography techniques [6,37,40,41]. Median tumour c _s prior to treatment (measured using 700 Hz) was ~ 4 m/s, corresponding an elastic modulus of ~ 48 kPa which falls within previously reported ranges of elastic moduli of the central slice of subcutaneous tumours measured using a commercial clinical 2D ultrasound system [35,36,42]. The relative change in c _s measured with our 3D-VSWE upon treatment with ZD6126 in MDA- MB-231 and U-87 MG tumour xenograft corroborates those measured with MR elastography (Ac _s 24h ZD6126 =-15± 2 %) in human colon carcinoma SW620 xenograft at 7T (1000Hz)[38].

Others have investigated ultrasound-based shear wave elastography for preclinical imaging. Some have used transient acoustic radiation forces impulses to generate shear waves to image tumours ex vivo [43] [44] and in vivo [35] [42] [45] [46] [47], and murine liver in vivo [G] but these studies have so far been limited to 2D imaging and analysis. Nabavizadeh et al [48] use a second transducer to generate shear waves that are less transient in nature than acoustic force impulses using 30 cycles bursts of ultrasound (duration 0.6 seconds). Using ultrafast imaging (1000) fps allows a 2D slice of tissue motion within a 2D slice.

Although our technique uses a moving ID array transducer (that collects 2D images) to build up a 3D volume, the analysis we perform is three-dimensional. Using a continuous vibrational source and repeatedly sampling the tissue motions at each transducer position, we build up a 3D map of tissue motion due to shear wave propagation. This enables us to identify the direction of the shear waves in three dimensions using our direction detector (as opposed to using directional filtering) and therefore measurements of shear wave speed are not subject to measurement bias from out of plane shear wave propagation [20].

Current MR elastography has higher spatial resolution than that achieved with our current system, however, the difference in acquisition time for one vibrational frequency (~ 1 min compared to 12 min for MRE), in addition to the established advantages of ultrasound over MRI in availability, footprint, and costs make our 3D-VSWE an attractive alternative for high-throughput longitudinal monitoring of tumour biomechanics in preclinical studies in small animals. Another advantage of our 3D-VSWE method is that it can be tailored in terms of shaker frequency and position to suit tumours of different type (cell line, stage of growth, size, type of treatment, time after treatment, etc.). For example, in both tumour types studied here, there was a trend for increased variability and smaller reduction in shear wave speed post treatment for 1000 Hz and 1200 Hz compared to 700 Hz. It is important to note that 700 Hz may not be optimal for all types of tumours, or tumour location. The biomechanical properties of the tumour will influence shear wave attenuation and therefore the ability to tune the vibrational frequency is essential. Stiffer tumour models may attenuate shear waves more rapidly and therefore we have developed our technique to calculate and display representative maps of C and GoF (Figure 25a) during data acquisition, helping the user select the position of the contactor and the frequency of vibration to obtain the best shear wave field data.

Following the use of external vibration sources in clinical reverberant field elastography [49] we have also looked at implementing two vibrational sources, instead of one, although our aims were somewhat different. Rather than explicitly trying to generate a reverberant field, our aim was to increase the amount of vibrational energy detectable through the tumour, to compensate for the rapid attenuation of shear waves at higher frequencies. Although we confirmed our hypothesis, we did not see a benefit in terms of repeatability or sensitivity of the technique to changes in shear wave speed and concluded that the use of one vibrational source only was appropriate for using with subcutaneous tumour models. This has a practical benefit, as one source is easier to incorporate in the imaging set-up.

Our method has limitations. In terms of spatial resolution, we identified that kernel dimension of 2 mm provided the best fit (smaller GoF error) for estimation of c _s as evidenced by the accurate representation of the spatial distribution of tissue shear wave speed. We have also shown that kernel sizes of 0.5 mm are feasible but further work is required to improve the signal to noise ratio in characterising the shear wave field at the smaller kernel sizes to maintain an acceptable degree of variation in shear wave speed estimates. Noise reduction methods could potentially be implemented to improve the signal to noise ratio: sampling the shear wave field for a longer period, or at an increased sampling rate, may reduce noise in the shear wave displacement-time data, which in turn can reduce measurement variation. By implementing plane-wave imaging, the above suggestions could be attempted without compromising speed of acquisition. Further improvements in spatial resolution could also be achieved by using a higher frequency transducer (our system currently uses 18.5 MHz) and wide aperture ultrasound focusing in more than one dimension, achievable with a 2D matrix array transducer [50]. The ability to acquire scans rapidly means that multiple frequencies may also be used to measure multi-frequency shear wave dispersion, which is sensitive to the underlying tissue architecture, molecular composition and is influenced by the solid/liquid state of the tissue[51]. Using 700, 1000 and 1200 Hz did not provide a reliable measure of shear wave dispersion, which may be due to the variable c _s values obtained at the higher frequencies. In addition to improvements in spatial resolution, next steps include the investigation of the range and number of frequencies that allow us to reliably characterise tumour viscous modulus.

5. Conclusions

We have developed a preclinical 3D vibrational elastography system for the rapid evaluation of tumour stiffness in vivo which can quantify and map regional heterogeneity in tissue shear wave speed throughout the whole tumour. With the speed, availability and low cost of ultrasound, our 3D-VSWE system represents an invaluable tool to enable the non-invasive assessment of cancer therapies. These include therapies targeting the ECM and mechano-transduction and to further understanding of how the biomechanical properties of tumours influence tumour progression, resistance to therapy, and metastatic potential.

Appendix A: See Figures 28-30 and corresponding brief descriptions.

Appendix B

Excluded datasets: For the repeatability study, one dataset at 1000 Hz using two shakers was excluded as it gave shear waves speeds more than 9 m/s which was ~100% greater than values obtained at 700 Hz and 1200 Hz. This was inconsistent with all other recorded datasets, in which 1000 Hz gave shear wave speeds higher than 700 Hz and lower than 1200 Hz. This dataset had the lowest conformance and greatest GoF error compared to the data obtained with two shakers at 700 Hz and 1200 Hz in the same imaging session. For the therapeutic challenge study, one mouse bearing an MDA-MB-231 tumour was excluded as endpoint histology shown none of the hallmarks of ZD6126 treatment such a central necrosis, reduced vascular density or presence of large number of cells undergoing cell death.

References

1. Northey, J. J.; Przybyla, L.; Weaver, V.M. Tissue force programs cell fate and tumor aggression. Cancer Discovery 2017, 7, 1224-1237, doi: 10.1158/2159-8290. cd-16-0733. 2. Northcott, J.M.; Dean, I.S.; Mouw, J.K.; Weaver, V.M. Feeling stress: the mechanics of cancer progression and aggression. Frontiers in Cell and Developmental Biology 2018, 6, doi: 10.3389/fcell.2018.00017.

3. Nagelkerke, A.; Bussink, J.; Rowan, A.E.; Span, P.N. The mechanical microenvironment in cancer: How physics affects tumours. In Proceedings of the Seminars in cancer biology, 2015; pp. 62-70.

4. Mohammadi, H.; Sahai, E. Mechanisms and impact of altered tumour mechanics. Nature cell biology 2018, 20, 766-774.

5. Dewhirst, M.W.; Secomb, T.W. Transport of drugs from blood vessels to tumour tissue. Nature reviews. Cancer 2017, 17, 738-750, doi: 10.1038/nrc.2017.93.

6. Li, J.; Jamin, Y.; Boult, J.K.R.; Cummings, C.; Waterton, J.C.; Ulloa, J.; Sinkus, R.; Bamber, J.C.; Robinson, S.P. Tumour biomechanical response to the vascular disrupting agent ZD6126 in vivo assessed by magnetic resonance elastography. British Journal of Cancer 2014, 110, 1727-1732.

7. Jamin, Y.; Boult, J.K.R.; Li, J.; Popov, S.; Garteiser, P.; Ulloa, J.L.; Cummings, C.; Box, G.; Eccles, S.A.; Jones, C.; et al. Exploring the biomechanical properties of brain malignancies and their pathologic determinants in vivo with magnetic resonance elastography. Cancer Res 2015, 75, 1216-1224, doi: 10.1158/0008-5472. CAN-14-1997.

8. Juge, L.; Doan, B.-T.; Seguin, J.; Albuquerque, M.; Larrat, B.; Mignet, N.; Chabot, G.G.; Scherman, D.; Paradis, V.; Vilgrain, V.; et al. Colon tumor growth and antivascular treatment in mice: complementary assessment with MR elastography and diffusion- weighted MR imaging. Radiology 2012, 264, 436-444, doi:doi: 10.1148/radiol.12111548.

9. Pepin, K.M.; Chen, J.; Glaser, K.J.; Mariappan, Y.K.; Reuland, B.; Ziesmer, S.; Carter, R.; Ansell, S.M.; Ehman, R.L.; McGee, K.P. MR elastography derived shear stiffness— a new imaging biomarker for the assessment of early tumor response to chemotherapy. Magnetic Resonance in Medicine 2014, 71, 1834-1840, doi: 10.1002/mrm.24825.

10. Feng, Y.; Clayton, E.H.; Okamoto, R.J.; Engelbach, J.; Bayly, P.V.; Garbow, J.R. A longitudinal magnetic resonance elastography study of murine brain tumors following radiation therapy. Physics in Medicine and Biology 2016, 61, 6121.

11. Garteiser, P.; Doblas, S.; Daire, J.-L.; Wagner, M.; Leitao, H.; Vilgrain, V.; Sinkus, R.; Van Beers, B.E. MR elastography of liver tumours: value of viscoelastic properties for tumour characterisation. European Radiology 2012, 22, 2169-2177, doi: 10.1007/s00330- 012-2474-6.

12. Reiss-Zimmermann, M.; Streitberger, K.-J.; Sack, I.; Braun, J.; Arlt, F.; Fritzsch, D.; Hoffmann, K.-T. High resolution imaging of viscoelastic properties of intracranial tumours by multi-frequency magnetic resonance elastography. Clinical Neuroradiology 2015, 25, 371-378, doi: 10.1007/s00062-014-0311-9.

13. Bohte, A.E.; Nelissen, J.L.; Runge, J.H.; Holub, O.; Lambert, S.A.; de Graaf, L.; Kolkman, S.; van der Meij, S.; Stoker, J.; Strijkers, G .; et al. Breast magnetic resonance elastography: a review of clinical work and future perspectives. NMR in Biomedicine 2018, 31, e3932, doi:doi : 10.1002/nbm.3932.

14. Riegler, J.; Labyed, Y.; Rosenzweig, S.; Javinal, V.; Castiglioni, A.; Dominguez, C.X.; Long, J.E.; Li, Q.; Sandoval, W.; Junttila, M.R. Tumor elastography and its association with collagen and the tumor microenvironment. Clinical Cancer Research 2018, 24, 4455-4467.

15. Evans, A.; Armstrong, S.; Whelehan, P.; Thomson, K.; Rauchhaus, P.; Purdie, C.; Jordan, L.; Jones, L.; Thompson, A.; Vinnicombe, S. Can shear-wave elastography predict response to neoadjuvant chemotherapy in women with invasive breast cancer? British journal of cancer 2013, 109, 2798-2802.

16. Evans, A.; Whelehan, P.; Thomson, K.; McLean, D.; Brauer, K.; Purdie, C.; Baker, L.; Jordan, L.; Rauchhaus, P.; Thompson, A. Invasive breast cancer: relationship between shear-wave elastographic findings and histologic prognostic factors. Radiology 2012, 263, 673-677.

17. Lampi, M.C.; Reinhart-King, C.A. Targeting extracellular matrix stiffness to attenuate disease: From molecular mechanisms to clinical trials. Science translational medicine 2018, 10.

18. Bamber, J.; Cosgrove, D.; Dietrich, C.F.; Fromageau, J.; Bojunga, J.; Calliada, F.; Cantisani, V.; Correas, J.-M.; D'onofrio, M.; Drakonaki, E. EFSUMB guidelines and recommendations on the clinical use of ultrasound elastography. Part 1 : Basic principles and technology. Ultraschall in der Medizin-European Journal of Ultrasound 2013, 34, 169- 184.

19. Fang, C.; Konstantatou, E.; Romanos, O.; Yusuf, G.T.; Quinlan, D.J.; Sidhu, P.S. Reproducibility of 2-Dimensional Shear Wave Elastography Assessment of the Liver: A Direct Comparison With Point Shear Wave Elastography in Healthy Volunteers. J Ultrasound Med 2017, 36, 1563-1569, doi: 10.7863/ultra.16.07018.

20. Zhao, H.; Song, P.; Urban, M.W.; Kinnick, R.R.; Yin, M.; Greenleaf, J.F.; Chen, S.

Bias observed in time-of- flight shear wave speed measurements using radiation force of a focused ultrasound beam. Ultrasound Med Biol 2011, 37, 1884-1892, doi: 10.1016/j.ultrasmedbio.2011.07.012.

21. Chang, S.; Kim, M.J.; Kim, J.; Lee, M.J. Variability of shear wave velocity using different frequencies in acoustic radiation force impulse (ARFI) elastography: a phantom and normal liver study. Ultraschall Med 2013, 34, 260-265, doi: 10.1055/S-0032-1313008. 22. Piscaglia, F.; Salvatore, V.; Mulazzani, L.; Cantisani, V.; Colecchia, A.; Di Donato, R.; Felicani, C.; Ferrarini, A.; Gamal, N.; Grasso, V.; et al. Differences in liver stiffness values obtained with new ultrasound elastography machines and Fibroscan: A comparative study. Dig Liver Dis 2017, 49, 802-808, doi: 10.1016/j.dld.2017.03.001.

23. Krouskop, T.A.; Dougherty, D.R.; Vinson, F.S. A pulsed Doppler ultrasonic system for making noninvasive measurements of the mechanical properties of soft tissue. J Rehabil Res Dev 1987, 24, 1-8.

24. Muthupillai, R.; Lomas, D.J.; Rossman, P.J.; Greenleaf, J.F.; Manduca, A.; Ehman,

R.L. Magnetic resonance elastography by direct visualization of propagating acoustic strain waves. Science 1995, 269, 1854-1857, doi: 10.1126/science.7569924.

25. Bayly, P.V.; Garbow, J.R. Pre-clinical MR elastography: Principles, techniques, and applications. J Magn Reson 2018, 291, 73-83, doi: 10.1016/j.jmr.2018.01.004.

26. Zeng, Q.; Honarvar, M.; Schneider, C.; Mohammad, S.K.; Lobo, J.; Pang, E.H.T.; Lau, K.T.; Hu, C.; Jago, J.; Erb, S.R.; et al. Three-Dimensional Multi-Frequency Shear Wave Absolute Vibro-Elastography (3D S-WAVE) With a Matrix Array Transducer: Implementation and Preliminary In Vivo Study of the Liver. IEEE transactions on medical imaging 2021, 40, 648-660, doi: 10.1109/tmi.2020.3034065.

27. Deffieux, T.; Gennisson, J.-L.; Bercoff, J.; Tanter, M. On the effects of reflected waves in transient shear wave elastography. IEEE transactions on ultrasonics, ferroelectrics, and frequency control 2011, 58, 2032-2035.

28. Blakey, D.C.; Westwood, F.R.; Walker, M.; Hughes, G.D.; Davis, P.D.; Ashton,

S.E.; Ryan, A. J. Antitumor activity of the novel vascular targeting agent ZD6126 in a panel of tumor models. Clinical Cancer Research 2002, 8, 1974-1983.

29. Loupas, T.; Powers, J.; Gill, R.W. An axial velocity estimator for ultrasound blood flow imaging, based on a full evaluation of the Doppler equation by means of a two- dimensional autocorrelation approach. IEEE transactions on ultrasonics, ferroelectrics, and frequency control 1995, 42, 672-688.

30. Hoyt, K.; Castaneda, B.; Parker, KJ. Two-dimensional sonoelastographic shear velocity imaging. Ultrasound in medicine & biology 2008, 34, 276-288.

31. Belevich, I.; Joensuu, M.; Kumar, D.; Vihinen, H.; Jokitalo, E. Microscopy Image Browser: A Platform for Segmentation and Analysis of Multidimensional Datasets. PLoS Biol 2016, 14, el002340, doi: 10.1371/journal.pbio.1002340.

32. Ge, X.; Quirk, J.D.; Engelbach, J. A.; Bretthorst, G.L.; Li, S.; Shoghi, K.I.; Garbow, J.R.; Ackerman, J. . Test-retest performance of a 1-hour multiparametric MR image acquisition pipeline with orthotopic triple-negative breast cancer patient-derived tumor xenografts. Tomography 2019, 5, 320-331. 33. Bankhead, P.; Loughrey, M.B.; Fernandez, J. A.; Dombrowski, Y.; McArt, D.G.; Dunne, P.D.; McQuaid, S.; Gray, R.T.; Murray, L.J.; Coleman, H.G.; et al. QuPath: Open source software for digital pathology image analysis. Sci Rep 2017, 7, 16878, doi: 10.1038/S41598-017-17204-5.

34. Li, J.; Zormpas-Petridis, K.; Boult, J.K.R.; Reeves, E.L.; Heindl, A.; Vinci, M.; Lopes, F.; Cummings, C.; Springer, C.J.; Chesler, L.; et al. Investigating the Contribution of Collagen to the Tumor Biomechanical Phenotype with Noninvasive Magnetic Resonance Elastography. Cancer Res 2019, 79, 5874-5883, doi: 10.1158/0008-5472. CAN-19-1595.

35. Chamming's, F.; Le-Frere-Belda, M.A.; Latorre-Ossa, H.; Fitoussi, V.; Redheuil, A.; Assayag, F.; Pidial, L.; Gennisson, J.L.; Tanter, M.; Cuenod, C.A.; et al. Supersonic Shear Wave Elastography of Response to Anti-cancer Therapy in a Xenograft Tumor Model. Ultrasound Med Biol 2016, 42, 924-930, doi: 10.1016/j.ultrasmedbio.2015.12.001.

36. Wang, H.; Mills, B.; Mislati, R.; Ahmed, R.; Gerber, S.A.; Linehan, D.; Doyley, M.M. Shear Wave Elastography Can Differentiate between Radiation-Responsive and Non- responsive Pancreatic Tumors: An ex Vivo Study with Murine Models. Ultrasound in medicine & biology 2020, 46, 393-404, doi: 10.1016/j.ultrasmedbio.2019.10.005.

37. Ma, Y.; Zhang, S.; Zang, L.; Li, J.; Li, J.; Kang, Y.; Ren, W. Combination of shear wave elastography and Ki-67 index as a novel predictive modality for the pathological response to neoadjuvant chemotherapy in patients with invasive breast cancer. Eur J Cancer 2016, 69, 86-101, doi: 10.1016/j.ejca.2016.09.031.

38. Li, J.; Jamin, Y.; Boult, J.K.; Cummings, C.; Waterton, J.C.; Ulloa, J.; Sinkus, R.; Bamber, J.C.; Robinson, S.P. Tumour biomechanical response to the vascular disrupting agent ZD6126 in vivo assessed by magnetic resonance elastography. Br J Cancer 2014, 110, 1727-1732, doi: 10.1038/bjc.2014.76.

39. Choi, H.; Sohn, Y.-M.; Seo, M. Comparison of 3D and 2D shear-wave elastography for differentiating benign and malignant breast masses: focus on the diagnostic performance. Clinical radiology 2017, 72, 878-886.

40. Shin, Y.; Huh, J.; Ham, S.J.; Cho, Y.C.; Choi, Y.; Woo, D.C.; Lee, J.; Kim, K.W. Test-retest repeatability of ultrasonographic shear wave elastography in a rat liver fibrosis model: toward a quantitative biomarker for preclinical trials. Ultrasonography 2021, 40, 126-135, doi: 10.14366/usg.19088.

41. Weis, J. A.; Flint, K.M.; Sanchez, V.; Yankeelov, T.E.; Miga, M.I. Assessing the accuracy and reproducibility of modality independent elastography in a murine model of breast cancer. J Med Imaging (Bellingham) 2015, 2, 036001, doi: 10.1117/1. MI.2.3.036001.

42. Elyas, E.; Papaevangelou, E.; Alles, E.J.; Erler, J.T.; Cox, T.R.; Robinson, S.P.; Bamber, J.C. Correlation of Ultrasound Shear Wave Elastography with Pathological Analysis in a Xenografic Tumour Model. Sci Rep 2017, 7, 165, doi: 10.1038/s41598-017- 00144-5.

43. Wang, H.; Mislati, R.; Ahmed, R.; Vincent, P.; Nwabunwanne, S.F.; Gunn, J.R.; Pogue, B.W.; Doyley, M.M. Elastography Can Map the Local Inverse Relationship between Shear Modulus and Drug Delivery within the Pancreatic Ductal Adenocarcinoma Microenvironment. Clin Cancer Res 2019, 25, 2136-2143, doi: 10.1158/1078-0432. CCR- 18-2684.

44. Wang, H.X.; Mills, B.; Mislati, R.; Ahmed, R.; Gerber, S.A.; Linehan, D.; Doyley, M.M. Shear Wave Elastography Can Differentiate between Radiation-Responsive and Non- Responsive Pancreatic Tumors: An Ex Vivo Study with Murine Models. Ultrasound in Medicine and Biology 2020, 46, 393-404, doi: 10.1016/j.ultrasmedbio.2019.10.005.

45. Baek, J.; Ahmed, R.; Ye, J.; Gerber, S.A.; Parker, K .; Doyley, M.M. H-scan, Shear

Wave and Bioluminescent Assessment of the Progression of Pancreatic Cancer Metastases in the Liver. Ultrasound Med Biol 2020, 46, 3369-3378, doi: 10.1016/j.ultrasmedbio.2020.08.006.

46. Riegler, J.; Labyed, Y.; Rosenzweig, S.; Javinal, V.; Castiglioni, A.; Dominguez, C.X.; Long, J.E.; Li, Q.; Sandoval, W.; Junttila, M.R.; et al. Tumor Elastography and Its Association with Collagen and the Tumor Microenvironment. Clin Cancer Res 2018, 24, 4455-4467, doi: 10.1158/1078-0432. CCR-17-3262.

47. Mieulet, V.; Garnier, C.; Kieffer, Y.; Guilbert, T.; Nemati, F.; Marangoni, E.; Renault, G.; Chamming's, F.; Vincent-Salomon, A.; Mechta-Grigoriou, F. Stiffness increases with myofibroblast content and collagen density in mesenchymal high grade serous ovarian cancer. Sci Rep 2021, 11, 4219, doi: 10.1038/s41598-021-83685-0.

48. Nabavizadeh, A.; Payen, T.; luga, A.C.; Sagalovskiy, I.R.; Desrouilleres, D.; Saharkhiz, N.; Palermo, C.F.; Sastra, S.A.; Oberstein, P.E.; Rosario, V.; et al. Noninvasive Young's modulus visualization of fibrosis progression and delineation of pancreatic ductal adenocarcinoma (PDAC) tumors using Harmonic Motion Elastography (HME) in vivo. Theranostics 2020, 10, 4614-4626, doi: 10.7150/thno.37965.

49. Ormachea, J.; Castaneda, B.; Parker, KJ. Shear Wave Speed Estimation Using Reverberant Shear Wave Fields: Implementation and Feasibility Studies. Ultrasound Med Biol 2018, 44, 963-977, doi: 10.1016/j.ultrasmedbio.2018.01.011.

50. Jensen, J. A. Medical ultrasound imaging. Prog Biophys Mol Biol 2007, 93, 153-165, doi: 10.1016/j.pbiomolbio.2006.07.025.

51. Garteiser, P.; Page, G.; d'Assignies, G.; Leitao, H.S.; Vilgrain, V.; Sinkus, R.; Van Beers, B.E. Necro-inflammatory activity grading in chronic viral hepatitis with three- dimensional multifrequency MR elastography. Scientific reports 2021, 11, 1-11. 52. Workman, P.; Aboagye, E.O.; Balkwill, F.; Balmain, A.; Bruder, G.; Chaplin, DJ.; Double, J. A.; Everitt, J.; Farningham, D.A.; Glennie, M .; et al. Guidelines for the welfare and use of animals in cancer research. Br J Cancer 2010, 102, 1555-1577, doi: 10.1038/sj.bjc.6605642. 53. Kilkenny, C.; Browne, W.J.; Cuthill, I.C.; Emerson, M.; Altman, D.G. Improving bioscience research reporting: The ARRIVE guidelines for reporting animal research. J Pharmacol Pharmacother 2010, 1, 94-99, doi: 10.4103/0976-500X.72351.