Stretchable ultrasonic arrays for the three-dimensional mapping of the modulus of deep tissue

Serial assessment of the biomechanical properties of tissues can be used to aid the early detection and management of pathophysiological conditions, to track the evolution of lesions and to evaluate the progress of rehabilitation. However, current methods are invasive, can be used only for short-term measurements, or have insufficient penetration depth or spatial resolution. Here we describe a stretchable ultrasonic array for performing serial non-invasive elastographic measurements of tissues up to 4 cm beneath the skin at a spatial resolution of 0.5 mm. The array conforms to human skin and acoustically couples with it, allowing for accurate elastographic imaging, which we validated via magnetic resonance elastography. We used the device to map three-dimensional distributions of the Young’s modulus of tissues ex vivo, to detect microstructural damage in the muscles of volunteers before the onset of soreness and to monitor the dynamic recovery process of muscle injuries during physiotherapies. The technology may facilitate the diagnosis and treatment of diseases affecting tissue biomechanics. A stretchable ultrasonic array conforming to the skin allows for the three-dimensional imaging of tissue modulus at depths of up to 4 cm.

Serial assessment of the biomechanical properties of tissues can be used to aid the early detection and management of pathophysiological conditions, to track the evolution of lesions and to evaluate the progress of rehabilitation.However, current methods are invasive, can be used only for short-term measurements, or have insufficient penetration depth or spatial resolution.
Here we describe a stretchable ultrasonic array for performing serial non-invasive elastographic measurements of tissues up to 4 cm beneath the skin at a spatial resolution of 0.5 mm.The array conforms to human skin and acoustically couples with it, allowing for accurate elastographic imaging, which we validated via magnetic resonance elastography.We used the device to map three-dimensional distributions of the Young's modulus of tissues ex vivo, to detect microstructural damage in the muscles of volunteers before the onset of soreness and to monitor the dynamic recovery process of muscle injuries during physiotherapies.The technology may facilitate the diagnosis and treatment of diseases affecting tissue biomechanics.
The mechanical properties of human tissues are vital for the structure and function of human physiological systems 1 .Frequent mechanical characterizations of various organs allow timely evaluation of tissue growth, metabolic state, immunologic function and hormone regulation [1][2][3] .Most importantly, the mechanical properties of diseased tissues can often reflect pathophysiological conditions.The monitoring of such properties can provide key information about disease progression, and guide intervention in a timely manner 4,5 .For instance, the Article https://doi.org/10.1038/s41551-023-01038-walgorithm is used to compare the radiofrequency data before and after compression, and calculate the displacements of the scattering sources with high sonographic sensitivity and accuracy (Supplementary Note 7 and Supplementary Fig. 12a,b) 29,30 .A least-squares strain estimator is used to transform displacements to strain while minimizing possible fluctuations (Supplementary Note 8 and Supplementary Fig. 12) 31 .A graphical user interface that consists of real-time display windows and control panels is designed for locating the device and quantifying the compression to the subject (Supplementary Note 9 and Supplementary Fig. 13).An inverse elasticity problem is solved to quantify the modulus distribution inside the tissue (Supplementary Note 10).
We choose a centre frequency of 3 MHz to balance the requirement of high spatial resolution 32 and frequency-dependent linear attenuation of ultrasound waves in tissues 33 .The characterized mean resonant and anti-resonant frequencies show small standard deviations, indicating the consistency across the entire array (Fig. 1b).Given the corresponding ultrasound wavelength of ~500 μm in soft tissues, we choose a pitch of 800 μm that is suitable for generating wave convergence 34 , producing high-quality images and minimizing crosstalk 35 .The pitch is altered upon conformation to a curvature.However, in this study, we continue using the initial planar pitch of 800 μm in the beamforming of a curved surface, due to the negligible changes of no more than 0.003% in the pitch within the adoptable curvature range of the array (Supplementary Note 11 and Supplementary Fig. 14).A scalable method is used to align the backing material and the transducer elements, which enhances the fabrication throughput and performance consistency, and avoids potential phase aberrations (Supplementary Fig. 15 and Methods) 35,36 .To individually address the 256 elements, six layers of activation electrodes and one layer of common ground in a serpentine shape are transfer-printed and routed to the same plane by vertical interconnect accesses (Supplementary Figs.16-18).A low-temperature bonding technique is used to preserve the high electromechanical coupling coefficient of the 256 elements (average 0.64, Fig. 1b inset), with performance comparable to elements in commercial ultrasound probes (0.58-0.69) (ref.37).The 0.022 of average dielectric loss of the array indicates the energy consumption caused by transducer heating is minimal (Supplementary Fig. 19).The time-and frequency-domain characterizations of pulse-echo response show ~50% bandwidth (Supplementary Fig. 19 and Methods).
Because of the well-designed pitch, the suppression of shear vibrations by the 1-3 composite, and the damping effect of silicone elastomer between the elements, crosstalk between the elements is below the standard -30 dB (Fig. 1c).A matching circuit effectively reduces the insertion loss to 16.98 dB at the resonant frequency (Supplementary Fig. 20), comparable to a commercial probe (17 dB) (ref.38).The low insertion loss leads to excellent sonographic sensitivity (Fig. 1d), which is pivotal for modulus imaging 39 .The high-performance 1-3 composite material, together with the new fabrication protocol and effective matching circuit, results in an average sonographic signal-to-noise ratio of 39 dB, comparable to a commercial probe (41 dB) (ref.39).The standard deviation of the sonographic signal-to-noise ratio measured underwater for 2 weeks is small (0.87 dB), because of the excellent waterproof property of the encapsulating silicone elastomer (Supplementary Fig. 21).The device can be reversibly deformed in various modes (Fig. 1e-g).The maximum biaxial stretchability of the device without affecting its electromechanical properties is ~40%, indicating its reliability for skin integration (Supplementary Figs.22 and 23) 40 .

Strategies for elastography
To determine the transmission mode for elastography, we simulated two-dimensional (2D) strain distributions in a bilayer phantom under three different transmission modes: single plane wave, mono-focus and coherent plane-wave compounding (Supplementary Fig. 24) 26 , and then analysed the SNR e and CNR e , the two most critical metrics for modulus mapping (Supplementary Notes 12 and 13) 41 .Among the stiffnesses of tumours are known to be different from those of healthy tissues 6 .Additionally, in some tumours, changes in stiffness can occur as they grow in certain developmental stages 7 , and these changes can happen quickly (Supplementary Fig. 1a and Supplementary Note 1) 5,[8][9][10][11][12] .Frequent inspections of the stiffness of these tumours are required for growth-stage assessment and therapy guidance 13 .Mechanical characterization is also critical in the diagnosis and rehabilitation of many musculoskeletal diseases and injuries.The monitoring of muscular moduli enables more proactive screening of the area at risk (Supplementary Fig. 1b) [14][15][16][17][18][19][20][21][22] .Surveillance of tissue moduli has also been demonstrated to help with early detection and tracking of cardiovascular diseases 23,24 .An ideal technology should provide non-invasive and three-dimensional (3D) mapping of deep tissues with accurate location, morphology and mechanical information 25 .However, existing methods are unable to address this critical need (Supplementary Figs.2-5 and Supplementary Note 2).
In this article, to fill this technological gap, we report a stretchable ultrasonic array with advances in device engineering and imaging algorithms (Supplementary Note 3).With new microfabrication protocols, we can achieve excellent transducer electromechanical coupling.The coherent compounding imaging strategy enables accurate displacement calculations, and therefore enhances elastographic signal-to-noise ratio (SNR e ) and contrast-to-noise ratio (CNR e ) in the entire sonographic window 26 .By solving an inverse elasticity problem, we can derive quantitative modulus distribution, which is a leap forward compared with the qualitative strain distribution obtained with conventional quasi-static elastography 27 (Supplementary Note 4).We show the reliability of this technology by testing on various artificial phantom models and ex vivo biological specimens, with quantitative validation by magnetic resonance elastography (MRE).In vivo studies on delayed-onset muscle soreness shows that this technology can track the recovery progress of muscular injury in a non-invasive serial manner, providing therapeutic guidance.These results suggest a convenient and effective approach to monitor tissue mechanical properties, facilitating the diagnosis and treatment of many diseases and symptoms.

Working principle, design and fabrication of the device
Measurements by traditional rigid ultrasound probes may suffer from poor acoustic coupling, because their rigid surfaces cannot accommodate the curvilinear shape of the human body (Supplementary Note 5).To solve this problem, we designed a stretchable ultrasonic array for modulus sensing.Figure 1a shows the schematic working mechanism of a 16-by-16 array (for the detailed process, see Supplementary Figs. 6 and 7 and Supplementary Movie 1).The array can conform to and acoustically couple with the human skin via a type of silicone-based couplant for high-quality imaging (Supplementary Fig. 5 and Supplementary Note 5).After the array is activated, ultrasound waves are sent into the tissue underneath the device.Scattering sources (for example, tissue interfaces) can reflect the ultrasound waves, which carry location information of these scattering sources.The reflected waves are then received by the transducer elements in the array as radiofrequency data.The collected radiofrequency data by each element are then enhanced by receive beamforming (Supplementary Fig. 8).After that, ~0.5-1% strain is applied to the tissue by quasi-static uniaxial compression.The low strain ensures that the tissues exhibit linear stress-strain behaviour-that is, their Young's modulus does not vary with load 28 .After compression, the maximum change in the array's pitch is only 0.02% (Supplementary Fig. 9).Therefore, the associated phase aberration is negligible for both transmit and receive beamforming, especially for deep-tissue inspections (Supplementary Note 6 and Supplementary Figs. 10 and 11).As a result, we can use the same time-delay profile for both post-and pre-compression receive beamforming (Supplementary Note 6).A normalized cross-correlation Article https://doi.org/10.1038/s41551-023-01038-wthree modes, the single plane wave mode produces the lowest SNR e and CNR e , because of the low transmission energy from non-focal scanning (Fig. 2a).The mono-focal mode provides better SNR e and CNR e , but only in regions near the focal point.The coherent plane-wave compounding mode consists of a series of single plane waves at different transmission angles, and these steered frames are coherently integrated to reconstruct a multi-angle compounding image, with multiplexed signal intensities across the entire region (Supplementary Note 13).The resulting SNR e and CNR e levels are substantially enhanced in all regions 26 .The performance of the compounding strategy highly depends on the step size and the number of the steering angles.In this work, we used 19 steering angles and a step size of 1°, which were found in experiments to give the best elastographic image quality (Fig. 2b,c and Supplementary Figs. 25 and 26).
A high SNR e is critical for high-quality elastographic imaging.The SNR e is influenced by the magnitude of applied strain and the normalized cross-correlation coefficient (NCC), which reflects the similarity of the radiofrequency signals before and after compression.The smallest Article https://doi.org/10.1038/s41551-023-01038-wdetectable strain with an SNR e of 6 dB is 0.0125%, which indicates the high elastographic sensitivity of the stretchable array (Fig. 2d).At low strain levels, the displacement signals are small, so the relatively large intrinsic electrical noise in the radiofrequency signals leads to large deviations in the elastographic images and thus a low SNR e .As the applied strain increases, the displacement signals increase proportionally, while the intrinsic electrical noise remains relatively constant, leading to an increase in the SNR e .The SNR e becomes the highest at ~1% strain.Beyond 1% strain, the radiofrequency signals distort in shape and amplitude 42 , resulting in a sharp fall in the NCC and therefore SNR e (Fig. 2d).The dynamic range (that is, the strain range with high SNR e ) is determined by a −3 dB cut-off from the maximum SNR e , which is ~0.36-2.34% in this work, corresponding to a NCC >0.8.
Both spatial resolution and contrast resolution are critical for elastographic imaging 43 .To characterize the spatial resolution, we imaged the modulus distribution and extracted the lateral and axial transition edges of an inclusion phantom (Supplementary Fig. 27 and  Methods).The first derivatives of the modulus distribution yield the point spread functions (PSFs) (Supplementary Fig. 28), whose full widths at half maximum (FWHMs) are used to estimate the elastographic spatial resolutions 25 , which are determined to be 0.56 mm in the lateral and 0.50 mm in the axial directions (Fig. 2e).Because of the synthetic focusing effect of the coherent plane-wave compounding mode, the spatial resolution remains high and consistent at different locations across the entire imaged region 26 .
Similar to the spatial resolution, the contrast resolution is defined as the modulus contrast with a corresponding CNR e of 6 dB.To quantify the contrast resolution, we perform tests on bilayer phantoms, composed of two homogeneous gelatin phantoms with different elastic moduli.The modulus of each layer is within ~10-100 kPa, covering the modulus range of all typical healthy and diseased tissues (Supplementary Fig. 29) 44 , which creates an interfacial modulus contrast ranging from 1.79 dB to 15.94 dB.The CNR e , determined by the ratio of the strain contrast to the standard deviation of strain distribution, decreases as the modulus contrast becomes lower (Fig. 2f and Supplementary Note 12).Logarithmic curve fitting, with a coefficient of determination >0.98, demonstrates the reliability of the experimental results (Supplementary Note 12) 43 .The contrast resolution of the device is determined to be ~1.79 dB (that is, a modulus ratio of 1.22).

Characterizations on phantom models
Following the clinical practice for tumour screening and diagnosis, we used four types of phantom to simulate different pathological tissue environments (Supplementary Fig. 27) 43,45 : a one-dimensional phantom consisting of two layers with different Young's moduli, which mimics muscles with an area of disease or injury 15,17 ; a 2D phantom consisting of a cylindrical cavity filled with fluid, which simulates cysts that frequently appear in central organs 46 ; a 2D phantom with a solid cylindrical inclusion; and a commercial 3D breast phantom with a spherical mass.The last two phantoms mimic tumours and nodes with various morphologies in the breast 47 .Young's moduli of all components have been characterized by either standard apparatus (for customized phantoms, Supplementary Fig. 29) or a clinical ultrasound machine (for the commercial phantom, see Supplementary Fig. 30).Many masses (except for the cyst phantom) have materials and constituents, and thus acoustic impedances, similar to the surrounding tissues.Thus, they exhibit a homogeneous echogenicity and a minimal sonographic contrast that can be hardly distinguished by B-mode imaging (Fig. 3a-d, first column) 48,49 .
The tests on phantom models focus on the axial displacement fields (Fig. 3a-d, second column).Compared with the lateral and elevational displacements, the axial displacements can reflect the movement of each scattering source more accurately, because it is parallel to the direction of ultrasonic wave transmission 50 .The strain is higher when the modulus of the component is lower (the third column in Fig. 3a-d and Supplementary Note 8).A common practice for computing strain is to take the spatial derivative of the displacement field 51 , which amplifies small unphysical fluctuations in the displacement field (Supplementary Fig. 12).To remedy this problem, we applied a least-squares strain estimator with a piecewise linear curve fitting method, which allows us to calculate the 2D strain distributions while smoothing unphysical fluctuations (Supplementary Fig. 12) 52 .The resulting strain distributions clearly reveal the inclusion.An exception is the cyst phantom, whose strain distribution map has a blue (B)-green (G)-red (R) region (Supplementary Figs. 31 and 32), a signature of cysts that has been used clinically to distinguish between cystic and solid lesions 46,53 .
We reconstructed 3D strain images of the phantoms by integrating 16 cross-sectional images obtained by the array.The 3D imaging results match those of a commercial ultrasound probe with a similar centre frequency (2.8 MHz) (Supplementary Figs.33-36 and Supplementary Note 14).The measurements are highly reproducible, reflecting the reliability of the stretchable array (Supplementary Fig. 37).Strain depends on applied loads, and thus strain mapping can be subjective and operator dependent.Additionally, strain maps cannot faithfully reveal quantitative modulus information if the load is non-uniform 54 .To avoid these issues, we quantify the spatial distribution of the shear modulus by solving an inverse elasticity problem 55 .Specifically, we formulate the inverse elasticity problem as a Article https://doi.org/10.1038/s41551-023-01038-wconstrained optimization problem.The objective is to seek a shear modulus distribution that produces a predicted displacement field that both satisfies the equilibrium equation of the 2D linear elasticity model 56 and matches the measured displacement field.We solve the optimization problem by a gradient-based minimization approach and compute the gradient efficiently using the adjoint method (Supplementary Fig. 38) 55 .The Young's modulus is three times the shear modulus 57,58 .For a measured displacement field, the derived modulus distribution can be accurately determined.Note that the non-uniform external compressions do not influence the reconstructed modulus distribution because the inverse approach does not assume a uniform stress field to a given cross-section (Supplementary Note 10 and Extended Data Fig. 1).
The modulus distribution maps visualize the morphology of the internal structures, which accurately match the design (the fourth column of Fig. 3a-d).Shadowing artefacts, usually caused by inclusions (Supplementary Fig. 39), do not appear in the mapping results here, because of the high transmission energy of the coherent compounding method and the excellent sonographic sensitivity of the stretchable ultrasonic array.The mean modulus contrasts between the stiff and soft components of the bilayer, inclusion and breast phantoms are 1.94, 1.50 and 2.21, with a difference of 4.59%, 6.81% and 12.78%, respectively, from those obtained by a standard apparatus.Note that, for the cyst phantom, the acquired values indicate only that the stiffness of the cyst is much lower than that of the surrounding matrix, and do not represent the exact modulus ratios (Supplementary Note 15) 59 .There is a slight underestimation of modulus contrast in all cases because the total variation regularization used to solve the inverse elasticity problem tends to sacrifice the contrast to generate less noisy images (Supplementary Fig. 40) 60 .The overall >87% accuracy of the results here is well above the average accuracy of 80% in the literature 27 , because of the outstanding transducer performance and the advanced coherent compounding approach.

Validation and serial monitoring on biological tissues
We validated the 3D imaging performance of the stretchable array against MRE on 4-cm-thick, multi-layered porcine abdominal tissues 61 .Each 1 × 16 linear array in the ultrasonic patch can map a 2D cross-sectional displacement field and the corresponding modulus distribution (Supplementary Fig. 41).A 3D elastographic image is reconstructed by integrating 16 slices of modulus maps with a 0.8 mm pitch (Fig. 4a), with an average discrepancy between the measured and predicted displacement fields of 3% (Fig. 4b and Methods).The 3D image illustrates the heterogeneous nature of the soft tissue.The same porcine abdominal tissue is then measured by MRE (Fig. 4c and Methods).The volumetric mapping result from the stretchable array is highly comparable to that from MRE in both morphology and modulus distribution.The average measured modulus contrasts of transversus abdominis muscle, obliquus internus abdominis muscle, fascia and intermuscular fat are approximately 1.92:1.67:1.30:1by the stretchable array, and 2.26:1.83:1.46:1by MRE; both are comparable to those in the literature 62,63 .
A method for serial surveillance needs to be stable for long-term reproducibility, with a high sensitivity for picking up dynamic changes in the target tissue.A longitudinal study is carried out by testing a commercial phantom with the stretchable ultrasonic array and a commercial probe over 8 weeks.We compared the strain contrast between the mass and matrix measured by the two methods and conducted Bland-Altman analysis.All data points are within 95% confidence interval, showing excellent agreement between the stretchable array and the commercial probe.Additionally, a small bias error of 0.02 shows the high accuracy of the stretchable array and a small precision error of 0.09 indicates the stability of the device for repeated measurements over the long term (Fig. 4d and Supplementary Fig. 42).To test the elastographic sensitivity of the stretchable array, we measure a piece of bovine gluteobiceps muscle 64 under controlled unilateral heating as its modulus gradually increases (Supplementary Fig. 43).The measurement lasts for 165 min.The results clearly show that, before heating, the tissue has a low modulus across the entire depth (Fig. 4e).When heated from the bottom, the region of the tissue close to the heat source starts to stiffen, because the actin in the myofibrils becomes firm and short, expelling liquid and making the structure dense 65 .As the heating continues, the high-modulus region gradually grows, while the boundary between the high-and low-modulus regions remains clearly defined.These recordings demonstrate the capability of the stretchable array for serial monitoring of deep tissue mechanics.

Multi-site mapping of the human body
In vivo measurements can further illustrate the clinical value of the stretchable ultrasonic array.Multiple sites on the human body where muscle injuries usually happen are selected in this study (Supplementary Note 16 and Supplementary Fig. 44). Figure 5a-d shows the results of strain mapping within 4 cm depth from the lateral side of the shoulder joint (Fig. 5a), the anterior forearm (Fig. 5b), the anterior thigh (Fig. 5c) and the posterior calf (Fig. 5d) with the anatomies labelled, juxtaposed with the corresponding B-mode images acquired by a commercial probe.The stretchable array can effectively resolve the mechanical properties of various tissue components (Supplementary Figs.45 and 46).
Testing on ten more volunteers confirmed the feasibility and reliability of the stretchable ultrasonic patch (Supplementary Note 16 and Supplementary Fig. 47).To verify the reliability of this device, we performed intra-session and inter-day testing on six participants 66 (Extended Data Fig. 2).To quantitatively evaluate those results, we calculated the mean values of the biceps brachii region in the strain images and used a two-way mixed-effect model 67 to derive the intraclass correlation coefficient 67 and standard error of measurement 68 (Supplementary Table 1), which indicate excellent relative and absolute reliabilities 66 , respectively, of using the stretchable ultrasonic device for mapping deep-tissue modulus.

Serial surveillance of delayed-onset muscle soreness
Over-exercise introduces injuries to the musculoskeletal system, associated with damages in the sarcolemma and others 69 .These disruptions lead to inflammation, stiffness increase and function impairment of the tissues 69 .Very often, the sensation of soreness does not occur until a few days later 70 .The delayed onset of body responses precludes timely treatments, and the injury often gets neglected and worsens.There are several studies that have reported modulus variation of the muscle after the eccentric exercise 19,71,72 .However, the initial test was performed 1 h after the participants had exercised, while by that time, muscle damage might have already occurred, making it impossible to pinpoint the exact time of the initial muscle damage.Additionally, the modulus changes were only sporadically tested over a period 19 , which was easy to miss the turning point in the change of the muscular modulus and failed to offer accurate trends of muscle recovery.Additionally, serial evaluation of the tissue can guide the rehabilitation strategy 71,72 .MRE is often used to evaluate the tissue stiffness to diagnose any injuries 71,72 .However, MRE is not viable for long-term testing, due to its large size, limited availability and high cost.
The stretchable ultrasonic array addresses these needs.A healthy volunteer was selected to perform eccentric elbow joint exercise to develop delayed-onset muscle soreness 73 (Fig. 5e, Supplementary Fig. 48 and Supplementary Note 17).Before exercise, well-defined anatomic components of the upper arm can be visualized in 3D (Supplementary Figs.49 and 50).After exercise, the muscles are either allowed to recover naturally, or treated with massotherapy or hyperthermia.In all experiments, the normalized modulus contrast of the biceps brachii is surveyed every day for 5 days to track the dynamic recovery process.Meanwhile, the intensity of the soreness is evaluated on the Article https://doi.org/10.1038/s41551-023-01038-wbasis of the pain visual analogue scale 19 (Supplementary Fig. 48).Wearing the device 1 h per day for 5 days does not induce any skin irritation (Supplementary Fig. 51).
Figure 5f shows the monitoring results, where normalized modulus contrast (top) and soreness intensity score (bottom) of each measurement are collected.The modulus of biceps brachii increases within 20 min after the exercise due to the muscle contracture induced by sarcolemmal disruption (Supplementary Note 17) 71,72 .However, the sensation of soreness does not arise until 1 day after exercise (Supplementary Note 17) 74 .As more time goes by, sarcolemmal disruption continues, and the muscle modulus keeps increasing.Meanwhile, the circulation system delivers supplies to repair the sarcolemma, and eventually tissue modulus drops (Supplementary Note 17) 75 .During natural recovery, the muscle modulus continues to increase for 2 days before dropping.When physiotherapies are applied, the increase sustains for only 1 day, with lower maximum moduli.The intensity of soreness has a latency but evolves in a trend similar to the modulus.The experiments have been performed on more participants, and the results are Di erence (stretchable device comparator)  Article https://doi.org/10.1038/s41551-023-01038-wsimilar (Extended Data Fig. 3).Blood cannot circulate well through the injured myofibrils, leading to involuntary tremors of the muscle during testing 76 , which induces variations in the measured modulus results.The coefficient of variation, defined as the ratio of the standard deviation to the mean, is smaller for physiotherapies than that of the natural recovery (Fig. 5g).These findings show the quicker healing

Article
https://doi.org/10.1038/s41551-023-01038-wpace of physiotherapies compared with natural recovery (Supplementary Note 18).The stretchable ultrasonic array demonstrates the feasibility of effectively monitoring the tissue stiffness under different muscle conditions.We expand the sample size and use maximal voluntary contraction, a standard clinical approach to evaluate muscle activity 77 , as a comparator to quantify the degree of the muscle damage after eccentric exercise (Supplementary Note 16).There is a notable drop in maximal voluntary contraction torque after doing the exercise followed by a gradual increase, which indicates the process of muscle damage followed by recovery (Extended Data Fig. 4) 78 .This trend corresponds to that of the modulus contrast variation, showing the reliability of the stretchable ultrasonic array to detect muscle activity (Extended Data Fig. 4).P value is the probability of one result equalling another under the assumption that the null hypothesis is correct.P values of both methods are all much less than 0.001, which further confirm this conclusion.

Discussion
The stretchable ultrasonic array can perform serial, non-invasive 3D mapping of the mechanical properties of deep tissues, which has not been realized by any existing devices.The array has high SNR e , CNR e , spatial resolution and contrast resolution, which can be attributed to the strategic array design, new microfabrication techniques and the effective ultrasound transmitting mode.Solving the inverse elasticity problem provides accurate and reliable quantitative information of modulus distribution in a 3D tissue.The stretchable ultrasonic array can detect muscle injuries before the subject's own feeling, and therefore enable timely intervention to prevent cumulative trauma disorders.Moving beyond muscle injuries, this technology has the potential to monitor the size and modulus of tumour in real time and inform therapeutic decisions, providing a new profiling method for both fundamental oncology research and clinical practice.Collectively, these findings demonstrate that the stretchable ultrasonic array is complementary to existing clinical monitoring modalities and can be used as a unique enabling technology for quantitatively sensing and treating many deep-tissue conditions.
Future developments may leverage advanced lithography, dicing and pick-and-place techniques to further optimize the array design and fabrication, which can reduce the pitch and extend the aperture to achieve a higher spatial resolution and a broader sonographic window.Additionally, the stretchable ultrasonic array is currently wired for data and power transmission.Those back-end tasks, undertaken by a desktop-based interfacing system, such as electronic control, pulser and receiver, and data processing, can be achieved by a flexible printed circuit board.With the advent of lower-power integrated circuits and flexible lithium-polymer battery technologies, we envision the entire hardware to be fully portable while maintaining its high performance.

Fabrication of the stretchable ultrasonic array
The fabrication of the stretchable ultrasonic array begins with preparing the multi-layered activation electrodes.To individually control each element of the array, 256 activation electrodes and 1 common ground are needed.It would be very challenging to place so many electrodes in a single layer while keeping the device layout compact and stretchable.We choose to integrate the electrodes in a multi-layered manner 79 .We use AutoCAD software to design the six layers of the activation electrodes to individually address each element (Supplementary Fig. 16), which minimizes the footprint of the entire device to a size much smaller than the stretchable ultrasonic device in our previous work 80 .The pitch of the array is shrunk from the previous 2 mm by 2 mm to the present 0.8 mm by 0.8 mm, and the kerf also has been decreased from the previous 0.8 mm to the present 0.2 mm, which enhances the active component density from the previous 39.06% to the present 58.05%.To route all elements to the same plane, vertical interconnect accesses based on laser ablation are adopted 79 .
For each single element, we use 1-3 composite as the active material given its superior electromechanical coupling properties, appropriate acoustic impedance 81 , and suppressed crosstalk between adjacent elements 35 .We choose silver-epoxy composites as the backing material.To make the backing layer more condensed and reduce the overall device thickness (~0.8 mm), uncured silver-epoxy mixture is centrifuged for 10 min at 3,000 r.p.m, which separates the extra liquid hardener, removes air bubbles in the silver-epoxy and improves the acoustic damping effect and the axial resolution.
The array with such a small footprint and large scale requires optimized fabrication techniques.We develop a new method of automatic alignment of the transducer elements to the bonding electrodes.We first bond a large piece of the backing layer to the 1-3 composite and then dice the bonded bilayer to the designed array configuration.A strong adhesion tape is used to fix the bilayer without delamination during dicing by the dicing saw (DAD3220, DISCO).The 2D array with designed element size and pitch is automatically formed.To prevent the array from tilting or moving when the array is bonded to the electrodes, a silicone elastomer (Ecoflex-0030, Smooth-On) is filled in the kerf and connects individual elements together.
Then conductive epoxy (Von Roll 3022 E-Solder) is used to bond the copper-based stretchable interconnects with the 1-3 composite.The device is put at room temperature for 8 h and at 40 °C for 2 h for tight bonding.High-temperature bonding in previous studies 35,36 generates serious damages to piezoelectric materials.First, dipoles in the piezoelectric materials are lost under high temperatures, which seriously damages its ability to transmit and receive ultrasonic waves, causing low signal-to-noise ratios of measurements.Although the dipoles can be re-aligned after applying an external electric field, the polarization is very time consuming.Additionally, an excessive electrical field will break down the piezoelectric element 36 .Second, high-temperature bonding will cause thermal damage to the epoxy in 1-3 composites.It softens and irreversibly deforms the epoxy so that the alignment of piezoelectric pillars in 1-3 composites is reduced.As a result, their electromechanical coupling performance decays.The low-temperature bonding method in this study avoids thermal damage of the epoxy and any possible depolarization of the piezoelectric materials in the 1-3 composite, which maximizes the piezoelectric performance of the elements.Before bonding the other electrode to the transducer array, the interface between the transducer array and the dicing tape is exposed to UV light (254 nm wavelength, PSD series Digital UV Ozone System, Novascan) for 10 min for surface de-adhesion.Then the array can be delaminated from the tape and bond with the other electrode using the same method.
As shown in Supplementary Fig. 27, gelatin is added into boiling de-ionized water.1-Propanol (Sigma-Aldrich) is dropped into the solution (9.2 g per 100 g de-ionized water) to facilitate the crosslinking of gelatin.Silicon dioxide particles (3% mass fraction of all ingredients) are added into the solution after the solution becomes clear.The particles serve as the markers for displacement tracking.After thorough mixing for 10 min, the suspension is poured into specific molds for cooling.When we make the bilayer phantoms for measurements of CNR e and those in Fig. 3a, the bottom layer is cured into a cuboid geometry first and then the suspension of the second layer is poured and cured Article https://doi.org/10.1038/s41551-023-01038-w on the cuboid.Before pouring, the solution needs to be cooled to ~40 °C to prevent melting and osmosis of the bottom phantom by the hot solution 82 .
The process of making the cyst and inclusion phantoms is different from making the bilayer phantom.To build a cylindrical cavity in the matrix, a tube with a 9.2 mm outer diameter is inserted into the uncured matrix solution.The tube is then removed after the matrix is cured.Solution or water is filled into the cavity to form the inclusion or cyst.
To calibrate the phantom modulus, eight pieces of homogeneous phantoms, including all concentrations of gelatin, are made.Compressive measurements are performed to characterize the Young's modulus of each component (Instron 5965).The testing rate is 0.05 min −1 , and the total compressive strain is ~2%.In this region, the phantoms obey the linear stress-strain behaviour and the Young's moduli are the slopes of the stress-strain curves (Supplementary Fig. 29).

Characterization of the device
We measure the frequency-dependent electrical impedance and phase angle curves following the conventional methods to characterize the electromechanical coupling performance of the ultrasound transducers (Fig. 1b).A network analyser (Hewlett-Packard 4195A) is used to do the testing with a frequency range of 1 MHz to 5 MHz.The effective electromechanical coupling coefficient K eff evaluates the electromechanical conversion efficiency of the transducer, which can be derived on the basis of the resonant (f r ) and anti-resonant (f a ) frequencies of the transducers 39 : Crosstalk is characterized by calculating the ratio between the measured and reference voltages.A peak-to-peak voltage of 5 V under a sinusoid burst mode is applied by a function generator (Keithley 3390) to excite the element.The reference voltage signal received by neighbouring elements is collected with a frequency increment step of 0.2 MHz.The voltage ratio is transformed to the logarithmic scale.
Insertion loss is used to reflect the sonographic sensitivity of the transducer, and these two variables are inversely proportional.The measurement is done in water using a functional generator with an output impedance of 50 Ω and an oscilloscope (Rigol DS1104) in a 1 MΩ coupling mode.A tone burst of a sine wave from 1.5 MHz to 4.5 MHz is generated to excite the element, and the ultrasound echoes reflected from a quartz crystal are captured by the same element.After compensating the 1.9 dB energy loss of transmitting into the quartz crystal and 2.2 × 10 −4 dB mm −1 MHz −2 loss due to the attenuation in water 83 , the insertion loss can be calculated by where V r and V t are the receiving and transmitting voltages, respectively, d is the distance between the transducer and the quartz crystal, and f r is the resonant frequency 83 .The matching circuit substantially improves the received signal amplitude, yielding a high sonographic sensitivity of the transducer (Fig. 1d and Supplementary Fig. 20).
To characterize its waterproof performance, the device is put underwater for 2 weeks (Supplementary Fig. 21).Pulse-echo signals are collected every day to analyse the variation in the sonographic signal-to-noise ratio.During the measurements, the transmitting distance and the ultrasound incident angle are kept consistent all the time to ensure the accuracy of the testing results.
The time-and frequency-domain characterizations of pulse-echo response and bandwidth are shown in Supplementary Fig. 19b,c.An experimental setup, including a multi-channel control system (Verasonics 256) and an aluminium block, is used to obtain the raw data.We put the aluminium block into water and ultrasonic waves are transmitted and received by the same transducer element.After getting the radiofrequency signal, a fast Fourier transform is applied to convert the signals from the time domain to the frequency domain.The frequency bandwidth of the signals at −6 dB is determined by where f u is the upper frequency, f l is the lower frequency and f c is the central frequency 84 .The signal-to-noise ratio of the transducer element can be calculated from Supplementary Fig. 19b, which is around 39 dB.Such a high signal-to-noise ratio is mainly attributed to the superior electromechanical coupling properties of the 1-3 composite, the low-temperature bonding technique and the electrical matching circuit.The bandwidth of the element is around 50% (Supplementary Fig. 19c), which still has room to improve.The relatively thin backing layer cannot completely dampen the excessive vibrations of the piezoelectric materials, leading to a long spatial pulse length.However, a thick backing layer will compromise the mechanical compliance of the device.Therefore, future improvements will focus on the development of new backing materials with high acoustic attenuation coefficients 85,86 while maintaining a small thickness of the backing layer.
The spatial resolution is characterized on the basis of the PSF of the transitional edges of an inclusion phantom 25 .The modulus distribution curves along the lateral and axial transition edges (that is, edges between the inclusion and matrix that are perpendicular and parallel to the ultrasound propagation direction, respectively) are extracted from the reconstructed 2D modulus mapping image (Supplementary Fig. 28).Then we apply the first derivative of the modulus distribution curves to derive the PSF, whose FWHM is defined as the spatial resolution.

Simulation of different imaging modes
The MATLAB software allows simulation of the three different ultrasound imaging modes: coherent plane-wave compounding, mono-focus and single plane wave.A model of 2 cm by 5 cm is built up for the simulations, where a bilayer structure with random points filling the entire region is set up.Here 1.5% and 0.5% of strain are added to the top and bottom layers with different moduli, respectively, and a new model with a compressed configuration is generated.Then, ultrasound waves are excited using the three transmission modes.Corresponding radiofrequency signals from pre-and post-compression are collected.The toolbox of MATLAB, Field II, is used to simulate the radiofrequency signals.These radiofrequency signals are beamformed to enhance the signal-to-noise ratio.Displacement fields are derived using the normalized cross-correlation algorithm.Finally, the spatial strain is mapped using the least-squares strain estimator (Supplementary Figs.7 and 24).

System setup and data collection
Our imaging system is mainly composed of two parts: the front end and the back end.The stretchable ultrasonic array as the front-end device transmits and receives ultrasound waves.The back-end device is a commercial multi-channel controller (Verasonics Vantage 256 system) that can be programmed to generate an arbitrary waveform.Cables are used to connect both ends for power supplying and data transmission.The MATLAB software (MathWorks) is used to code to control the back-end hardware and drive the front-end stretchable patch, involving the processes of plane wave compounding excitation, echo signals collection and storage, and data processing.As shown in Supplementary Fig. 7, two receive buffers are designed by MATLAB to collect and store the radiofrequency signals.The first and second receive buffers are used to store the radiofrequency signals before and after compression, respectively.Then these stored radiofrequency data are retrieved and processed to form the displacement fields.For the in vivo data collection, all human tests were approved under University of California, San Diego Institutional Review Board protocol 801085, which was obtained retrospectively due to the delays caused by coronavirus disease disruptions.

Phantom verifications by the commercial probe and ex vivo testing protocols
The commercial probe for verifying and reconstructing the 3D phantom images is a 64-element phased array transducer (P4-2v, Verasonics) with 2.8 MHz resonant frequency that is almost on a par with that of the stretchable ultrasonic array (f r = 3 MHz).During the experiments, a 3D linear stage (Newport) is used to fix the commercial ultrasound probe, apply pressure to get 2D strain images and move to the next position.Each moving step is 0.8 mm, the same as the elevational pitch of the stretchable patch.Sixteen slices of the 2D images are post-processed to reconstruct the 3D images (Supplementary Figs.34b and 36b).3D rendering is accomplished with the Amira software (Visualization Sciences Group).
All ex vivo studies are approved by the Institutional Animal Care and Use Committee at the University of California, San Diego.A frozen food-grade porcine abdominal tissue is used for ex vivo testing by the stretchable ultrasonic array and MRE.The porcine abdominal tissue is selected because of the well-defined anatomy and modulus distributions 61 .When the device is in conformal contact with the tissue surface, ~0.5-1% strain is applied, which can be followed in real time by a customized graphical user interface.If needed, the large strain with high NCCs can potentially be achieved by combining successive steps of small strain 87 .Then, 16 slices of 2D displacement fields are collected sequentially under the compressive strain.The inverse elasticity problem model yields the modulus mapping images when the predicted displacement fields match the measured ones (Supplementary Fig. 41).More details of the inverse elasticity problem calculation are in Supplementary Note 10.
The MRE test (General Electric Discovery MR750 3.0T) is performed right after the ultrasound test.A hermetic bag is used to preserve the porcine abdominal tissue during travelling and testing to prevent any dehydration-induced changes in the modulus.A mechanical vibration paddle is tightly bonded to the sample tissue to generate shear waves in the tissue.Fifteen scanning slices with 0.9 mm × 0.9 mm spatial resolutions of each are obtained.The total measuring time is about 30 min.The shear modulus of each slice can be displayed in ImageJ software.Specifically, Young's modulus equals to three times of shear modulus in soft tissues 58 .The modulus ratios of all slices are then calculated to compare with the results from the stretchable ultrasonic arrays.The Amira software is applied to combine all slices (reconstructed in y and z directions) of ultrasound elastographic images and MRE results along the x direction to generate the 3D volumetric image.Cubic spline interpolation is used to compensate and smooth the image in the x direction.
The possible sacrifice of contrast due to the variation regularization and the 2D finite element model may have resulted in the slightly higher modulus contrast of MRE than ultrasound elastography.Besides, the 2D elasticity model used in solving the inverse elasticity problem may introduce unavoidable discrepancies when applied to cross-sections of anisotropic 3D samples 60 .Additionally, different imaging equipment manufacturers using proprietary data processing algorithms may also bring systematic deviations 88 .

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Fig. 1 |
Fig. 1 | Working principle, design and fabrication of the stretchable ultrasonic array.a, Schematics of the stretchable ultrasonic array laminated on a soft tissue.The device consists of a 16-by-16 array of transducer elements connected in-parallel by a seven-layer electrode and encapsulated with waterproof and biocompatible silicone elastomer.The exploded view of each element is in Supplementary Fig. 6, showing the structure and components of the element.This layout allows activation of each element individually with a particular time-delay profile and capture of the reflected echoes from scattering sources, pre-(left) and post-compression (right).The beamformed pre-and post-compression signals are cross-correlated to derive the displacement,

FittingFig. 2 |
Fig. 2 | Strategies for elastography.a, Simulations showing the SNR e and CNR e of the coherent plane wave compounding, mono-focus and single plane wave transmission modes in a bilayer phantom.These data are extracted at seven different depths along the central line of strain images in Supplementary Fig. 24.b, The SNR e as a function of the step size of the steering angle.Inset is a schematic diagram showing the process of compounding: a set of plane waves are transmitted with different steering angles and a final compounded image is built by coherently adding all beamformed images.c, The SNR e and reconstruction time with different numbers of steering angles.d, The SNR e and NCC as a function of applied strain.In NCC, data were normalized by the product of the standard deviations of pre-and post-compression signals.The shaded area with >0.8 NCC by a −3 dB strain filter shows the dynamic range.e, Lateral and axial resolutions of the stretchable ultrasonic array based on the FWHMs of PSFs.f, Quantification of the contrast resolution based on the relationship between the CNR e and modulus contrast of phantoms.

Fig. 3 |
Fig. 3 | Characterizations on phantom models.a-d, Four types of phantom are used to simulate different pathological tissue environments: a bilayer phantom where the top layer is 2.03 times stiffer than the bottom layer (a), a cyst phantom consisting of a cylindrical cavity filled with fluid (b), phantoms with a cylindrical inclusion that is 1.61 times stiffer than the matrix (c) and a spherical mass that is 2.54 times stiffer than the matrix (d).B-mode images of all phantoms acquired

Fig. 4 |
Fig. 4 | Validation and serial monitoring on biological specimens.a, A 3D quantitative elastographic image of a porcine abdominal tissue by the stretchable ultrasonic array.In the tissue, two muscle layers intersect with the intermuscular fat, and the fascias are embedded in the transversus abdominis muscle.The 3D elastogram embodies homogeneous, lower modulus fat layers made of glycerol and fatty acids molecules, and relatively heterogeneous muscle groups with higher modulus made from criss-crossing and confluence fibres.b, The average 3% discrepancy of the 16 pairs of measured and predicted displacement fields, with a high degree of correspondence between the

Fig. 5 |
Fig. 5 | Multi-site mapping and serial surveillance of delayed-onset muscle soreness in human.a-d, Optical images, B-mode images and corresponding strain mapping results of a lateral shoulder joint (a), an anterior forearm (b), an anterior thigh (c) and a posterior calf (d).Key anatomical structures are labelled in the strain images.e, Schematics showing the exercise and physiotherapy protocols.f, Serial monitoring results of normalized modulus contrast and soreness intensity of the biceps brachii muscle before and after the eccentric exercise.We map the modulus distribution of the muscle and then calculate the Articlehttps://doi.org/10.1038/s41551-023-01038-w Articlehttps://doi.org/10.1038/s41551-023-01038-wPublisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.Extended Data Fig. 1 | Simulation results of uniform and nonuniform compressions.(a) A synthetic specimen with an inclusion that has a shear modulus 10 times higher than that of the surrounding matrix.The synthetic "measured" displacement fields based on (b) uniform compression and (c) non-uniform compression on the specimen.(d) and (e) show the strain distributions from uniform and non-uniform compression, respectively.The strain distributions vary upon different applied loads, which indicates that the strain-based elastography is qualitative only reflecting a relative stiffness of each component.The reconstructed modulus distributions obtained from solving inverse elasticity problems based on (f) uniform compression and (g) non-uniform compression.(h) Quantitative analysis of modulus contrasts and their deviations from the ground truth in (a).Extended Data Fig. 3 | Serial surveillance of delayed-onset muscle soreness in multiple subjects.Serial monitoring results of normalized modulus contrast of the biceps brachii muscle before and after the eccentric exercise.All six subjects did two rounds of experiments.In the first round, all subjects took the natural recovery after doing the exercise.In the second round, subjects in (a), (b) took natural recovery, subjects in (c), (d) took massotherapy, and subjects in (e), (f) took hyperthermia.For each test, we mapped the modulus distribution of the muscle.Then, we calculated the mean and standard deviations of the biceps brachii area.In each physiotherapy session, data are normalized by the modulus contrast of the biceps brachii muscle before exercise.Extended Data Fig. 4 | Validation by a clinical standard method on a larger sample size.(a) Changes in the maximal voluntary contraction torque (black) and normalized modulus contrast (blue) of the biceps brachii muscle measured by the stretchable ultrasonic array before and after the eccentric exercise of 16 subjects.The points and error bars of the blue curves indicate the mean and standard deviation of modulus contrast of the biceps brachii muscle of every test.For each test, we mapped the modulus distribution of the muscle.Then, we calculated the mean and standard deviations of the biceps brachii area.In each physiotherapy session, data were normalized by the modulus contrast of the biceps brachii muscle before exercise.As a clinical standard approach, the maximal voluntary contraction can validate the muscle strength and muscle damage.(b) P-values of the clinical standard method and the stretchable ultrasonic array, which were calculated by single-sided paired t-test.Error bars are standard deviations of the data of 16 subjects (n = 16).