Title

题目

Measurement of biomechanical properties of transversely isotropic  biological tissue using traveling wave expansion

利用行波展开法测量横观各向同性生物组织的生物力学特性

01

文献速递介绍

纤维嵌入结构在自然界中普遍存在。从脑白质（罗曼诺等人，2012年）到肌肉（帕帕佐格鲁等人，2006年）等生物组织，都可以被视为嵌入了纤维的基质，这对于神经信号传导和肌肉力量至关重要。嵌入的纤维赋予了生物组织各向异性的生物力学特性，这对于全面理解其生理作用以及疾病的影响至关重要（富兰克林和弗伦奇-康斯坦特，2008年；萨尔布等人，2016年）。纤维嵌入结构的破坏，比如脑白质中轴突的髓鞘脱失，会损害生物和生理功能。研究表明，改善少突胶质细胞的健康状况和髓鞘的完整性，有望成为延缓阿尔茨海默病发展和减缓其进程的一个有前景的靶点（德普等人，2023年）。准确测量各向异性有助于深入了解形态与功能之间的复杂相互作用，为个性化治疗方案提供参考，并推动医学和生物工程解决方案的创新（卢等人，2023年；孙等人，2024年）。此外，在体内对组织各向异性进行表征，对于监测疾病进展和评估治疗效果至关重要（卡彭特等人，2015年；米托等人，2018年；邵等人，2019年；孙等人，2024年）。它在表征生物材料和组织工程方面也起着关键作用（宋等人，2023年；于等人，2019年；张等人，2023年），其目标是复制或改善天然组织的力学性能（加尔布塞拉等人，2008年；李和戈尔，2004年）。因此，对纤维嵌入组织的生物力学特性进行表征，对于理解组织功能以及指导开发有效的治疗方法和干预措施至关重要。 对复杂纤维组织内的各向异性粘弹性特性进行体内表征是一项重大挑战，这源于各向异性材料所固有的复杂力学行为。传统的基于波的弹性成像技术通过波速的变化来推断各向异性。典型的应用包括使用多角度布里渊显微镜（埃尔托尼等人，2022年）来确定角膜的各向异性特征，以及使用三维超声旋转剪切波成像（奈特等人，2022年）来评估肌肉组织的各向异性。然而，这些方法需要对一致的组织结构进行假设，并且要生成纯净的单色波场。这些对于测量具有复杂纤维网络的深部组织（如脑白质）构成了挑战。磁共振弹性成像（MRE）（冯等人，2023年；加泰泽尔等人，2012年；森坂等人，2018年；施等人，2022年）能够捕捉三维波场，其优势在于可以适应大规模且复杂的组织结构的异质性。尽管如此，波场中包含了多种极化方式和波速。从丰富的波场数据中准确、高效地提取各向异性粘弹性信息仍然具有挑战性（安德森等人，2016年）。 由于大多数具有纤维嵌入结构的软组织（如脑白质）具有明显的纤维取向，通常会采用横观各向同性（TI）材料模型。一般的横观各向同性模型由五个参数组成（古尔德和冯，2018年）。通过假设生物组织的不可压缩性，应用最广泛的模型是具有三个参数的近似不可压缩横观各向同性（NITI）模型（冯等人，2013年；特韦滕等人，2015年）。许多研究已经利用NITI模型对纤维嵌入材料的参数估计进行了探讨（特韦滕等人，2017年）。方向滤波方法（卡拉等人，2019年；特韦滕等人，2015年）利用方向和极化滤波器来分离波场、估计波速并合成弹性参数，然而其在异质组织估计中的性能受到滤波器内核的限制（罗曼诺等人，2012年、2005年；特韦滕等人，2017年）。直接代数反演方法（郭等人，2016年；辛库斯等人，2005年；宋等人，2012年）通常假设纤维取向与坐标对齐以简化计算，这可能不适用于复杂的纤维结构。机器学习算法也已被应用（侯等人，2022年），但仅限于由声辐射力引起的局部波场。仅依靠波前形状进行估计的方法（帕帕佐格鲁等人，2006年、2005年）仅适用于简单的波场。采用有限元方法的迭代反演（巴巴埃等人，2021年；乔蒂等人，2023年；麦加里等人，2022年、2021年）可以缓解复杂模型中的精度挑战。然而，对整个大脑进行反演的时间通常超过12小时，并且参数调整也很复杂。目前仍然需要一种能够在模型精度和计算效率之间取得平衡，用于估计纤维嵌入材料各向异性粘弹性特性的反演方法。 在这项研究中，我们引入了一种基于行波展开（TWE）的近似不可压缩横观各向同性（TWE-NITI）方法，以基于NITI模型估计软生物组织的各向异性生物力学特性。这种方法将行波展开的解析优势与波的特性相结合，用于直接且高效地进行NITI反演，同时保持模型精度。通过模拟实验、离体肌肉实验以及在体人脑实验的验证，证明了TWE-NITI方法在体内各向异性生物力学分析方面的潜在价值。

Abatract

摘要

The anisotropic mechanical properties of fiber-embedded biological tissues are essential for understanding their development, aging, disease progression, and response to therapy. However, accurate and fast assessment of mechanical anisotropy in vivo using elastography remains challenging. To address the dilemma of achieving both accuracy and efficiency in this inverse problem involving complex wave equations, we propose a computational framework that utilizes the traveling wave expansion model. This framework leverages the unique wave characteristics of transversely isotropic material and physically meaningful operator combinations. The analytical solutions for inversion are derived and engineering optimization is made to adapt to actual scenarios. Measurement results using simulations, ex vivo muscle tissue, and in vivo human white matter validate the framework in determining in vivo anisotropic biomechanical properties, highlighting its potential for measurement of a variety of fiber-embedded biological tissues.

嵌入纤维的生物组织的各向异性力学特性对于理解其生长发育、衰老过程、疾病进展以及对治疗的反应至关重要。然而，使用弹性成像技术在体内准确且快速地评估力学各向异性仍然颇具挑战性。 为了解决在这个涉及复杂波动方程的反问题中同时实现准确性和效率的两难困境，我们提出了一个利用行波展开模型的计算框架。该框架利用了横观各向同性材料独特的波动特性以及具有物理意义的算子组合。我们推导出了用于反演的解析解，并进行了工程优化以适应实际情况。 通过模拟实验、离体肌肉组织实验以及在体人类白质的测量结果，验证了该框架在确定体内各向异性生物力学特性方面的有效性，突显了其在测量各种嵌入纤维的生物组织方面的潜力。

Method

方法

2.1. Wave characteristics in NITI model

For materials embedded with fibers, anisotropic properties are observed relative to the fiber orientation, namely transversely isotropy. To describe the mechanical behaviour of the biological structures with this feature, a nearly incompressible transversely isotropic (NITI) model (Feng et al., 2013), in terms of substrate shear modulus μ, shear anisotropy ϕ and tensile anisotropy ζ, can be used (Fig. 1a). The polarization and wave speed of the mixed wave field of transversely isotropic biological tissue have complex relationships with the anisotropy parameters, wave propagation direction and fiber direction (Fig. 1b, Supplementary Note S1).If the fiber orientation is a = [a1, a2, a3] T, wave propagation direction is n= [n1, n2, n3] T, the angle between a and nis θ, the slow and fast shear waves of the NITI model can be obtained from the Kelvin–Christoffel equations generated from Eq. (S4) (Carcione, 2022; Tweten et al., 2015). If the slow shear wave speed is cs, the fast shear wave speed is cf , the polarization direction of the slow shear wave is ms = [ms1 , ms2 , ms3 ] T and polarization direction of the fast shear wave is mf = [ mf1 , mf2 , mf3 ]T* (Fig. 1b), the following equations describe the relationship between the elastic parameters of the NITI model and the wave speed (Tweten et al., 2015): ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩c2 s = μρ{ 1 + ϕcos2θ }c*2 f = μ ρ { 1 + ϕcos22θ + ζsin2 2θ }ms = n × amf = n× ms(1)

2.1. NITI模型中的波特性 对于嵌入纤维的材料，相对于纤维取向会表现出各向异性特性，也就是横观各向同性。为了描述具有这种特征的生物结构的力学行为，可以使用近似不可压缩横观各向同性（NITI）模型（冯等人，2013年），该模型涉及基质剪切模量μ、剪切各向异性ϕ和拉伸各向异性ζ（图1a）。横观各向同性生物组织的混合波场的极化方向和波速与各向异性参数、波传播方向以及纤维方向之间存在复杂的关系（图1b，补充说明S1）。 如果纤维取向为a = [a1, a2, a3]T，波传播方向为n = [n1, n2, n3]T，a与n之间的夹角为θ，那么NITI模型的慢剪切波和快剪切波可以从由公式（S4）推导出来的开尔文 - 克里斯托费尔方程中得到（卡乔内，2022年；特韦滕等人，2015年）。如果慢剪切波速度为cs，快剪切波速度为cf，慢剪切波的极化方向为ms = [ms1, ms2, ms3]T，快剪切波的极化方向为mf = [mf1, mf2, mf3]T（图1b），那么以下公式描述了NITI模型的弹性参数与波速之间的关系（特韦滕等人，2015年）： [ \begin{cases} c_s^2 = \frac{\mu}{\rho}(1 + \phi\cos^2\theta) \ c_f^2 = \frac{\mu}{\rho}(1 + \phi\cos^22\theta + \zeta\sin^2 2\theta) \ \mathbf{m}s = \mathbf{n} \times \mathbf{a} \ \mathbf{m}f = \mathbf{n} \times \mathbf{m}_s \end{cases} ]  (1)

Conclusion

结论

In this study, we proposed a computational framework for solving the inverse problem of magnetic resonance elastography in anisotropic biological tissue using the Traveling Wave Expansion (TWE) model. The proposed TWE-NITI effectively resolves the contradiction between computational complexity and accuracy in estimating complex anisotropic parameters. The performance of the method has been validated through simulations, ex vivo muscle experiments, and in vivo human brain experiments. The results demonstrate its potential for applications in diagnosis and tissue engineering.

在本研究中，我们提出了一个计算框架，用于运用行波展开（TWE）模型来解决各向异性生物组织中磁共振弹性成像的反问题。所提出的TWE-NITI方法有效地解决了在估计复杂各向异性参数时计算复杂度与准确性之间的矛盾。 该方法的性能已通过模拟实验、离体肌肉实验以及在体人脑实验得到了验证。结果表明，它在疾病诊断和组织工程等方面具有应用潜力。

Results

结果

3.1. Validation of anisotropic inversion using simulated wavefield patches

A dataset simulated using the TWE model (Details in 2.5) was employed to evaluate the performance of the TWE-NITI method in estimating anisotropic parameters. Within a 4 × 4 × 4 cubic patch with a resolution of 1 mm, a series of single-frequency (30Hz) wavefields containing different numbers (M = 1,2,3) of traveling waves in random propagation directions were generated (Fig. 2a). These simulations were designed to systematically vary the fiber orientations, defined by an orientation vector function with θ1 and θ2 (Fig. 2b). Specifically, fiberorientations were first sampled with a fixed θ1 (60) and varying θ2 (trajectory 1, yellow), and then with a fixed θ2 (60) and varying θ1 (trajectory 2, gray). Each case was sampled with 200 points for simulation purposes.Subsequently, a comprehensive examination and assessment of the inversion performance using the TWE-NITI method was conducted. If the wavefield with only one wave propagation direction M = 1 was used, the estimation results of the anisotropic parameters were not satisfactory (Fig. 2c). However, with 2 or 3 different directions of traveling waves M = 2 and M = 3, the relative errors of parameter estimation reduced to 4.6% and 3.3%, respectively. This demonstrates that with two or more different traveling wave directions (M ≥ 2), the TWENITI method is capable of accurately recovering all of the anisotropic, viscoelastic parameters of the NITI model, using wavefields at a one single frequency. While inversion can indeed be performed using a single frequency, multifrequency approaches are essential to enhance the signal-to-noise ratio (SNR) and improve robustness to noise in practical applications, as will be discussed in the subsequent sections. These results are in agreement of the minimum requirements(Tweten et al., 2017) for NITI inversion, showing the potential of using TWE-NITI method for anisotropic inversion. Inversion verification has also been conducted using the wavefield generated using COMSOL (Supplementary Note S4, Figure. S4). This ensures the robustness and reliability of TWE-NITI.

使用基于行波展开（TWE）模型模拟得到的数据集（详见2.5部分）来评估TWE-NITI方法在估计各向异性参数方面的性能。在一个边长为4×4×4毫米、分辨率为1毫米的立方体区域内，生成了一系列单频（30赫兹）波场，这些波场包含数量不同（M = 1, 2, 3）、传播方向随机的行波（图2a）。 这些模拟旨在系统地改变由θ1和θ2定义的取向向量函数所确定的纤维取向（图2b）。具体来说，首先在固定θ1（取值为60）且θ2变化的情况下对纤维取向进行采样（轨迹1，黄色），然后在固定θ2（取值为60）且θ1变化的情况下对纤维取向进行采样（轨迹2，灰色）。为了进行模拟，每种情况都采样了200个点。 随后，对使用TWE-NITI方法的反演性能进行了全面的检查和评估。如果仅使用具有一个波传播方向（M = 1）的波场，各向异性参数的估计结果并不理想（图2c）。然而，当使用具有2个或3个不同传播方向的行波（M = 2和M = 3）时，参数估计的相对误差分别降低到了4.6%和3.3%。 这表明，当有两个或更多不同的行波传播方向（M ≥ 2）时，TWE-NITI方法能够利用单一频率的波场准确地恢复NITI模型的所有各向异性粘弹性参数。 虽然确实可以使用单一频率进行反演，但在实际应用中，多频率方法对于提高信噪比（SNR）以及增强对噪声的鲁棒性至关重要，这将在后续章节中进行讨论。这些结果符合NITI反演的最低要求（特韦滕等人，2017年），显示了使用TWE-NITI方法进行各向异性反演的潜力。 还使用了通过COMSOL生成的波场进行了反演验证（补充说明S4，图S4）。这确保了TWE-NITI方法的稳健性和可靠性。

Figure

图

Fig. 1. An overview of TWE-NITI method for extracting anisotropic viscoelastic parameters. (a) The three parameters (μ, ϕ and ζ) of the NITI model characterize the anisotropic behaviour of fiber-embedded structures. (b) The polarization ( ms, mf) and wave lengths ( λ**s, λ**f ) of the slow and fast waves within the NITI medium exhibit a high degree of complexity, which is closely associated with the propagation direction n, the fiber orientation a, and the NITI parameters μ, ϕ, and ζ. This sub-figure was adapted from Fig. 2 of (Tweten et al., 2015). (c) The TWE model facilitates the decomposition of the complex NITI wave field U into slow and fast wave components ( Us, Uf ) , each associated with a single direction of propagation, thereby simplifying the analysis process. (d) Operators with specific physical significance are applied to all traveling wave components. (e) The TWE-NITI method, based on the analysis of wavefield using the TWE model, establishes a relationship between the wavefield, fiber orientations, and the anisotropic parameters of the NITI model. And this method is capable of estimating the anisotropic parameters of the NITI model from the wavefield of fiber-embedded structures acquired by MRE. The red arrows indicate the design process of the computational method, and the blue indicate the application process of the computational method

图1. 用于提取各向异性粘弹性参数的TWE-NITI方法概述  (a) NITI模型的三个参数（μ、ϕ 和 ζ）表征了纤维嵌入结构的各向异性行为。(b) NITI介质中慢波和快波的极化方向（ms，mf）以及波长（λs，λf）表现出高度的复杂性，这与传播方向 n、纤维取向 a 以及NITI参数 μ、ϕ 和 ζ 密切相关。此子图改编自（特韦滕等人，2015年）的图2。(c) TWE模型有助于将复杂的NITI波场 U 分解为慢波和快波分量（Us，Uf），每个分量都与单一的传播方向相关联，从而简化了分析过程。(d) 具有特定物理意义的算子被应用于所有行波分量。(e) 基于使用TWE模型对波场的分析，TWE-NITI方法建立了波场、纤维取向和NITI模型的各向异性参数之间的关系。并且该方法能够从通过磁共振弹性成像（MRE）获取的纤维嵌入结构的波场中估计NITI模型的各向异性参数。红色箭头表示计算方法的设计过程，蓝色箭头表示计算方法的应用过程。

Fig. 2. Performance of the TWE-NITI method in estimating the viscoelastic parameters. (a) Three types of wavefields were generated in a 4 × 4 × 4 cubic patch with a resolution of 1 mm at 30 Hz, which contain different numbers (M=1,2,3) of traveling waves in random propagation directions. (b) Fibers were simulated with an orientation vector function defined by θ1 and θ2. The fiber orientations were first sampled with a fixed θ1 (60) and varying θ2 (trajectory 1, yellow). Then, the fiber orientations were sampled with a fixed θ2 (60) and varying θ1 (trajectory 2, gray). Each case was sampled with 200 points for simulation. (c) Comparisons of the anisotropic viscoelastic parameters estimated using TWE-NITI method and the ground truth. Plots were distinguished by trajectories of fiber orientations and M, with real (Re(⋅)) and imaginary (Im(⋅)) parts separated for clarity. The angle of the polar coordinates corresponds to the angle defined in (b). Radial unit of the polar coordinates is in kPa.

图2. TWE-NITI方法在估计粘弹性参数方面的性能  (a) 在一个边长为4×4×4毫米、分辨率为1毫米的立方体区域内，在30赫兹的频率下生成了三种类型的波场，这些波场包含数量不同（M* = 1, 2, 3）、传播方向随机的行波。(b) 使用由θ1和θ2定义的取向向量函数来模拟纤维。首先在固定θ1（取值为60）且θ2变化的情况下对纤维取向进行采样（轨迹1，黄色）。然后，在固定θ2（取值为60）且θ1变化的情况下对纤维取向进行采样（轨迹2，灰色）。每种情况都采样200个点用于模拟。(c) 将使用TWE-NITI方法估计得到的各向异性粘弹性参数与真实值进行比较。绘图根据纤维取向的轨迹和M值进行区分，为了清晰起见，将实部（Re(⋅)）和虚部（Im(⋅)）分开显示。极坐标的角度对应于(b)中所定义的角度。极坐标的径向单位为千帕（kPa）。

Fig. 3. Performance evaluation of the TWE-NITI method in estimating the viscoelastic properties of fiber-embedded phantom. (a) A cube (240 × 240 × 60 mm3 ) with an isotropic spherical core and radiated fibers was created. Shear excitations were applied to the sides of the cube. The central slice (pink coloured) was selected for demonstration. (b) Wavefields of the central slice at frequencies of 30, 40, 50, and 60Hz. (c) Comparisons of the ground truth and estimated viscoelastic parameters at the central slice. (d) Comparisons of the estimated viscoelastic parameters (mean and standard deviation) for regions with isotropic and anisotropic properties

图3. TWE-NITI方法在估计嵌入纤维的仿体粘弹性特性方面的性能评估  (a) 创建了一个尺寸为240×240×60立方毫米的立方体，其内部有一个各向同性的球形核心以及呈辐射状分布的纤维。对该立方体的各侧面施加剪切激励。选取了中间切片（粉红色部分）用于展示。(b) 中间切片在30赫兹、40赫兹、50赫兹和60赫兹频率下的波场情况。(c) 中间切片处真实的粘弹性参数与估计所得的粘弹性参数之间的比较。(d) 对具有各向同性和各向异性特性的区域所估计的粘弹性参数（均值和标准差）进行比较

Fig. 4. Estimation of anisotropic parameters of muscle tissue using benchtop and MRE test. (a) Experimental setup of muscle tissue testing with MRE. Two actuators were placed at two different locations of the muscle tissue to produce different wave propagation patterns. Arrows indicate to different vibration directions. (b) Benchtop test of muscle tissue. Two samples of muscle tissue were measured with shear and compression so that the viscoelastic properties were acquired with respect to parallel and perpendicular to the fiber orientation. (c) The muscle tissue used for testing and its corresponding T1-weighted (T1 W) image and fiber orientation map generated from diffusion tensor image (DTI). ROIs were selected to include the two samples used for rheological tests. Displacement wavefields are shown for the three directions at 30 Hz. (d) Estimated viscoelastic parameters from MRE. (e) Comparisons of the mean and standard deviation values of the anisotropic parameters measured from MRE and benchtop test.

图4. 使用台式测试和磁共振弹性成像（MRE）测试估计肌肉组织的各向异性参数  (a) 利用磁共振弹性成像对肌肉组织进行测试的实验装置。将两个激振器放置在肌肉组织的两个不同位置，以产生不同的波传播模式。箭头表示不同的振动方向。(b) 肌肉组织的台式测试。对两个肌肉组织样本进行了剪切和压缩测量，从而获取了平行和垂直于纤维取向的粘弹性特性。(c) 用于测试的肌肉组织及其相应的T1加权（T1 W）图像，以及由扩散张量图像（DTI）生成的纤维取向图。选择的感兴趣区域（ROI）包含了用于流变测试的两个样本。图中展示了在30赫兹频率下三个方向的位移波场。(d) 由磁共振弹性成像估计得到的粘弹性参数。(e) 比较由磁共振弹性成像和台式测试所测得的各向异性参数的均值和标准差数值。

Fig. 5. MRE-based estimation of anisotropic parameters of brain tissue in vivo. (a) Wavefields were acquired at two different testing modes. In Test 1, the actuator plate was closely fitted to the entire head of the volunteer, directly applying vibration to the brain. In Test 2, the actuator plate was positioned on the neck of the volunteer, with the neck acting as a conduit to transmit the vibration to the brain. The corresponding wavefields at 30 Hz are shown for the two tests. (b) Magnitude map of MRE, FA and fiber orientation map generated from DTI. (c) Real and imaginary parts of the anisotropic parameters μ, μϕ and μζ estimated by the two tests, and the fused results.

图5. 基于磁共振弹性成像（MRE）对体内脑组织各向异性参数的估计  (a) 在两种不同的测试模式下获取波场。在测试1中，激振器板紧密贴合在志愿者的整个头部，直接向大脑施加振动。在测试2中，激振器板放置在志愿者的颈部，颈部充当将振动传递到大脑的通道。图中展示了两种测试在30赫兹频率下相应的波场。(b) 磁共振弹性成像的幅值图、各向异性分数（FA）图以及由扩散张量成像（DTI）生成的纤维取向图。(c) 两种测试所估计出的各向异性参数μ、μϕ和μζ的实部和虚部，以及融合后的结果。

Fig. 6. Statistical estimation of anisotropy parameters in the brains of 10 healthy volunteers using TWE-NITI. The mean and standard deviation values of the NITI parameters for gray matter (GM) and white matter (WM), along with the experimental data points, are illustrated.

图6：使用TWE-NITI方法对10名健康志愿者大脑中的各向异性参数进行的统计估计。图中展示了灰质（GM）和白质（WM）的NITI模型参数的平均值和标准差数值，同时也给出了实验数据点。

Table

表

Table 1 The mean value of the NITI parameters measured from a single test and those obtained from the fusion of two tests across the entire brain. Units are in kPa.

表1 从单次测试中测得的NITI模型参数的平均值，以及通过对全脑的两次测试结果进行融合后得到的NITI模型参数的平均值。单位为千帕（kPa）。

Table 2 The mean and standard deviation of the NITI parameters for gray and white matter, obtained from 10 healthy volunteers. Units are in kPa

表2 从10名健康志愿者处获得的灰质和白质的NITI模型参数的平均值和标准差。单位为千帕（kPa）。