1. Part 4: Viscoelastic Properties of Soft • Tissues in a Living Body Measured by MR • Elastography Gen Nakamura • Department of Mathematics, Hokkaido University, Japan • (Supported by Japan Science and Technology Agency) • Joint work with Yu, Jiang • ICMAT, Madrid, May 12, 2011

2. Magnetic Resonance Elastography, MRE • A newly developed non-destructive technique • (Muthupillai et al., Science, 269, 1854-1857, 1995, Mayo Clinic.) • Measure the viscoelasticity of soft tissues in a living body • Diagnosis: • the stage of liver fibrosis • early stage cancer: breast cancer, pancreatic cancer, prostate cancer, etc. • neurological diseases: Alzheimer’s disease, hydrocephalus, multiple sclerosis, etc. • Nondestructive testing (high frequency rheometer): • biological material, polymer material

3. MRE System in Hokkaido Univ. (JST Proj.) Japan Science and Technology Agency (JST) (1) External vibration system Electromagnetic vibrator (2) Pulse sequence with motion-sensitizing gradients (MSG) GFRP Bar 2~4 m (3) phantom Object External vibration system Micro-MRI Wave image Storage modulus (4) Inversion algorithm

4. MRE phantom: agarose or PAAm gel 100mm 65mm • --- time harmonic external vibration (3D vector) • --- frequency of external vibration (50～250Hz) • --- amplitude of external vibration (≤ 500 μm ) 70mm 10 mm hard soft

5. Viscoelastic wave in soft tissues • Time harmonic external vibration • Interior viscoelastic wave • --- amplitude of viscoelastic wave ( : real part, : imaginary part) viscoelastic body after some time

6. MRE measurements: phase image MRIsignal real part: R imaginary part: I ２D FFT magnitude image phase image MRE measurement

7. MRE measurements: phase image components in vertical direction (unit: )

8. Data analysis for MRE viscoelasticity of soft tissues or phantom interior wave displacement Step 1: modeling Step 3: recovery （inverse problem）　　 Step 2: numerical simulation （forward problem）　 viscoelasticity models for soft tissues or phantom (PDE)

9. Viscoelasticity models for soft tissues • Time: • : bounded domain; • : Lipschitz continuous boundary; • Displacement: • General linear viscoelasticity model:

10. Viscoelasticity models for soft tissues • Stress tensor: • Density: • Small deformation (micro meter) ⇒ linear strain tensor

11. Constitutive equation • Voigt model: • Maxwell model: • Zener model:

12. Viscoelasticity tensors • full symmetries: • strong convexity (symmetric matrix ):

13. Time harmonic wave • Boundary: • : open subsets of with , Lipschitz continuous; • Time harmonic boundary input and initial condition: • Time harmonic wave (exponential decay): • Jiang, et. al., submitted to SIAM appl. math.. (isotropic, Voigt) • Rivera, Quar. Appl. Math., 3(4), 629-648, 1994. • Rivera, et. al., Comm. Math. Phys. 177(3), 583–602, 1996.

14. Time harmonic wave • Stationary model: • Sobolev spaces of fractional order 1/2 or 3/2 • an open subset with a boundary away from and • the set of distributions in the usual fractional Sobolev space compactly supported in • This can be naturally imbedded into

15. Constitutive equation (stationary case) • Voigt model: • Maxwell model: • Zener model:

16. Modified Stokes model • Isotropic+ nearly incompressible • Asymptotic analysis ⇒ modified Stokes model: • Jiang et. al., Asymptotic analysis for MRE, submitted to SIAP • H. Ammari et. al., Quar. Appl. Math., 2008: isotopic constant elasticity

17. Storage modulus andloss modulus • Storage ・ loss modulus () • Voigt model • Maxwell model • Zener • Angular frequency: • Shear modulus: • Shear viscosity: • Measured by rheometer

18. Modified Stokes model • 2D numerical simulation (Freefem++) • Plane strain assumption mm

19. Curl operator • Modified Stokes model: • Constants : Curl operator: filter of the pressure term

20. Pre-treatment: denoising • Mollifier(Murio, D. A.: Mollification and Space Marching) • Smooth function defined in the nbd of • : a bounded domain • : an extension of to • Function : a nonnegative function over such that and

21. Denoising

22. Recovery of storage modulus • Constants: • Mollification: • Curl operator: • Numerical differentiation method • Numerical differentiation is an ill-posed problem • Numerical differentiation with Tikhonov regularization Unstable!!!

23. Recovery of storage modulus • Constants: • Mollification: • Curl operator: • Numerical Integration Method • : test region (2D or 3D) • : test function Unstable!!!

24. Recovery of storage modulus • Constants: • Mollification: • Curl operator: • Modified Integral Method • : test region (2D or 3D) • test region size: about one wavelength

25. Recovery from no noise simulated data Inclusion: small large outside • Exact value: 3.3 kPa 3.3 kPa 7.4 kPa • Mean value: 3.787 kPa 3.768 kPa 7.436 kPa • Stddev: 0.147 0.060 0.003 • Relative error: 0.1476 0.1418 0.00049

26. Recovery from noisy simulated data • 10% relative error Inclusion: small large outside • Exact value: 3.3 kPa 3.3 kPa 7.4 kPa • Mean value: 4.636 kPa 3.890 kPa 7.422 kPa • StdDev: 0.328 0.129 0.322 • Relative error: 0.4048 0.1788 0.00294

27. Recovery from experimental data Layered PAAm gel: hard (left) soft (right) Mean value: 31.100 kPa 10.762 kPa StdDev: 0.535 0.201 250 Hz 0.3 mm kPa cm

28. Recovery of storage modulus G’ Layered PAAm gel: hard (left) soft (right) Mean value: 31.100 (25.974) kPa 10.762 (8.988) kPa Standard deviation: 0.535 (6.982) 0.201 (4.407) modified method (old method (polynomial test function)) 250 Hz 0.3 mm kPa kPa cm cm

29. Recovery of storage modulus G’ Independent of frequencies (1 ~ 250 Hz) hard soft Rheometer: 32.5456 kPa 9.2472 kPa MRE, 250 Hz: 31.100 kPa 10.762 kPa Relative error: 0.0444 0.1638 Rheometer : ARES-2KFRT, TA Instruments Frequency: 0.1 ~ 10 Hz Strain mode: 5%

30. Thank you for your attentions!