Optical microelastography via a 2D boundary condition-free nonlinear inversion approach

The mechanical phenotype of a cell, including its viscoelastic properties, is recognized as a label-free biomarker for diagnosing cellular states. Optical microelastography (OME) assesses intracellular mechanical heterogeneity by mapping the shear modulus distribution within cells using time-harmonic elastic waves observed within an optical image plane. However, reconstructing viscoelastic properties at the microscale is challenging not only because of inherent scale limitations, but also because, in OME, the complex 3D wave motion is only tracked within a single 2D plane. To address this challenge, a 2D boundary-condition-free nonlinear inversion (2D-NoBC-NLI) method is introduced to reconstruct viscoelastic properties from noisy 2D displacement fields. Numerical simulations of a homogeneous sphere, a heterogeneous sphere, and an asymmetric nucleated cell were designed to assess the robustness of 2D-NoBC-NLI reconstructions. Experiments were conducted on homogeneous, 75 µm-diameter polyacrylamide (PAAm) microbeads, which were expected to yield uniform viscoelasticity maps. With optimum parameter conditions, the proposed 2D-NoBC-NLI approach achieved mean reconstruction errors ranging from 1 to 13% across all simulated models. Within homogeneous PAAm microbeads, the method demonstrated frequency dependency of viscoelastic parameters consistent with previous measurements. The proposed nonlinear inversion algorithm enables storage and loss moduli imaging without out-of-plane motion data, and without using simplifying 2D approximations. This technique supports 2D elastography imaging and may enable OME-based cell mechanobiology studies through spatially resolved viscoelastic property mapping.

Highlights

•Introduce a new non-linear inversion for viscoelasticity imaging in the context of optical microelastography without out-of-plane data.

•Numerical simulations show viscoelastic reconstructions for different geometries and noise conditions.

•Experimental feasibility demonstrated on 75 µm viscoelastic microbeads at multiple frequencies.

•The technique represents a new tool for viscoelasticity imaging in the context of cell mechanobiology.

1 Introduction

The mechanical phenotype of a cell is increasingly recognized as a label-free biomarker that reflects its state, and can be utilized for both diagnostic and therapeutic purposes [1, 2]. Precise and reliable microscopic mechanical characterization of cells has become essential for a comprehensive understanding of biological behavior and disease progression [3, 4]. Several techniques have been developed to assess cell mechanics including atomic force microscopy, microfluidics technologies, and micropipette aspiration [5]. Each method is tailored to characterize the intracellular or surface mechanics at a local or cell wide level, with specific spatial and temporal parameters to consider [6]. Consequently, measured values can vary significantly—often by several orders of magnitude—due to factors such as deformation rates, applied mechanical stress, and the specific measurement technique and tools employed, including probe geometry, contact area, and length scale, all of which impact viscoelastic property estimation [7]. Moreover, the time resolution at which measurements are taken can also influence the estimated properties [8].

Because cellular mechanics are inherently time-dependent, viscoelastic characterization capturing both storage and loss moduli is required to fully describe mechanical cell behavior. These parameters not only reflect underlying cytoskeletal remodeling but are also directly linked to fundamental processes including migration, differentiation, and disease progression [9]. Compared to elasticity alone, a full viscoelastic assessment provides superior characterization, with impact to better understand cancer metastasis and tissue fibrosis with a mechanical perspective, underscoring its potential as a practical biomarker for diagnosis, prognosis, and therapeutic response [10, 11]. For example, Ma et al. demonstrated that the distinction between normal and senescent endothelial cells can be clearly established using viscoelastic parameters, which provide greater sensitivity than elasticity measurements alone [12]. In addition, cells display pronounced intracellular heterogeneity in their mechanical behavior. The cytoplasm, nucleus, and cortical regions each exhibit distinct viscoelastic properties arising from differences in cytoskeletal organization, nuclear structure, and intracellular crowding [13–15]. Therefore, mapping intracellular viscoelastic heterogeneity provides a mechanical fingerprint of the cell’s state, which is indispensable for advancing our understanding of cellular biomechanics and disease mechanisms.

Elastography is an imaging technique that retrieves mechanical properties of soft tissues through a reconstruction based on induced deformations. Traditionally applied at the macroscopic level, recent advances integrating high-resolution optical techniques have refined its spatial resolution to the microscale [16, 17], enabling the characterization of intracellular mechanical properties that were previously unobservable [18]. Optical microelastography (OME) is one such technique to assess intracellular mechanical properties, based on optical microscopy [19]. It utilizes high-frequency (typically 15–60 kHz) induced elastic waves that propagate through the cell and are detected by a high-speed camera (e.g., 100–300 kfps) integrated into an optical microscopy platform, enabling detailed characterization of intracellular viscoelasticity. Grasland-Mongrain et al. [19] demonstrated the application of OME for ultrafast imaging of cell elasticity, highlighting its potential for high-resolution mechanical mapping, whereas Flé et al. [20] used OME to measure viscoelastic properties of mouse oocytes, underscoring its potential utility in reproductive biology.

In the field of elastography, one of the key challenges is the reconstruction of viscoelastic properties from partial or corrupted displacement measurements. The subzone nonlinear inversion (NLI) reconstruction method has been widely used in magnetic resonance elastography (MRE) to address this challenge [21], and was recently adapted for cell OME [20]. The method involves dividing the region of interest into smaller overlapping subzones. Within each subzone, the NLI approach minimizes the difference between displacements obtained from finite element (FE) simulations and the experimental data. By combining reconstructed material properties across all subzones, the impact of measurement noise is reduced. However, the FE forward problem within each subzone requires boundary conditions, which are constructed from measurement data, allowing displacement data noise to influence the viscoelastic property reconstruction. Additionally, the two-dimensional (2D) nature of optical microscopy data under-represents the three-dimensional (3D) nature of the cell’s mechanical response to external vibration [22].

To redress these dimensional limitations, elastography reconstruction methods based on 2D displacement data typically rely on plane-strain or plane-stress approximations [23, 24], which are only valid under specific conditions. For instance, plane-strain is applicable in structures with a large out-of-plane dimension, or structures where the measurement plane corresponds to a plane of symmetry, or cases where the out-of-plane motion is restricted. Similarly, plane stress generally applies to thin plates. These approximations are rarely valid in biological contexts and their use can significantly impact reconstruction accuracy [25]. Errors and artifacts in 2D reconstructions arise not only from measurement noise but also from systematic inaccuracies of the underlying 2D approximations themselves, and have been shown to roughly 20% in previous studies [23, 24]. In the case of OME, cells generally have complex geometric and physical structures, with diverse asymmetric morphologies and heterogeneous mechanical properties, rendering 2D approximations inadequate. Consequently, the application of plane-strain and plane-stress approximations is inappropriate and can introduce substantial inaccuracies and systematic artefacts in viscoelastic property reconstruction.

Alternatively, recent advances in NLI have introduced boundary-condition-free (NoBC) approaches, which utilize coupled adjoint field (CAF) formulations to eliminate the need for known boundary conditions in each subzone when solving FE forward problems [26–28]. Eliminating the use of measurement data as boundary conditions in the forward FE problem provides a more robust property estimate, particularly in the presence of noise. Building on these CAF formulations, we propose and introduce the 2D-NoBC-NLI approach to achieve reliable OME reconstruction without 2D geometry approximations. This method uses a loosely enforced incompressibility assumption within the NoBC-NLI framework to reconstruct viscoelastic properties from 2D data using a fully 3D model. The approach is particularly effective for complex incompressible materials at the microscale, as it does not require prior knowledge of morphology or out-of-plane motion. Notably, it enables the reconstruction of the intricate geometry of cells from 2D displacement fields available in OME. Additionally, it is compatible with any 2D imaging modality, offering a versatile solution for analyzing complex geometries.

In this study, the 2D-NoBC-NLI framework was implemented in OME to reconstruct viscoelastic properties at the microscale by mapping the complex valued shear modulus using elastic wave-induced 2D displacements extracted from a 3D geometry. Validation was performed using numerical simulations on homogeneous, heterogeneous, and asymmetric microscale models to understand the behavior of the reconstruction process in various 3D configurations. The sensitivity of the method to NLI inversion parameters was further assessed using a homogeneous spherical model to evaluate their effect on shear modulus reconstructions, followed by tests on heterogeneous and asymmetric models. Noise was then introduced into the heterogeneous and asymmetric simulation models to examine the robustness of the reconstruction process under experimental conditions. Finally, the feasibility of the technique was demonstrated through experiments on 75-micron diameter homogeneous polyacrylamide (PAAm) microbeads across a range of actuation frequencies. Results demonstrate that the 2D-NoBC-NLI method provides assessment of viscoelastic properties at the microscale, confirming its feasibility in the context of OME.

2 Methods

2.1 Experimental set up (OME)

A micropipette, fabricated from borosilicate glass capillaries (World Precision Instruments, 1B100-6, United States), was tapered using a vertical pipette puller (David Kopf Instruments, DKI700C, United States), polished, and bent with a microforge (Narishige Scientific Instruments, MF-2, Japan). This micropipette was attached to a piezoelectric transducer (Thorlabs, PK2FQP1, United States) to generate elastic waves within polyacrylamide (PAAm) microbeads. Details of the bead preparation protocol are provided in Section 2.6 (Experimental validation). The piezoelectric transducer was connected to a voltage amplifier (Amplifier Research, 75A250, United States) and a signal generator (Agilent Technologies, 33250A, United States) to induce 20–60 kHz harmonic actuations. The micropipette slightly touched the bead surface to produce shear motions, as in [19, 20]. This frequency range was selected to provide actuation conditions compatible with 2D-NoBC-NLI reconstructions, as discussed in Section 2.6.

To stabilize the suspended microbead during testing, a second holding micropipette connected to a microinjector (Narishige Scientific Instruments, IM21) was employed. Both micropipettes were attached to micromanipulators to facilitate precise handling and positioning. The system was mounted on an inverted microscope (Olympus, IX71, Japan) equipped with a ×40 magnification objective lens. The microscope was coupled to an ultrafast camera (Photron Limited, Fastcam SA-Z, Japan), providing an effective resolution of 0.5 µm per pixel at an image acquisition rate of 250 kHz. An overview of the experimental setup is displayed in Figure 1a.

Figure 1

Figure 1. (a) Experimental set up. (b) Displacements calculated over time compared to the resting state. (c) 2D-NoBC-NLI reconstruction scheme (a.u., arbitrary units; equations are defined in Section 2.3; the wave propagation is perpendicular to the particle displacement shown).

In-plane wave displacement tracking was done using the Lagrangian speckle model estimator (LSME), which is based on the Lucas-Kanade optical flow method [29]. This approach estimates a 2D displacement field by comparing each frame in a time-series of images with an image of the microbead at rest, before the onset of vibration (Figure 1b). In-plane displacements were estimated in x and y directions, using an 18 × 18 pixel² sliding window, with 90% overlap in both directions. The wave propagates perpendicular to the direction of particle displacements, which is shown along x and y directions. The shear storage (G′) and loss (G″) moduli were reconstructed using the 2D-NoBC-NLI method (Figure 1c).

2.2 2D-NoBC-NLI

As briefly introduced, the complex shear modulus was reconstructed from the harmonic displacement field utilizing the subzone NLI approach [21], which divides the total region of interest into a set of overlapping subzones. Each subzone is modeled using FE analysis to simulate the viscoelastic response to the excitation. Material properties are iteratively updated in each subzone to minimize the mismatch between the experimentally measured displacement field and the solution obtained from the FE forward model. The mismatch is minimized within the domain of the region of interest, $\mathrm{\Omega }$, through the error function $\mathrm{\varphi }$:

$\mathrm{\varphi }=1/2{\int }_{\mathrm{\Omega }}{\left(\mathbf{u}\left(\mathrm{\theta }\right)-{\mathbf{u}}_{\mathit{m}}\right)}^H\left(\mathbf{u}\left(\mathrm{\theta }\right)-{\mathbf{u}}_{\mathit{m}}\right)d\mathrm{\Omega }(1)$

where the bold font denotes symbols representing tensors and vectors. Here, $\mathbf{u}\left(\mathrm{\theta }\right)$ represents displacements calculated via the FE forward model and the current estimate of material properties, $\mathrm{\theta }$, whereas ${\mathbf{u}}_{\mathit{m}}$ are measured displacements. The superscript H indicates the Hermitian (complex conjugate) transpose. To minimize this function and determine material mechanical properties, the conjugate gradient (CG) method was used, which relies on the gradient of the objective function with respect to the material properties, ${\mathrm{\varphi }}_{\mathrm{\theta }}^{\prime }$, to converge toward an optimal solution for $\mathrm{\theta }$.

For the forward problem used to calculate the displacement field $\mathbf{u}$ as a function of material properties $\mathrm{\theta }$, NLI solves a FE discretization of the time-harmonic Navier’s equation for heterogeneous, isotropic, and viscoelastic materials. The governing equations for elastic wave propagation in such a medium under harmonic excitation, written in the displacement–pressure formulation for nearly incompressible materials [26], are presented in Appendix A1. The weak form of the equilibrium system, detailed in [26] is given by the functional A, where we define the test functions, $\mathrm{W}=\left(\mathbf{w},\mathrm{q}\right)$, consisting of the test displacement field, $\mathbf{w},$ and a test pressure field, q, which, by definition, are set to be zero on $\mathrm{\Gamma }$ :

$\begin{array}{ll}\hfill & \hfill \mathrm{A}\left(\mathrm{W},\mathrm{U}\left(\mathrm{\theta }\right);\mathrm{\theta }\right)={\int }_{\mathrm{\Omega }}G\left(\nabla \mathbf{u}+{\nabla }^T\mathit{u}\right):\nabla \mathbf{w}-p\mathbb{I}:\nabla \mathbf{w}-{\omega }^2\rho \mathbf{u}·\mathbf{w}d\mathrm{\Omega }\\ \hfill & +{\int }_{\mathrm{\Omega }}-\left(\nabla ·\mathbf{u}\right)\mathrm{q}-\frac{pq}{K}d\mathrm{\Omega }\hfill \end{array}(2)$

In Equation 2, $\mathbf{u}$ is the 3D vector displacement field in meters [m], $p$ is a scalar pressure field in Pascals [Pa], and $\mathbb{I}$ is the identity matrix representing the isotropic component of the stress tensor. Material properties include the complex valued shear modulus $G$, in Pascal [Pa], the bulk modulus $K$, in Pascal [Pa], and the density of the material $\rho$, in [kg.m^-3], such that $\mathrm{\theta }=\left\{G,K,\rho \right\}$. In this time-harmonic case, $\omega$ is the angular frequency of the harmonic elastic wave excitation [rad.s^-1]. These equations are applied within the domain Ω, which represents the region of the material under study. The traditional subzone NLI method [21] requires the full three-component measured displacement vector ${\mathbf{u}}_{\mathrm{m}}$ to enforce the boundary conditions and to ensure a well-posed forward problem, as detailed in Appendix A1. However, in many practical settings, only two in-plane displacement components are available. To address this limitation, an alternative formulation has been developed to enable solution of the forward and adjoint problems using only in-plane displacement data.

The aim is to enforce the constraint defined by the weak form of the governing equation in Equation 2, while simultaneously minimizing the data misfit defined in Equation 1, without relying on explicit boundary displacement conditions or the out-of-plane motion measurement. To this end, we define the Lagrangian functional L, combining the objective function describing the displacement error, and constraints of the viscoelastic forward problem given by A, as expressed by:

$\mathrm{L}\left(\mathrm{W},\mathrm{U}\left(\mathrm{\theta }\right);\mathrm{\theta }\right)=\mathrm{\varphi }\left(\mathrm{U}\left(\mathrm{\theta }\right)\right)+\mathrm{A}\left(\mathrm{W},\mathrm{U}\left(\mathrm{\theta }\right);\mathrm{\theta }\right)(3)$

This Lagrangian framework enables solving the inverse problem as a constrained optimization problem, where material properties $\mathrm{\theta }$ are identified to minimize the mismatch between measured and simulated displacement fields, while ensuring that the simulated displacement field $\mathbf{u}\left(\mathrm{\theta }\right)$ satisfies governing equations. To obtain the optimal solution, directional derivatives of the Lagrangian with respect to U, W, and $\mathrm{\theta }$ are computed as follows:

${\mathrm{L}}^{\prime }={\mathrm{L}}_{\mathrm{U}}^{\prime }\left(\delta \mathrm{U}\right)+{\mathrm{L}}_{\mathrm{W}}^{\prime }\left(\delta \mathrm{W}\right)+{\mathrm{L}}_{\mathrm{\theta }}^{\prime }\left(\delta \theta \right)(4)$

Optimal conditions for the system defined by Equation 3 are obtained by equating each term of Equation 4 to zero, resulting in three equations that form the basis of the coupled adjoint-based gradient computation and optimization algorithm. These equations are:

1. ${\mathrm{L}}_U^{\prime }\left(\delta \mathrm{U}\right)=0$

$\begin{array}{ll}\hfill & \hfill {\int }_{\mathrm{\Omega }}\left[{\left(\mathbf{u}\left(\mathrm{\theta }\right)-{\mathbf{u}}_{\mathbf{m}}\right)}^{\mathrm{H}}-\nabla ·\left[\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^{\mathrm{T}}\mathbf{w}\right)\right]-\nabla \mathrm{q}-{\mathrm{\omega }}^2\mathrm{\rho }\mathbf{w}\right]·\delta \mathrm{u}\\ \hfill & -\left[\left(\nabla ·\mathbf{w}\right)+\frac{\mathrm{q}}{\mathrm{K}}\right]\mathrm{\delta }\text{pd}\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\left[\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^{\mathrm{T}}\mathbf{w}\right)·\mathrm{\delta }\mathbf{u}\right]·\text{nd}\mathrm{\Gamma }=0.\hfill \end{array}(5)$

2. ${\mathrm{L}}_{\mathrm{W}}^{\prime }\left(\delta \mathrm{W}\right)=0$

$\begin{array}{ll}\hfill & \hfill {\int }_{\mathrm{\Omega }}\left[-\nabla ·\left[\mathrm{G}\left(\nabla \mathbf{u}+{\nabla }^{\mathrm{T}}\mathbf{u}\right)\right]+\nabla \mathrm{p}-{\mathrm{\omega }}^2\mathrm{\rho }\mathbf{u}\right]·\mathrm{\delta }\mathbf{w}-\left[\left(\nabla .\mathbf{u}\right)+\frac{\mathrm{p}}{\mathrm{K}}\right]\mathrm{\delta }\text{qd}\mathrm{\Omega }\\ \hfill & +{\int }_{\mathrm{\Gamma }}\left[\mathrm{G}\left(\nabla \mathbf{u}+{\nabla }^{\mathrm{T}}\mathbf{u}\right)-\mathrm{p}\mathbb{I}\right]·\mathrm{\delta }\mathbf{w}·\text{nd}\mathrm{\Gamma }=0.\hfill \end{array}(6)$

3. ${\mathrm{L}}_{\theta }^{\prime }\left(\delta \theta \right)=0$

$\begin{array}{ll}\hfill & \hfill {\mathrm{L}}_{\mathrm{\theta }}^{\prime }\left(\delta \theta \right)={\mathrm{L}}^{\prime }={\mathrm{\varphi }}_{\mathrm{\theta }}^{\prime }\left(\mathrm{U}\left(\mathrm{\theta }\right)\right)-{\mathrm{A}}_{\mathrm{\theta }}^{\prime }\left(\mathrm{W},\mathrm{U}\left(\mathrm{\theta }\right);\mathrm{\theta }\right)=0\to {\mathrm{\varphi }}_{\mathrm{\theta }}^{\prime }\left(\mathrm{U}\left(\mathrm{\theta }\right)\right)\\ \hfill & ={\int }_{\mathrm{\Omega }}\left(\nabla \mathbf{u}+{\nabla }^{\mathrm{T}}\mathbf{u}\right):\nabla \mathbf{w}\mathrm{d}\mathrm{\Omega }.\hfill \end{array}(7)$

Here, we see that the gradient of the objective function with respect to the material properties ${\mathrm{\varphi }}_{\mathrm{\theta }}^{\prime }$ is given directly by Equation 7, and can be calculated given the fields $\mathbf{w}$ and $\mathbf{u}$, provided by solutions of Equations 5, 6, respectively. Equation 6 describes the standard FE forward problem for the unknown displacement field, which has been described in detail in numerous sources (see, for example, [25]). Equation 5 describes the so-called adjoint problem, described in detail by Tan et al. [30]; its strong form counterpart is given in Appendix A2. This adjoint problem is well-posed and can be solved for the adjoint field $\mathbf{w}$, given measured and calculated 3D displacement fields. However, the forward problem (Appendix A1) still requires explicit Dirichlet boundary data on u. To circumvent this, we adopted a coupled adjoint field (CAF) formulation in which forward and adjoint systems are solved simultaneously (Kurtz et al. [26]). In this framework, known boundary conditions on w ($\mathbf{w}=0$ $\text{on\hspace{0.5em}}\mathrm{\Gamma }$) substitute for unknown conditions on u in a combined, simultaneous solution for u and w. This combined system can be expressed as shown below, where the three components of the elastic equilibrium equation (Equation 8 a–c) are expanded to illustrate the development that follows:

$\widehat{x}·\left(\left[\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^T\mathbf{w}\right)-q\mathbb{I}\right)+{\omega }^2\rho \mathbf{w}\right]-\left(\mathbf{u}-{\mathbf{u}}_{\mathit{m}}\right)\right)=0in\mathrm{\Omega },\left(\mathrm{a}\right)$

$\widehat{y}·\left(\left[\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^T\mathbf{w}\right)-q\mathbb{I}\right)+{\omega }^2\rho \mathbf{w}\right]-\left(\mathbf{u}-{\mathbf{u}}_{\mathit{m}}\right)\right)=0in\mathrm{\Omega },\left(\mathrm{b}\right)$

$\widehat{z}·\left(\left[\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^T\mathbf{w}\right)-q\mathbb{I}\right)+{\omega }^2\rho \mathbf{w}\right]-\left(\mathbf{u}-{\mathbf{u}}_{\mathit{m}}\right)\right)=0in\mathrm{\Omega },\left(\mathrm{c}\right)$

$\nabla ·\mathbf{w}+\frac{q}{k}=0in\mathrm{\Omega },\left(\mathrm{d}\right)(8)$

$\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{u}+{\nabla }^T\mathbf{u}\right)-p\mathbb{I}\right)+{\omega }^2\rho \mathbf{u}=0in\mathrm{\Omega },\left(\mathrm{e}\right)$

$\nabla ·\mathbf{u}+\frac{p}{k}=0in\mathrm{\Omega },\left(\mathrm{f}\right)$

$\mathbf{w}=0on\mathrm{\Gamma }.\left(\mathrm{g}\right)$

In Equation 8, $\widehat{x},\widehat{y},$ and $\widehat{z}$ denote unit vectors along the Cartesian directions. Equation 8 c indicates that the knowledge of the out-of-plane measured displacements in the z direction, $\widehat{z}·\left({\mathbf{u}}_{\mathbf{m}}\right)$, is required. Here, the inherently 2D nature of the measurement data does not provide this measured motion component, as described in [20]. To circumvent this limitation to 2D data, we introduce a dimension reduction approach to the coupled equilibrium Equation 8, where we replace the condition on out-of-plane motions used for the CAF formulation with an alternative condition based on the divergence free nature of the displacement field for incompressible materials, i.e., $\left(\frac{\mathrm{\partial }{\mathbf{u}}_{\mathrm{z}}}{\mathrm{\partial }\mathrm{z}}+\frac{\mathrm{\partial }{\mathbf{u}}_{\mathrm{x}}}{\mathrm{\partial }\mathrm{x}}+\frac{\mathrm{\partial }{\mathbf{u}}_{\mathrm{y}}}{\mathrm{\partial }\mathrm{y}}\right)=0$. This constraint is fully consistent with the assumption of near-incompressibility introduced earlier, which was used to ensure numerical stability in the formulation. Together, these elements represent complementary strategies for modeling incompressible materials. In the 2D dimension reduction formulation proposed here, we apply this condition to the out-of-plane component of the calculated displacement field, ${\mathbf{u}}_{\mathit{z}}$, to the in-plane components of measured displacements ${\mathbf{u}}_{\mathbf{\text{mx}}}$ and ${\mathbf{u}}_{\mathbf{\text{my}}}$. Thus, Equation 8 is modified to become:

$\widehat{x}·\left(\left[\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^T\mathbf{w}\right)-q\mathbb{I}\right)+{\omega }^2\rho \mathbf{w}\right]-\left(\mathbf{u}-{\mathbf{u}}_{\mathit{m}}\right)\right)=0in\mathrm{\Omega },$

$\widehat{y}·\left(\left[\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^T\mathbf{w}\right)-q\mathbb{I}\right)+{\omega }^2\rho \mathbf{w}\right]-\left(\mathbf{u}-{\mathbf{u}}_{\mathit{m}}\right)\right)=0in\mathrm{\Omega },$

$\widehat{z}·\left(\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^T\mathbf{w}\right)-q\mathbb{I}\right)+{\omega }^2\rho \mathbf{w}\right)-\left(\frac{\partial {\mathbf{u}}_{\mathbf{z}}}{\partial z}+\frac{\partial {\mathit{u}}_{\mathbf{m}\mathbf{x}}}{\partial x}+\frac{\partial {\mathit{u}}_{\mathbf{m}\mathbf{y}}}{\partial y}\right)=0in\mathrm{\Omega },$

$\nabla ·\mathbf{w}+\frac{q}{k}=0in\mathrm{\Omega },(9)$

$\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{u}+{\nabla }^T\mathit{u}\right)-p\mathbb{I}\right)+{\omega }^2\rho \mathbf{u}=0in\mathrm{\Omega },$

$\nabla ·\mathbf{u}+\frac{p}{k}=0in\mathrm{\Omega },$

$\mathbf{w}=0on\mathrm{\Gamma }.$

By solving the novel formulation introduced in Equation 9 for u and w, which enables reconstruction using only 2D in-plane displacement data, the gradient of the objective function can then be calculated directly via Equation 7 and used to minimize the objective function in Equation 1, via non-linear conjugate gradient methods. Note that for the elastography problem described here, derivatives with respect to material properties $\mathrm{\theta }$ are only considered for the two components of the complex shear modulus $G$, which are described by spatially distributed parameter fields, which are optimized to minimize $\mathrm{\varphi }$ and form the resulting elastography image. The density and bulk modulus are considered constant within the material, and set respectively to 1,000 kg/m³ and 2.2 × 10⁹ Pa. Gaussian spatial filtering with a width of 80% of the mesh resolution was applied at the end of the process to further stabilize the reconstruction of $G$.

For the FE method used in the 2D-NoBC-NLI introduced here, we applied approximately 12.5 nodes per wavelength, with the same mesh resolution used for displacement fields (u and w) and material properties, $G$. For the out-of-plane mesh resolution, we used 1/5 of the in-plane mesh resolution. The subzone size was set to be close to one mechanical wavelength, to respect the effective elastography diffraction limit described in [31]. To ensure converged minimization of the objective function, $\mathrm{\varphi }$, we performed 1,000 iterations using the conjugate gradient method.

2.3 Numerical validation (finite element model)

To investigate reconstruction capabilities of the 2D-NoBC-NLI method in the context of OME, synthetic displacement datasets were generated by solving the viscoelastic forward problem (Appendix A1) using FE methods in 3D cell-mimicking geometries. Three computational models were used, as detailed in Table 1.

Table 1

Table 1. Parameters used in numerical models.

The first model was a 75 µm homogeneous sphere, mimicking microbead experiments (Figure 2a). The second model consisted of a 75 µm heterogeneous sphere with a concentric 30 µm spherical inclusion, mimicking a spherical cell with its nucleus (Figure 2b). The third model was developed to mimic an asymmetric cell, as shown in Figure 2c. The base diameter was 75 μm, the height was 35 μm, and the nucleated inclusion was a 30-µm sphere, making the model asymmetric in the z-direction. Transverse displacements in selected 2D planes for these models are shown in Figures 2d–f, respectively. Note that all three of these models were subjected to asymmetric loading, where prescribed displacements were applied in the zone indicated on the >x surface, while fixed, zero displacement conditions were prescribed on the zone indicated on the <x surface.

Figure 2

Figure 2. (a) A homogeneous sphere viscoelastic FE model. (b) A heterogeneous sphere viscoelastic FE model with an inclusion. (c) An asymmetric cell mimic viscoelastic FE model with an inclusion. The y-direction displacements in reconstructed planes are shown in (d) for the homogeneous sphere model, (e) for the heterogeneous sphere model, and (f) for the asymmetric cell mimic model. BC: boundary condition.

Dirichlet boundary conditions were applied to replicate the experimental setup. For the first two spherical models, two circular projections were positioned at opposite ends of the outer sphere. One projection, with a diameter of 20 μm, represented the holding micropipette and was assigned a fixed boundary condition, indicated by the pink color in Figures 2a,b. The other projection, also 20 µm in diameter, mimicked the vibrating micropipette and was assigned a prescribed displacement boundary condition, indicated by the green color in Figures 2a,b, with a harmonic displacement amplitude of 1 µm in the transverse direction y. For the asymmetric cell-mimicking model, a single circular projection at the side, shown in green in Figure 2c, served as the vibrating boundary condition, while the bottom of the model, shown in pink, was fixed to simulate its attachment to a substrate (as in cell culture).

Simulations were carried out using Comsol 5.5, LiveLink for Matlab (Comsol Inc., Sweden) and Matlab R2020b (The MathWorks, United States) software. The material was considered to be linear viscoelastic, isotropic, and nearly incompressible, and material property values for storage and loss moduli were computed using an interpolated power-law equation fitted to experimental multifrequency data at the prescribed actuation frequency. A model-free linear viscoelastic approach was employed, where the complex shear modulus G* = G′ + iG″ was used directly as input, without relying on a predefined rheological model. For the background material, mean experimental values at 40 kHz of G′ = 1,140 Pa and G′′ = 437 Pa were used based on experimental measurements on PAAm microbeads (see Section 2.6). Inclusion moduli were set as 25% lower or higher than the background, corresponding to G′ = 855 or 1,425 Pa and G′′ = 327 or 546 Pa, to represent moderate mechanical contrasts, similar to those typically observed between the nucleus and cytoplasm of biological cells [32–34]. The density and Poisson’s ratio were chosen as 1,000 kg/m³ and 0.499, respectively, to simulate near-incompressibility. All models were meshed using a tetrahedral mesh with a maximum element size of λ/10, where λ is the wavelength of the mechanical elastic wave in the background material, given by $\mathrm{\lambda }=\frac{\sqrt{2\left({G^{\prime }}^2+{G^{″}}^2\right)}}{\rho \left(G^{\prime }+\sqrt{\left({G^{\prime }}^2+{G^{″}}^2\right)}\right)}.\frac{1}{f}$. The number of elements used for the FE calculations are provided in Table 1. The complex 2D displacement data from the FE simulation were interpolated to match the experimental data resolution, as described in Section 2.1, and transformed into the time domain to replicate the format of experimental displacements used in the inverse reconstruction process. Finally, relative errors between reconstructed property results and ground truth values for G were computed for storage moduli using: $\frac{\left|G_{recon}^{\prime }-G_{true}^{\prime }\right|}{G_{true}^{\prime }}\times 100$, and for loss moduli: $\frac{\left|{\mathrm{G}}_{recon}^{″}-G_{true}^{″}\right|}{G_{true}^{″}}\times 100$. These were calculated pixel-wise and then averaged over the entire area of interest.

2.4 Effect of inversion parameters

To assess the sensitivity of the reconstruction process to key parameters that could influence 2D-NoBC-NLI reconstructions, values were systematically varied for the homogeneous sphere model, including the subzone size, out-of-plane mesh resolution, number of FE nodes per wavelength (NPW), and number of shear wavelengths (NSW) within the material. Specifically, the zone size was varied from 50% to 100% of the wavelength, and out-of-plane mesh resolutions were changed from 5% to 100% of the in-plane mesh resolution. The effect of the FE mesh resolution, indicated by the number of FE NPW, on relative errors of storage and loss moduli, as well as the computation time relative to the computation time of 8 NPW, were analyzed. Additionally, the NSW within the material was assessed for its impact on reconstruction accuracy, with the material size to shear wavelength ratio changing from 300% to 100%. This was achieved by increasing the sphere size while keeping constant the actuation frequency, the wavelength-to-subzone ratio, and other inversion parameters. We examined how variations in these parameters impacted the overall reconstruction accuracy. By adjusting these factors, we aimed to gain a deeper understanding of their effects on the precision of 2D-NoBC-NLI reconstructed viscoelastic properties. To assess the reconstruction near optimal inversion parameters for asymmetric and heterogeneous situations, parameters were applied to both heterogeneous sphere and asymmetric cell mimic models with either softer or stiffer inclusions relative to the background. In these cases, the excitation frequency was increased to 60 kHz to ensure enough shear wavelengths within each region, while keeping storage and loss moduli unchanged. This adjustment enabled reconstructions to be performed at close to optimal inversion conditions.

2.5 Noise stability analysis

The robustness of the 2D-NoBC-NLI reconstruction was also evaluated by introducing Gaussian noise into the frequency-domain displacement data of FE models at 40 kHz, chosen as the midpoint of the experimental bandwidth (20–60 kHz), providing a representative operating condition for the study. The analysis was performed with softer inclusions relative to the background. Noise was applied with zero mean and standard deviations ranging from 0%, 5%, and 10% of the mean displacement amplitude, reflecting displacement uncertainty associated with motion estimation techniques. Gaussian white noise was chosen because electronic/thermal noise and displacement tracking errors in elastography are commonly approximated by Gaussian processes, making it a simple and widely accepted model of uncertainty in the displacement field [35].

2.6 Experimental validation (polymer microbeads)

To experimentally validate the 2D-NoBC-NLI method, OME was performed on polyacrylamide (PAAm) microgel beads. The beads were produced following the method described in [36], with minor modifications. Briefly, a polydimethylsiloxane-based flow-focusing microfluidic chip with channel dimensions of 40 (width) × 60 (height) µm² was used to produce PAAm beads with a mean diameter of 75 ± 5 µm. The PAAm pre-gel mixture was prepared using acrylamide (40% w/v, Sigma-Aldrich, A4058, Germany) as the monomer, bis-acrylamide (2% w/v, Sigma-Aldrich, M1533) as the cross-linker, and ammonium persulfate (0.05% w/v, Cytiva, GE17-1,311-01, Germany) as the free radical initiator. The PAAm pre-gel mixture also contained 31 × 10⁹ latex nanoparticles/mL (Sigma-Aldrich, LB6, mean size of 0.6 µm) and 15% (v/v) OptiPrep (Sigma-Aldrich, D1556) to prevent particle sedimentation during droplet formation. The total pre-gel volume was 545 μL, with a total monomer concentration of 7.9% and a cross-linker-to-monomer concentration ratio of 2.6%. The continuous phase consisted of 2% dSURF surfactant as an emulsion stabilizer (Fluigent, France) in Novec 7,500 oil (3M, United States), supplemented with 0.4% (v/v) N,N,N′,N′-tetramethylethylenediamine (TEMED) as a catalyst (Sigma-Aldrich, T9281, CAS 110-18-9). Following in-drop polymerization, beads were washed and resuspended in 1× phosphate-buffered saline (pH 7.4, Gibco, United States). Beads were stored at 4 °C and shipped between laboratories (Erlangen to Montreal) under controlled temperature conditions. Beads were vortexed to resuspend them uniformly, and allowed to equilibrate to the ambient temperature before experiments. The latex nanoparticles embedded within PAAm microbeads allowed optical contrasted images for elastic wave tracking.

Experiments were conducted over a frequency range of 20–60 kHz, corresponding to wavelengths of approximately 70 μm at 20 kHz and 35 μm at 60 kHz in PAAm. This ensured that one to two shear wavelengths were present within each bead, aiding reliable 2D-NoBC-NLI viscoelastic reconstructions. Since the reconstructed shear modulus is frequency dependent, results were analyzed using a power-law relationship, G′ ∝ f^α and G′′ ∝ f^β. The power-law model was chosen because it effectively captures the frequency-dependent shear modulus observed in biological and polymeric materials, which exhibit a broad distribution of relaxation times [8, 37–39]. Coefficients α and β were determined using least-squares regressions. Reported results are based on measurements over 10 microbeads. These experimentally measured viscoelastic parameters were also used as input reference values for the finite-element simulations described in Section 2.3.

3 Results

3.1 Finite element model

3.1.1 Effect of 2D-NoBC-NLI inversion parameters

The homogeneous sphere FE model with a cell-scale diameter exhibited a storage modulus reconstruction error of 1% and a loss modulus error of 9%. The effect of inversion parameters on reconstruction accuracy for this model is summarized in Figure 3, with optimal parameter ranges highlighted in gray. As shown in Figure 3a, increasing the zone size to wavelength ratio resulted in enhanced accuracy for both storage and loss moduli. Mean relative errors for storage and loss moduli were minimized for zone sizes larger than 70% of the shear wavelength, consistent with [31]. Relative errors of storage and loss moduli remained constant when the out-of-plane mesh resolution was between 20% and 100% of the in-plane mesh resolution (Figure 3b). By increasing the number of FE nodes per wavelength in the homogeneous sphere model, relative errors of storage and loss moduli were reduced overall (Figure 3c). However, the computation time increased with the NPW. Figure 3d shows relative errors of storage and loss moduli as a function of the NSW within the homogeneous sphere model. The material size to shear wavelength ratio changed from 3 to 1, and it is observed that if fewer than 2 wavelengths are present in the material, the relative error for the storage modulus increased from 2% to 24%, while the error for the loss modulus increased from 43% to 84%.

Figure 3

Figure 3. Relative errors of storage (blue) and loss (red) moduli as functions of (a) zone size to wavelength ratio, (b) out-of-plane to in-plane mesh resolution, (c) number of finite element nodes per wavelength (NPW), and (d) number of shear wavelengths (NSW) within the homogeneous sphere model.

Based on the parameter sweeps shown in Figure 3, the following 2D-NoBC-NLI settings were selected as optimal and used in all subsequent reconstructions (Table 2).

Table 2

Table 2. Optimal inversion parameters.

Figure 4 illustrates storage and loss moduli in representative slices of the heterogeneous sphere and asymmetric cell mimic models using optimal inversion parameters. Panels (a) and (d) show ground truth values for storage and loss moduli, while (b) and (e) display the corresponding reconstructed storage and loss moduli. Panels (c) and (f) display relative reconstruction error maps for storage and loss moduli, computed as described in Section 2.3. Across the entire region of interest, relative errors for the heterogeneous sphere model were 3.5% for the storage modulus and 10.4% for the loss modulus; for the asymmetric cell-mimic model, they were 5.0% for the storage modulus and 12.9% for the loss modulus.

Figure 4

Figure 4. (a) Ground truth storage and loss moduli for the heterogeneous sphere model, (b) reconstructed storage and loss moduli, and (c) relative reconstruction errors map. (d–f) Same as (a–c) but for the asymmetric cell mimic model. The average error across the entire heterogeneous region was 3.5% (storage) and 10.4% (loss) for the heterogeneous sphere model, and 5.0% (storage) and 12.9% (loss) for the asymmetric cell mimic model.

Figure 5 illustrates storage and loss moduli in representative slices of the heterogeneous sphere and asymmetric cell mimic models with a stiffer inclusion relative to the background. Panels (a) and (d) show ground truth storage and loss moduli, panels (b) and (e) present corresponding reconstructions, and panels (c) and (f) display relative reconstruction error maps. For the heterogeneous sphere model with a stiffer inclusion, average errors across the heterogeneous region were 4.7% (storage) and 20.0% (loss). For the asymmetric cell mimic model, average errors were 9.5% (storage) and 22.6% (loss).

Figure 5

Figure 5. (a) Ground truth storage and loss moduli for the heterogeneous sphere model, (b) reconstructed storage and loss moduli, and (c) relative reconstruction errors map. (d–f) Same as (a–c) but for the asymmetric cell mimic model. Average errors across the entire heterogeneous region were 4.7% (storage) and 20.0% (loss) for the heterogeneous sphere model, and 9.5% (storage) and 22.6% (loss) for the asymmetric cell mimic model.

Detailed relative errors for background and inclusion regions, for both softer and stiffer inclusions, are summarized in Table 3.

Table 3

Table 3. Relative reconstruction errors (%) in storage and loss moduli for heterogeneous sphere and asymmetric cell mimic models.

3.1.2 Noise level analysis

Reconstruction errors were examined for both heterogeneous sphere and asymmetric cell mimic models under different noise conditions. Both models contained a softer inclusion relative to the background. Note that these reconstructions were performed at 40 kHz, which yields a smaller NSW in each region compared to the 60 kHz reconstructions shown in Figure 4. In Figure 6 (heterogeneous sphere model), panel (a) shows y-direction displacement fields at 0%, 5%, and 10% noise; panel (b) presents corresponding storage modulus relative error maps; and panel (c) displays loss modulus relative error maps. Figure 7 shows the same results for the asymmetric cell mimic model.

Figure 6

Figure 6. (a) Displacement in the y-direction for the heterogeneous sphere model with a softer inclusion relative to the background at 0%, 5%, and 10% added Gaussian noise, (b) relative error maps of the storage modulus, and (c) relative error maps of the loss modulus corresponding to these noise levels. Across noise levels, storage modulus relative errors were 6.0% (0% noise), 13.8% (5% noise), and 20.4% (10% noise); loss modulus errors were 23.0%, 21.0%, and 21.6%, respectively.

Figure 7

Figure 7. (a) Displacement in the y-direction for the asymmetric cell-mimic model with a softer inclusion relative to the background at 0%, 5%, and 10% added Gaussian noise, (b) relative error maps of the storage modulus, and (c) relative error maps of the loss modulus corresponding to these noise levels. Across noise levels, storage modulus relative errors were 6.6% (0% noise), 14.0% (5% noise), and 17.6% (10% noise); loss modulus errors were 23.2%, 22.5%, and 21.4%, respectively.

For the heterogeneous sphere model (Figure 6), at 0% noise, background errors were 3.0% (storage) and 21.0% (loss), whereas inclusion errors were 19.0% (storage) and 51.0% (loss). With increasing noise levels (0%, 5%, and 10%), relative errors of the storage modulus in the background rose from 3.0% to 24.0%, whereas in the inclusion they decreased from 22.0% to 3.0%. Errors for the loss modulus remained approximately 17% in the background and declined slightly from 51.0% to 39.0% in the inclusion. Across noise levels, overall storage modulus relative errors were 6.0%, 13.8%, and 20.4%, while overall loss modulus relative errors remained around 22%.

Similarly, for the asymmetric cell-mimic model (Figure 7), at 0% noise background relative errors were 6.0% (storage) and 14.0% (loss), whereas inclusion relative errors were 10.0% (storage) and 51.0% (loss). With increasing noise (0%, 5%, 10%), relative errors of storage modulus in the background increased to 18.6%, whereas in the inclusion they decreased to 7.0%; relative errors of loss modulus in the background remained around 11%, and in the inclusion they declined from 50.0% to 44.0%. Across noise levels, the overall storage modulus relative errors were 6.6%, 14.0%, and 17.6%, whereas overall loss modulus relative errors remained around 22%.

3.2 Microbead experiments

Two representative examples of reconstructed storage and loss moduli from the estimated 2D-displacement field in a PAAm microbead actuated at 20 kHz and 60 kHz are shown in Figure 8. Each panel includes a bright-field microscopy image (a), 2D displacements along y and x directions (b), and reconstructed storage and loss moduli (c) at 20 kHz (top) and 60 kHz (bottom). Mean values were G′ (20 kHz) = 600 ± 3 Pa, G″ (20 kHz) = 407 ± 4 Pa, G′ (60 kHz) = 3451 ± 14 Pa, and G″ (60 kHz) = 1848 ± 21 Pa (Figure 8c).

Figure 8

Figure 8. (a) Optical images of the PAAm microbead in bright field microscopy. (b) Displacements within the microbead in y and x directions. (c) Reconstructed storage and loss moduli of the microbead. Top panel results are at 20 kHz and bottom panel results are at 60 kHz.

Figure 9 presents storage (G′) and loss (G″) moduli for measurements between 20 and 60 kHz. Results were modeled with power-law relationships (G′ ∝ f^α and G′′ ∝ f^β). Power-law coefficients are α = 1.234 and β = 1.125. Mean coefficients of variation (CV) for reconstructed samples were 0.7% for the storage modulus and 0.9% for the loss modulus.

Figure 9

Figure 9. Storage and loss moduli of PAAm microbeads (N = 10) measured with the 2D-NoBC-NLI reconstruction method between 20 and 60 kHz.

4 Discussion

This study proposed a two-dimensional, boundary condition free, nonlinear inversion approach to solve the inverse problem of microelastography in the presence of incomplete, 2D planar data. The technique was evaluated using a finite element-based simulation study and through experiments on PAAm microbeads. Results suggest potential for applications in cellular elastography, and mapping of viscoelastic properties in complex, heterogeneous, and asymmetric morphologies. The 2D-NoBC-NLI method was validated in experimentally realistic scenarios and was robust in the presence of noise. Compared to previous work, this method offers significant advantages in 2D imaging scenarios, enabling the quantitative measurement of mechanical properties that were previously only measurable qualitatively [40]. Unlike previous OME approaches, which assumed either out-of-plane symmetry [20] or plane strain and stress conditions [40], the proposed method can reconstruct quantitative viscoelastic properties without such assumptions.

Optimization of inversion parameters is crucial for specific applications. In this study, the optimization process was demonstrated using the homogeneous sphere model, which served as a representative case. As demonstrated, the effect of the out-of-plane mesh size on reconstruction accuracy indicated an accuracy becoming stable for out-of-plane mesh resolutions between 20% and 100% of the in-plane mesh resolution. A resolution of 20% was chosen in this study due to faster convergence speed. Furthermore, previous findings indicated that having more than eight nodes per wavelength is critical for robust reconstructions [41]. As observed, increasing the NPW from 8 to 14 in the homogeneous sphere model reduced relative errors of both storage and loss moduli. While increasing the NPW improves resolution and reduces the relative error, it comes at the cost of increased computation time compared to using 8 NPW (Figure 3c).

Additionally, the NLI is inherently influenced by the ratio of the zone size to wavelength [42]. It has been demonstrated that the subzone size can affect reconstruction accuracy. The mechanical property characterization is dependent on the portion of the mechanical shear wavelength present within the reconstruction domain, and when less than half of the mechanical shear wavelength is present, accurate reconstruction is only possible in certain conditions [31]. In this study, we investigated the effect of these ratios on the accuracy of the reconstruction in the in silico model. Increasing the zone size resulted in enhanced accuracy in both storage and loss moduli. As observed (Figure 3a), the optimal zone size-to-wavelength ratio of 0.7 minimized relative errors and stabilized the reconstruction process.

Another limitation at the microscale is that boundary effects can influence the propagation of elastic waves and the type of waves present in the material, potentially affecting the accuracy of the reconstruction if a specific wave type is assumed. The NLI forward problem inherently accounts for different wave types in the wave propagation model, including guided waves, shear waves, and surface waves, by directly solving the governing elastic equilibrium equations. This approach ensures that the method is robust in a variety of experimental conditions, as it does not rely on any assumptions about the wave type. Nonetheless, due to the limitation to in-plane data, the 2D-NoBC-NLI method remains sensitive to the number of shear wavelength available within the material or material size-to-wavelength ratio. By maintaining this optimal ratio and varying the number of shear wavelength within the material from 1 to 3, we found that results remained stable when more than two shear wavelengths were present within the entire mimicked cells. For smaller numbers of wavelengths, source-related artifacts in the simplified sphere simulations became more pronounced, particularly within the first wavelength inside the material. For quantitative reconstruction, at least half a wavelength must be present within the material, as previously established for the zone size criterion [31].

Based on these findings (Figure 3; Table 2), the optimal inversion parameters were determined to be 12 NPW. The out-of-plane mesh resolution was defined as 20% of the in-plane mesh resolution. The subzone size was selected to be as close as possible to one mechanical wavelength. By applying these optimal inversion parameters and ensuring close to two shear wavelengths within each region for the heterogeneous and asymmetric cell mimic models, the 2D-NoBC-NLI reconstruction achieved good accuracy.

Under these conditions, the relative error for the storage modulus in the softer-than-background case was reduced to below 10% for both models, in both the background and inclusion. For the loss modulus, the error was reduced to below 20% in both the background and inclusion. In the stiffer-than-background case, the storage modulus error remained below 10% in both regions, while the loss modulus error was below 23%. This demonstrates the method’s capability to achieve accurate reconstruction at the microscale, even under complex material conditions and geometric asymmetry, using only 2D plane displacements.

Numerical simulation results demonstrated the ability of the 2D-NoBC-NLI to reconstruct heterogeneous viscoelasticity maps without prior information. These in silico experiments showed that even in the presence of up to 10% Gaussian noise, storage modulus reconstructions maintained an error below 20%, in both heterogeneous and asymmetric models. The loss modulus had a higher percentage of errors (up to 24%) due to the low sensitivity to this property in regions comprising limited wave content [31]. However, despite the presence of high noise levels, the stability of the loss modulus was still observed, with errors remaining below 20% in the background and stable across different noise levels. This behavior arises because noise partially compensates for the inherent overestimation of the loss modulus near boundaries and in regions with limited wavelength content, leading to an apparently more stable reconstruction under noisy conditions. In our study, we found that storage modulus errors were generally due to underestimation, particularly in the background regions, whereas loss modulus errors tended to be due to overestimation in both regions.

It should be noted that errors in reconstructed loss moduli were generally larger than those in storage moduli. This limitation is well recognized in elastography [24], as the loss modulus estimation relies on the phase component of the displacement field, which is inherently more sensitive to noise and boundary effects, and errors are thus higher than for the storage modulus. The effect is particularly pronounced when the number of nodes is limited relative to available wavelengths, leading to higher biases and reduced symmetry in reconstructions. In our results, this manifested as elevated loss modulus errors, especially near boundaries and wave source regions. Similar observations have been reported in prior elastography studies [41, 43]. Nonetheless, the present method achieved robust reconstructions, and the accuracy of the loss modulus recovery could be further improved by refining or coarsening the mesh resolution depending on available wavelengths [41]. The error for the storage modulus within the inclusion appeared to be influenced by underestimation in the background and boundary regions, leading to artificially lower error values as the background error increased. This effect was particularly pronounced in the inclusion because it contains a limited number of pixels, making it more sensitive to boundary underestimation.

There was also an asymmetry/asymmetrical behavior in reconstruction errors. This asymmetry primarily arises from source-related artifacts near the vibration actuator, and from boundary effects near edges of the domain. In cases where the inclusion was stiffer than the background, wave reflection and refraction at the material interface could further contribute to the asymmetry and increase the error due to the amplified boundary artifacts, which distort the reconstruction accuracy. The maximum local errors were mainly observed near inclusion boundaries and source regions, where wave interactions and strain gradients were highest. Despite these limitations, the proposed 2D-NoBC-NLI framework achieved robust reconstructions across heterogeneous domains, with storage modulus errors remaining low and loss modulus errors within ranges tolerated in elastography studies. For example, Tomita et al. reported that in a homogeneous viscoelastic cubic models, the maximum error between the recovered and true storage moduli reached 22.7% in 2D inversion [44]. Similarly, Zhang et al. observed that in a 2D inversion of an inclusion model with 3% Gaussian noise, the reconstructed inclusion exhibited an underestimation of about 20% for the storage modulus, whereas the loss modulus showed greater instability [24]. These findings confirm that the magnitude of errors observed in our study is consistent with prior reports, further validating the reliability of the proposed approach. Most importantly, the method avoided any out-of-plane assumptions, offering a methodological advance that enabled viscoelastic mapping in microscale materials under realistic experimental conditions. This ability to accurately model complex three-dimensional behavior without relying on simplifications makes the approach particularly suitable for studying heterogeneous and asymmetric materials, such as biological cells, where traditional methods may fall short.

For heterogeneous materials, spatial resolution plays a crucial role in characterization. Shorter wavelengths (higher frequencies) can theoretically improve resolution and provide more accurate mapping of material properties. Improving the resolution of optical images and displacement maps also contributes to better material property resolution. However, achieving these improvements often comes at the cost of higher computational demands and increased attenuation and noise, particularly at higher frequencies and frame rates. In low-SNR conditions, down sampling displacement maps may improve the stability and robustness of the reconstructed material properties [41], though this comes with a reduction in spatial resolution. Maintaining a minimum NPW of 8 is essential for reliable reconstructions. Under optimal SNR conditions, the achievable spatial resolution aligns with the displacement mesh or acquisition resolution and can approach the optical resolution of the system. For quantitative reconstruction of heterogeneous regions in real-world noisy conditions, it is important that at least half a wavelength is present within the region of interest [31]. Meeting this criterion enables reliable estimation of viscoelastic properties. While increasing the number of nodes per wavelength improves reconstruction accuracy, it also increases computation time. In general, successful experimental applications of the 2D-NoBC-NLI method should consider a balance between the spatial resolution, excitation frequency, computational time, and robustness. By optimizing these parameters, it is possible to achieve high-quality viscoelastic property reconstructions, even for complex heterogeneous materials.

Experimental measurements in microbeads offer advantages over traditional rheology methods [6] due to the very short acquisition times (≈0.4 ms) needed to provide quantitative spatial distributions of the complex viscoelastic shear modulus. This rapid data acquisition makes the method suitable for dynamic cellular processes such as cytoskeletal remodeling. For example, each measurement cycle for the 60 kHz excitation took 17 μs, and for all 25 measurements on a single bead, the total acquisition time corresponds to 0.4 ms. These multiple acquisitions provided more robust/reliable displacement data for the inverse problem.

At high frequencies, both storage and loss moduli increased (Figure 9), confirming the frequency-dependent behavior of PAAm microbeads, consistent with previously reported findings for PAAm gels [39, 45]. The discrepancy between the equivalent Young’s modulus in our results and that of a previously published study by our team [36], which measured purely elastic properties using the same microbead fabrication method, is approximately 20% at the lowest tested frequency (20 kHz). These discrepancies may stem not only from differences in the pre-gel mixture and bead size, but also from the experimental measurement method employed (atomic force microscopy or real-time deformability cytometry), and the frequency at which mechanical properties were measured [36]. Indeed, at low frequencies (several Hz), mechanical properties of entangled macromolecular networks are predominant, while at higher frequencies the behavior of individual polymer dynamics becomes prevailing [37, 38]. In this way, OME, which operates at these high frequencies, is able to reveal the behavior of cytoskeletal filaments and probes local properties of the cytoskeleton rather than the collective dynamics of the entire cytoskeletal network [19]. While advantageous to improve the spatial resolution, a high-frequency excitation may limit wave propagation distance due to increased attenuation. The limit frequency at which elastic waves propagate could be determined by the ratio of the loss modulus (viscosity behavior) to the storage modulus (elastic behavior) [46, 47]. A lower ratio leads to a higher frequency limit. PAAm microbeads, due to their dominant elastic behavior, can propagate waves at high frequencies without significant attenuation for the size of particles considered in this study. Microbead experiments performed here revealed a consistent and uniform distribution of mechanical properties, with CV = 0.7% and 0.9% for storage and loss moduli, in accordance with the inherent mechanical homogeneity of samples, as previously described [36]. The homogeneity of 2D-NoBC-NLI results presented here (Figure 8) for PAAm particles highlights the ability of this method to reconstruct viscoelastic properties of microscale materials in true experimental conditions. Notably, the displacement field behavior in the microbeads mirrors that of the homogeneous sphere simulations, with shear wavelengths spanning one to two wavelengths, highlighting the precision of our model and reinforcing the method’s power in capturing complex material behaviors.

Importantly, the proposed method is not limited to microbeads of 75 µm in diameter. It can be applied to any microscale material, provided that the excitation wavelength is comparable to the specimen size. This can be achieved by selecting an optimal excitation frequency that provides sufficient shear wavelengths while maintaining an acceptable signal amplitude. It is important that the amplitude is sufficiently high to capture the wave’s displacement (vibration) inside the material, despite attenuation. Therefore, experimental designs should balance both high and low-frequency applicable values, considering that a higher spatial resolution at high frequency implies more wave attenuation than at lower frequencies. The wavelength to specimen size criterion ensures that the approach remains valid for smaller cells and for irregular morphologies, consistent with our simulations. Ensuring sufficient shear-wave content improves reconstruction stability and accuracy, as boundary effects dominate when fewer than one or two wavelengths are present. This criterion is important across various elastography frameworks to avoid reconstruction bias. When the shear wavelength is comparable to the material size, direct methods, such as wavelength- or shear-speed-based approaches are particularly sensitive to geometric effects, leading to significant reconstruction bias. Model-based inversion methods, like nonlinear inversion, reduce these effects by solving the full elastic wave equations. However, even with advanced methods, it remains crucial to have at least one to two effective wavelengths within the region of interest for reliable reconstructions [48].

As a proof of concept, it is also worth noting that at the lowest excitation frequency of 20 kHz tested in microbead experiments, the shear wavelength approximately corresponded to the diameter of the specimen. Despite this challenging condition, reconstructed results remained consistent with those obtained using another established method [36], demonstrating the robustness of the proposed approach. Based on our experimental results with PAAm microbeads, we found that frequencies in the range of 20–60 kHz were optimal for our application and for reliably reconstructing material properties. Beyond OME, the 2D-NoBC-NLI approach is also applicable to other elastography techniques, such as MR elastography, ultrasound, and optical elastography. The method’s ability to handle complex materials without relying on simplified 2D assumptions makes it versatile for a variety of 2D imaging modalities.

Despite promising results, the 2D-NoBC-NLI approach has some limitations that need further investigation. Source artifacts (i.e., due to the vibrating micropipette) affect the quality of reconstructions, and the relatively higher percentage of errors observed for the loss modulus in this region requires further refinement. This led to an underestimation of the shear modulus. While the viscoelastic model used for reconstruction provides a more realistic portrayal of material properties compared to purely elastic models, it still has limitations, such as potential errors introduced by incompressibility assumptions. These issues need to be further examined, particularly in microscale experiments at high frequencies, including comparison with other reconstruction methods using the same frequency range. Addressing these challenges will be essential to improving the overall accuracy and robustness of the method, particularly in applications involving complex biological cell morphologies.

5 Conclusion

In conclusion, this study investigated high-frequency OME using the new 2D-NoBC-NLI method as an image reconstruction approach, demonstrating its ability to robustly reconstruct viscoelastic material properties in 3D geometries from 2D displacement data without the use of 2D geometry approximations. The feasibility of the technique was validated with data obtained from FE simulations and experimental results on PAAm microbeads. The coupled adjoint-based optimization technique effectively determined homogeneous and heterogeneous viscoelastic properties, and the sensitivity analysis highlighted the importance of optimizing factors such as zone size, out-of-plane mesh size, and wavelengths within the material to improve the reconstruction accuracy. Overall, the 2D-NoBC-NLI technique showed great potential for broad applications in complex and noisy 2D displacement fields, particularly in 2D imaging modalities, including mechanical cell phenotyping and mapping of cellular viscoelastic properties.

Data availability statement

The raw data supporting the conclusions of this article can be made available by the authors, without undue reservation.

Author contributions

SjG: Validation, Visualization, Methodology, Formal Analysis, Investigation, Software, Conceptualization, Writing – review and editing, Writing – original draft. HN: Writing – review and editing, Investigation, Formal Analysis, Data curation. GF: Supervision, Writing – review and editing. BC: Resources, Methodology, Writing – review and editing. RG: Writing – review and editing, Resources. SlG: Investigation, Writing – review and editing, Resources. JG: Supervision, Writing – review and editing, Funding acquisition, Resources. GC: Conceptualization, Supervision, Funding acquisition, Writing – review and editing, Project administration. EV: Methodology, Supervision, Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. The authors acknowledge the support from the Canadian Institutes of Health Research (PJT-190229) and the Natural Sciences and Engineering Research Council of Canada (RGPIN-2019-06660).

Acknowledgements

We would also like to express our sincere gratitude to the late Jochen Guck, who passed away during the course of this work. His generosity and support made a lasting impact on this research, and his legacy will continue to inspire and motivate our work.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Urbanska M, Guck J. Single-cell mechanics: structural determinants and functional relevance. Annu Rev Biophys (2024) 53:367–95. doi:10.1146/annurev-biophys-030822-030629

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Guck J. Some thoughts on the future of cell mechanics. Biophys Rev (2019) 11:667–70. doi:10.1007/s12551-019-00597-0

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Di Carlo D. A mechanical biomarker of cell state in medicine. SLAS Technol (2012) 17:32–42. doi:10.1177/2211068211431630

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Eroles M, Rico F. Advances in mechanical biomarkers. J Mol Recognit (2023) 36:e3022. doi:10.1002/jmr.3022

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Akinyi T, Grasland-Mongrain P, Bhatt M, Catheline S, Cloutier G. Chapter ten - recent advances in imaging of cell elasticity. In: I Pajic-Lijakovic, and EH Barriga, editors. Viscoelasticity collect. Cell migr. Academic Press (2021). p. 257–96.

Google Scholar

6. Hao Y, Cheng S, Tanaka Y, Hosokawa Y, Yalikun Y, Li M. Mechanical properties of single cells: measurement methods and applications. Biotechnol Adv (2020) 45. doi:10.1016/j.biotechadv.2020.107648

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wu P-H, Aroush DR-B, Asnacios A, Chen W-C, Dokukin ME, Doss BL, et al. Comparative study of cell mechanics methods. Nat Methods (2018) 15:491–8. doi:10.1038/s41592-018-0015-1

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Deng L, Trepat X, Butler JP, Millet E, Morgan KG, Weitz DA, et al. Fast and slow dynamics of the cytoskeleton. Nat Mater (2006) 5:636–40. doi:10.1038/nmat1685

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Zhou H, Liu R, Xu Y, Fan J, Liu X, Chen L, et al. Viscoelastic mechanics of living cells. Phys Life Rev (2025) 53:91–116. doi:10.1016/j.plrev.2025.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Mierke CT. Viscoelasticity, like forces, plays a role in mechanotransduction. Front Cell Dev. Biol. (2022) 10:. doi:10.3389/fcell.2022.789841

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Liu Z, Ling SD, Liang K, Chen Y, Niu Y, Sun L, et al. Viscoelasticity of ECM and Cells—Origin, measurement and correlation. Mechanobiol Med (2024) 2:100082. doi:10.1016/j.mbm.2024.100082

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ma D, Ren J, Qian X. Enhanced viscoelasticity as a biomechanical signature in senescent vascular endothelial cells. Cell Rep. Phys. Sci. (2025) 6. doi:10.1016/j.xcrp.2025.102844

CrossRef Full Text | Google Scholar

13. Vahabikashi A, Park CY, Perkumas K, Zhang Z, Deurloo EK, Wu H, et al. Probe sensitivity to cortical versus intracellular cytoskeletal network stiffness. Biophys J (2019) 116:518–29. doi:10.1016/j.bpj.2018.12.021

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Mandal K, Asnacios A, Goud B, Manneville J-B. Mapping intracellular mechanics on micropatterned substrates. Proc Natl Acad Sci (2016) 113:E7159–E7168. doi:10.1073/pnas.1605112113

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Eskandari MA, Fischer J, Veyret N, Marx D, Betz T. Active intracellular mechanics: a key to cellular function and organization. J Cell Sci. (2025) 138:jcs263874. arXiv preprint arXiv:2501.14538. doi:10.1242/jcs.263874

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Leartprapun N, Adie SG. Recent advances in optical elastography and emerging opportunities in the basic sciences and translational medicine. Biomed Opt Express (2023) 14:208–48. doi:10.1364/BOE.468932

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Li H, Flé G, Bhatt M, Qu Z, Ghazavi S, Yazdani L, et al. Viscoelasticity imaging of biological tissues and single cells using shear wave propagation. Front Phys (2021) 9. doi:10.3389/fphy.2021.666192

CrossRef Full Text | Google Scholar

18. Kabakova I, Zhang J, Xiang Y, Caponi S, Bilenca A, Guck J, et al. Brillouin microscopy. Nat Rev Methods Primer (2024) 4:8. doi:10.1038/s43586-023-00286-z

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Grasland-Mongrain P, Zorgani A, Nakagawa S, Bernard S, Paim LG, Fitzharris G, et al. Ultrafast imaging of cell elasticity with optical microelastography. Proc Natl Acad Sci (2018) 115:861–6. doi:10.1073/pnas.1713395115

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Flé G, Van Houten E, Rémillard-Labrosse G, FitzHarris G, Cloutier G. Imaging the subcellular viscoelastic properties of mouse oocytes. Proc Natl Acad Sci (2023) 120:e2213836120. doi:10.1073/pnas.2213836120

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Van Houten EE, Paulsen KD, Miga MI, Kennedy FE, Weaver JB. An overlapping subzone technique for MR-based elastic property reconstruction. Magn Reson Med (1999) 42:779–86. doi:10.1002/(sici)1522-2594(199910)42:4<779::aid-mrm21>3.0.co;2-z

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Doyley MM. Model-based elastography: a survey of approaches to the inverse elasticity problem. Phys Med Biol (2012) 57:R35–R73. doi:10.1088/0031-9155/57/3/R35

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Bishop J, Samani A, Sciarretta J, Plewes DB. Two-dimensional MR elastography with linear inversion reconstruction: methodology and noise analysis. Phys Med Biol (2000) 45:2081–91. doi:10.1088/0031-9155/45/8/302

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhang Y, Oberai AA, Barbone PE, Harari I. Solution of the time-harmonic viscoelastic inverse problem with interior data in two dimensions. Int J Numer Methods Eng (2012) 92:1100–16. doi:10.1002/nme.4372

CrossRef Full Text | Google Scholar

25. Van Houten EE, Miga MI, Weaver JB, Kennedy FE, Paulsen KD. Three-dimensional subzone-based reconstruction algorithm for MR elastography. Magn Reson Med (2001) 45:827–37. doi:10.1002/mrm.1111

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kurtz S, Wattrisse B, Van Houten EEW. Minimizing measurement-induced errors in viscoelastic MR elastography. IEEE Trans Med Imaging (2024) 43:1138–48. doi:10.1109/TMI.2023.3329293

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Seidl DT, Oberai AA, Barbone PE. The coupled adjoint-state equation in forward and inverse linear elasticity: incompressible plane stress. Comput Methods Appl Mech Eng (2019) 357:112588. doi:10.1016/j.cma.2019.112588

CrossRef Full Text | Google Scholar

28. Diaz MI, Aquino W, Bonnet M. A modified error in constitutive equation approach for frequency-domain viscoelasticity imaging using interior data. Comput Methods Appl Mech Eng (2015) 296:129–49. doi:10.1016/j.cma.2015.07.025

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Porée J, Garcia D, Chayer B, Ohayon J, Cloutier G. Noninvasive vascular elastography with plane strain incompressibility assumption using ultrafast coherent compound plane wave imaging. IEEE Trans Med Imaging (2015) 34:2618–31. doi:10.1109/TMI.2015.2450992

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Tan L, McGarry MDJ, Van Houten EEW, Ji M, Solamen L, Weaver JB, et al. Gradient-based optimization for poroelastic and viscoelastic MR elastography. IEEE Trans Med Imaging (2017) 36:236–50. doi:10.1109/TMI.2016.2604568

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Van Houten E, Geymonat G, Krasucki F, Wattrisse B. General guidelines for the performance of viscoelastic property identification in elastography: a Monte-Carlo analysis from a closed-form solution. Int J Numer Methods Biomed Eng (2023) 39. doi:10.1002/cnm.3741

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Lammerding J. Mechanics of the nucleus. Compr Physiol (2011) 1:783–807. doi:10.1002/j.2040-4603.2011.tb00343.x

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Walther BK, Sears AP, Mojiri A, Avazmohammadi R, Gu J, Chumakova OV, et al. Disrupted stiffness ratio alters nuclear mechanosensing. Matter (2023) 6:3608–30. doi:10.1016/j.matt.2023.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ofek G, Natoli RM, Athanasiou KA. In situ mechanical properties of the chondrocyte cytoplasm and nucleus. J Biomech (2009) 42:873–7. doi:10.1016/j.jbiomech.2009.01.024

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Dietrich O, Raya JG, Reeder SB, Reiser MF, Schoenberg SO. Measurement of signal-to-noise ratios in MR images: influence of multichannel coils, parallel imaging, and reconstruction filters. J Magn Reson Imaging (2007) 26:375–85. doi:10.1002/jmri.20969

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Girardo S, Träber N, Wagner K, Cojoc G, Herold C, Goswami R, et al. Standardized microgel beads as elastic cell mechanical probes. J Mater Chem B (2018) 6:6245–61. doi:10.1039/c8tb01421c

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Rigato A, Miyagi A, Scheuring S, Rico F. High-frequency microrheology reveals cytoskeleton dynamics in living cells. Nat Phys (2017) 13:771–5. doi:10.1038/nphys4104

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Krajina BA, Tropini C, Zhu A, DiGiacomo P, Sonnenburg JL, Heilshorn SC, et al. Dynamic light scattering microrheology reveals multiscale viscoelasticity of polymer gels and precious biological materials. ACS Cent Sci (2017) 3:1294–303. doi:10.1021/acscentsci.7b00449

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Abidine Y, Laurent V, Michel R, Duperray A, Palade L, Verdier C. Physical properties of polyacrylamide gels probed by AFM and rheology. EPL Europhys Lett (2015) 109. doi:10.1209/0295-5075/109/38003

CrossRef Full Text | Google Scholar

40. Mei Y, Feng X, Jin Y, Kang R, Wang X, Zhao D, et al. Cell nucleus elastography with the adjoint-based inverse solver. Comput Methods Programs Biomed (2023) 242. doi:10.1016/j.cmpb.2023.107827

PubMed Abstract | CrossRef Full Text | Google Scholar

41. McGarry MDJ, Van Houten EEW, Johnson CL, Georgiadis JG, Sutton BP, Weaver JB, et al. Multiresolution MR elastography using nonlinear inversion. Med Phys (2012) 39:6388–96. doi:10.1118/1.4754649

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Anderson AT, Johnson CL, McGarry M, Paulsen KD, Sutton BP, Van Houten EE, et al. Inversion parameters based on convergence and error metrics for nonlinear inversion MR elastography, presented at the Proc. Intl. Soc. Mag. Reson. Med (2017) p. 1139.

Google Scholar

43. Manduca A, Bayly PJ, Ehman RL, Kolipaka A, Royston TJ, Sack I, et al. MR elastography: principles, guidelines, and terminology. Magn Reson Med (2021) 85:2377–90. doi:10.1002/mrm.28627

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Tomita S, Suzuki H, Kajiwara I, Nakamura G, Jiang Y, Suga M, et al. Numerical simulations of magnetic resonance elastography using finite element analysis with a linear heterogeneous viscoelastic model. J Vis (2018) 21:133–45. doi:10.1007/s12650-017-0436-4

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Krivega ES, Kotova SL, Timashev PS, Efremov YM. Mechanical characterization of soft biomaterials: which time and spatial scale to choose? Soft Matter (2024) 20:5095–104. doi:10.1039/d4sm00530a

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Torres J, Callejas A, Gomez A, Rus G. Optical micro-elastography with magnetic excitation for high frequency rheological characterization of soft media. Ultrasonics (2023) 132:107021. doi:10.1016/j.ultras.2023.107021

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Laloy-Borgna G, Zorgani A, Catheline S. Micro-elastography: toward ultrasonic shear waves in soft solids. Appl Phys Lett (2021) 118:113701. doi:10.1063/5.0039816

CrossRef Full Text | Google Scholar

48. Daire J-L, Sinkus R, Van Beers B, Vilgrain V. Elasticity imaging via MRI: basics, overcoming the waveguide limit, and clinical liver results. Curr. Med. Imaging Rev. (2012) 8:56–63. doi:10.2174/157340512799220544

CrossRef Full Text | Google Scholar

Appendix A

A1 Forward (displacement–pressure) problem

The time harmonic elastic wave propagation in a heterogeneous, isotropic, nearly incompressible viscoelastic medium is formulated in mixed displacement–pressure form. The governing equations in strong form are given as follows:

$\widehat{\mathrm{x}}·\left(\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{u}+{\nabla }^{\mathrm{T}}\mathbf{u}\right)-\mathrm{p}\mathbb{I}\right)+{\mathrm{\omega }}^2\mathrm{\rho }\mathbf{u}\right)=0in\mathrm{\Omega },\left(\mathrm{a}\right)$

$\widehat{\mathrm{y}}·\left(\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{u}+{\nabla }^{\mathrm{T}}\mathbf{u}\right)-\mathrm{p}\mathbb{I}\right)+{\mathrm{\omega }}^2\mathrm{\rho }\mathbf{u}\right)=0in\mathrm{\Omega },\left(\mathrm{b}\right)$

$\widehat{\mathrm{z}}·\left(\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{u}+{\nabla }^{\mathrm{T}}\mathbf{u}\right)-\mathrm{p}\mathbb{I}\right)+{\mathrm{\omega }}^2\mathrm{\rho }\mathbf{u}\right)=0in\mathrm{\Omega },\left(\mathrm{c}\right)$

$\nabla ·\mathbf{u}+\frac{\mathrm{p}}{\mathrm{K}}=0in\mathrm{\Omega },\left(\mathrm{d}\right)(\mathrm{A1})$

$\mathbf{u}={\mathbf{u}}_{\mathbf{m}}on\mathrm{\Gamma }.\left(\mathrm{e}\right)$

In Equation A1, $\mathbf{u}$ is the 3D vector displacement field in meters [m], $\widehat{x}$, $\widehat{y},$ and $\widehat{z}$ are unit vectors along the Cartesian directions, $p$ is a scalar pressure field in Pascals [Pa], and $\mathbb{I}$ is the identity matrix representing the isotropic component of the stress tensor. Material properties include the complex valued shear modulus $G$, in Pascal [Pa], the bulk modulus $K$, in Pascal [Pa], and the density of the material $\rho$, in [kg.m^-3], such that $\mathrm{\theta }=\left\{G,K,\rho \right\}$. In this time-harmonic case, $\omega$ is the angular frequency of the harmonic elastic wave excitation [rad.s^-1]. These equations are applied within the domain Ω, which represents the region of the material under study, typically with Dirichlet boundary conditions $\mathbf{u}={\mathbf{u}}_{\mathbf{m}}$ specified on the boundary Γ.

A2 Adjoint problem

The strong form of the adjoint problem for the fields $\left(\mathbf{w},\mathrm{q}\right)$ is expressed as:

$\widehat{\mathrm{x}}·\left(\left[\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^{\mathrm{T}}\mathbf{w}\right)-\mathrm{q}\mathbb{I}\right)+{\mathrm{\omega }}^2\mathrm{\rho }\mathbf{w}\right]-\left(\mathbf{u}-{\mathbf{u}}_{\mathbf{m}}\right)\right)=0,in\mathrm{\Omega }$

$\widehat{\mathrm{y}}·\left(\left[\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^{\mathrm{T}}\mathbf{w}\right)-\mathrm{q}\mathbb{I}\right)+{\mathrm{\omega }}^2\mathrm{\rho }\mathbf{w}\right]-\left(\mathbf{u}-{\mathbf{u}}_{\mathbf{m}}\right)\right)=0,in\mathrm{\Omega }$

$\widehat{\mathrm{z}}·\left(\left[\nabla ·\left(\mathrm{G}\left(\nabla \mathbf{w}+{\nabla }^{\mathrm{T}}\mathbf{w}\right)-\mathrm{q}\mathbb{I}\right)+{\mathrm{\omega }}^2\mathrm{\rho }\mathbf{w}\right]-\left(\mathbf{u}-{\mathbf{u}}_{\mathbf{m}}\right)\right)=0,in\mathrm{\Omega }$

$\nabla ·\mathbf{w}+\frac{\mathrm{q}}{\mathrm{k}}=0,in\mathrm{\Omega }(\mathrm{A2})$

$\mathbf{w}=0,on\mathrm{\Gamma }$

System Equation A2 is well-posed and can be solved for the adjoint field w, given measured and calculated 3D displacement fields, $\mathbf{u}$ and ${\mathbf{u}}_{\mathbf{m}}$. Here, $\mathbf{w}$ denotes the complex 3D adjoint displacement field and q the adjoint pressure; all other parameters are defined as in Appendix A1.