Multiple variational image assimilation for accessible micro-elastography

The mechanical properties of organic materials are relevant to biological processes at many scales. For example, pathological and physiological changes in tissue can often be detected by monitoring their elastic properties. At another scale, the stiffness of the substrate has a direct influence on cell function and differentiation. However, current elastographic methods do not yet cover all spatial scales and require complex and expensive experimental setups. In this work, we present a method that could be easily implemented in biology laboratories by using any available fluorescent microscope, regardless of its resolution, and a piezoelectric actuator. To this end, we have developed a variational framework that integrates multiple measurements into a PDE-constrained minimization of a customized image motion descriptor. Preliminary results with a linearly elastic PDE model show that this method is able to extract reliable maps of the rigidity modulus under the noisy imaging conditions characteristic of conventional microscopy and within the deformation ranges of commercially available actuators.


Introduction
Variations in the mechanical properties of biological tissue reflect functional changes at the cell scale [1]. This relationship can be exploited to indirectly characterize pathological and physiological changes. For example, at the organ scale, cirrhotic livers become stiffer; at an intermediate scale, tumors are harder than their surrounding tissue; and at the cell scale, the stiffness of myocites is in phase with the cardiac beat [2]. Therefore, an accurate quantification of the elastic parameters within biological materials at several spatial scales is essential to provide good diagnoses.
Precisely, the aim of elastography is to yield stiffness maps of soft tissues by use of biomedical imaging. The most prominent techniques in elastography are ultrasound and magnetic resonance imaging. These methods are focused on extracting tissue features at the macro scale for clinical diagnosis, maxing out at a spatial resolution in the order of the millimeter [3]. However, finer scales are relevant to many biological processes and could also yield more precise diagnostic tools. Respective examples are the influence of the elastic extracellular matrix on cell migration and differentiation, and the amelioration of cancer-recurrence ratios with more precise tumor excisions. Elastography has progressed into the micro scale by turning to optical imaging, including techniques such as optical coherence tomography [4] and Brillouin microscopy [5]. These techniques usually yield strain measurements and have only been recently complemented with indirect stress sensors [6] or shear waves [7] to yield tissue elasticity. At the smallest scale, atomic force microscopy is able to extract measurements of nanometer precision but requires point by point direct mechanical probing with a cantilever. In any case, all these methods 2 1234567890 ''"" require complex and expensive experimental set-ups that depend on precise calibrations and impede widespread adoption, notably for non-clinical research. In this paper, we propose a computational method to extract elasticity measurements in biological materials that only requires (i) a standard microscopy set-up, available in most biology laboratories, and (ii) a piezoelectric actuator [6] or, possibly, a chip with a vacuum chamber [8] (both commercially available). Moreover, the scale of the resulting measurements is directly linked to the resolutive power of the microscope and thus can range across several scales, potentially bridging the gap between the nano and the micro scale [9] with the advent of super-resolution. We note, however, that our system is model-dependent and, even though in theory any type of imaging can be used, in the case of fluorescence microscopy some materials may require a label in the form of either fluorescent proteins or beads.
More concretely, the experimental set-up (Fig. 1a) consists of a piezoelectric (PZT) actuator that moves a glass slide to compress the sample with great precision, either from the top or from the side, generating one or more boundary conditions. If necessary, the sample can be rotated to generate additional measurements. We note that these actuators and set-ups have been used to induce compression at different scales, generating displacements over a wide range of magnitudes within different samples [6], [10]. For example, the PZT ring actuator by Piezomechanik in [6] is able to induce up to 50% strain in step-wise micrometer-scale displacements on different synthetic and biological materials. Besides, more expensive models are able to offer nanometer accuracy. The strain induced within the sample can thus be regulated at almost any scale and is usually extracted by taking pictures before and after the compression is applied. These can be taken in a range of resolutions by using ultrasounds, optical coherence tomography or, as in the present case, conventional fluorescence microscopy.
Our technique consists in extracting the internal motion of the sample observed in this way under the different prescribed boundary displacements, while simultaneously inferring the elasticity parameter of a physical model of the sample's material. Namely, we minimize a (custom-modified) optical flow functional [11] describing the movement of the pixel intensity within the images (Sec. 3), constrained by a system of PDEs [12] modeling the relation between the movement of the observed material and its elastic parameter of interest (Sec. 2). The constrained minimization is then linearized and embedded in an image pyramid scheme, and solved by a combination of the Finite Element Method and gradient descent based on the adjoint method (Sec. 5). We exemplify this method with an isotropic linearly elastic model [13]. This is already a good approximation because the displacements are prescribed by the user and because the tissue integrity cannot be compromised.  The coupled constrained minimization and the integration of several measurements (Sec. 4) make the method able to work with image noise, specially farther from the compression source (where the strain is smaller), attenuate the error propagation of decoupled methods, and offer implicit measurements of the stress. We present results (Sec. 6) on simulations reproducing the real imaging process. That is, on real microscopy images, under real noise conditions, and using real actuator boundary conditions that induce strains as in the literature.

The forward problem
We consider our sample to be a linearly elastic material with domain Ω ⊂ R d , d ∈ {2, 3}, and boundary ∂Ω. Its potential energy E(u) is a function of the displacement vector field u : Ω → R d : where . . stands for the Frobenius inner product. In more detail, E is a function of two symmetric d 2 -dimensional second rank tensors: the strain ε(u), describing the relative displacement under a deformation; and the stress σ(ε(u)), representing the internal forces within the continuous material caused by the deformation. In our case, we work with small strain theory, ε(u) := 1/2 ∇u + ∇u T ; and we choose isotropic conditions for the constitutive equation, σ(u) := λ tr (ε(u))I + 2με(u) where μ > 0, λ > 0 are the Lamé parameters that characterize the material and I the d-dimensional identity matrix. E is also a function of the body and traction forces b, T. We define the functional derivative of E with respect to u in the direction du as usual: We can then use it to derive the elastic equations of equilibrium from the perspective of energy minimization, d u E = 0, by perturbing the displacement u with w = du. In the present case of linear elasticity, it is enough to use the linearity of operators ε, σ and their symmetries to derive: where a(u, By choosing the weighting and trial spaces as u ∈ H 1 g (Ω) and w ∈ H 1 0 (Ω), where g and 0 are the respective boundary conditions, eq. (3) becomes the standard weak formulation of linear elasticity. We note that we will not consider body forces and that the traction stresses will be weighted out by the Dirichlet conditions, i.e. u = g on ∂Ω.
While the framework presented throughout this paper can be used in 3D imaging, this data is still not widely available in biology. Here and in Sec. 3, we adapt it to 2D imaging, which is the condition appropriate to our experiments. Therefore, we consider a plane-stress model, i.e. σ i,z = 0. In this case, the strain is uniform over the thickness of the sample and the planar components of the displacement u x,y satisfy eq. (3) in d = 2. Thus, we use u and x to refer to the planar d = 2 case from here on. In addition, we have where u z is the out-of-plane component. Given a map of the varying rigidity modulus μ(x) within the sample, x ∈ Ω, the forward problem consists in finding u(x) that satisfies eq. (3).

The data term
The problem is to find μ(x) from the material deformation u(x). However, this displacement field is not directly available and has to be retrieved from images. To this end, two images are acquired. One before applying a known g on the sample boundary ∂Ω, and one after. We define I (x, t) as the function representing the evolution of the image intensity from the first image I (x, 0) to the second I (x, 1). Optical flow techniques [14] usually extract the motion o between the images by (i) assuming the intensity is conserved, i.e. dI/dt = o · ∇I + ∂ t I = 0, and (ii) imposing that the motion is regular in some sense.
Nonetheless, the motion o of the intensity within the image does not directly reflect the original motion in the sample described by u and u z . For example, material can go out of the 2D imaging plane. To amend this, we consider the whole imaging process, from the emission of light to the final orthographic projection (Fig. 1b). We begin by assuming that it is the original brightness B(x, z, t) emitted by the fluorophores in the 3D material that is conserved, i.e. dB/dt = 0. When a laser is shone, the microscope collects the 3D brightness of the sample within the section s of the laser sheet and projects it orthographically onto the 2D imaging plane [15]. The resulting intensity is proportional to the integral of the brightness: I(x, t) ∝ s B(x, z, t)dz. From this and the conservation equation for the brightness we derive a new conservation equation for the intensity: dI/dt = v · ∇I + ∂ t I = ∂ z u z I up to first order (whereas we would recover dI/dt = 0 for a 3D-stack acquisition), where v := s Bu/ s B is the brightness-weighted displacement. Notice that the material derivative of the intensity has gained an extra source term that accounts for out-of-plane movement.
Using the plane-stress assumption and eq. (5), we have ∂ z u = 0, v = u, and thus the conservation equation is exact and relates the motion u of the sample directly with the intensity: dI/dt = −λ∇ · u/(2μ). This is an ordinary differential equation and can be integratedà la [16], assuming a constant displacement field: I (x + u, 1) = I (x, 0) exp (−λ∇ · u/(2μ)). The problem of motion estimation is then to find u that minimizes Since this problem is non-convex, we linearize the inner expression by Taylor expansion. However, the resulting functional is only valid for small displacements of the order of one pixel. Therefore, we embed the problem into a multi-resolution scheme [17], creating a Gaussian pyramid of image pairs (see Sec. 5). The problem is then solved progressively, from the coarsest image pairs of the pyramid (where the biggest displacements are only pixel-sized) to the finest, by warping back the second image at each scale with the displacement solution of the previous scale.
The minimization problem is also not well-posed. To help this, the seminal optical flow is regularized by adding a smoothness constraint to the displacement field. In our case, we can regularize the problem in a more physically sound way: where α > 0 is the parameter that weights the regularization term, relative to the data term D, in the sum D + R.

The inverse problem
Our interest is to find a μ that minimizes D(I, u; μ) + R subject to a(u, w; μ) = 0 ∀ w ∈ H 1 0 (Ω). In the incompressible case, i.e. λ → +∞, inverting μ from u is a well-posed problem for certain boundary conditions [18]. However, several complications arise in practical situations. First, biological materials are often compressible. Second, the high levels of noise and the limited sampling in images make the problem very sensitive. Third, in order to conserve the integrity of the sample, only small displacements can be applied as boundary conditions. As a consequence, the displacement field over the sample can be hard to resolve, specially far from the source and in the presence of noise. For these reasons, we introduce a new functional that incorporates n observations of the sample under different boundary conditions g i that yield n different image pairs I i (x, 0), I i (x, 1). With this strategy, we aim at completing any lack of information stemming from the above-mentioned difficulties by accumulating independent observations. We also have to introduce the corresponding combined weak elastic equation: where U = (u 1 , · · · , u n ) ∈ U := × i H 1 g i (Ω) and W = (u 1 , · · · , u n ) ∈ W := × i H 1 0 (Ω). The final problem thus reads: We use a quasi-Newton method to solve the optimization problem (10) and compute the derivatives with the adjoint method. By using the chain rule to calculate the necessary derivatives with respect to μ, we can separate partial and total derivatives. While the partial derivatives ∂ U D, ∂ μ D, ∂ μ R are easy to calculate from eqs. (6) and (7) (either analytically or via finite differences), computing the remaining term is more costly. Indeed, note that finding d μ U requires solving ∀ W ∈ W, where we use ∂ U A(U, ·; ·) · d μ U = A(d μ U, ·; ·) because a(u i , ·; ·) is linear and thus A(U, ·; ·) as well. Eventually, equation (12) would have to be solved for each of the many possible directions of μ, thus being slow. As a faster alternative, we propose to use the adjoint method.
In this case, we have Using the adjoint method, we require only one solution of (13) (here the right-hand side does not change as opposed to eq. (12)) and we can evaluate the gradient (11) in the different directions of μ at the cost of dot products with the changing partial derivative ∂ μ A.

Numerical implementation
The Galerkin formulation of the problems is obtained by introducing the mixed finite spaces U h ⊂ U and W h ⊂ W. The resulting discretization is solved via the Finite Element Method (FEM) [19] with mixed P1 Lagrange elements. In practice, the variational problem A = 0 is assembled into its sparse matrix using the Fenics library [19,20], and the MUMPS library [21] 6 1234567890 ''"" is used to solve the linear system in parallel. In order to calculate the gradient of the functional D + R with respect to μ, the adjoint system was originally solved by taking the continuous PDE approach shown in (13), however, we are currently using the very convenient discrete approach offered by the dolfin-adjoint library [22,23]. The minimization of the functional is performed by gradient descent, using the limited memory BFGS algorithm [24], a quasi-Newton method implemented in the PETSc library [25].
Regarding the image processing, the image pyramid I k i (x, 0), I k i (x, 1), k ∈ {1 . . m} is created by iteratively applying Gaussian filters and downsampling each original pair of images I i (x, 0) and I i (x, 1), i ∈ {1 . . n}; the spatial gradient ∇I k i (x, 1) is also calculated. This and the image warping are done using the standard Python libraries scipy and numpy [26,27]. At each scale of the pyramid a mesh for the FEM is created automatically by using Mshr and CGAL [28].
The final joint algorithm goes as follows. Starting by the coarsest scale and going towards the finest, at each step k do: warp all the i image pairs I k i (x, 1) with the previous displacement solution U k−1 , i.e. I k i (x + u k−1 i , 1), and calculate their gradient ∇I k i (x + u k−1 i , 1). Use these images to build the linearized functional D k . Minimize the functional by performing gradient descent with initial value μ k−1 until the convergence criterion is met and a new μ k is obtained, then do a final solve of the weak equation to obtain U k and move on to the next scale.

Examples and discussion
In order to test the methodology, we use different synthetic examples derived from real data, to which we apply simulated actuations. This allows us to work with fully controlled conditions. In particular, we use original microscopy images I i (x, 0) of polyacrylamide gels, a well-studied linearly elastic material used widely in biology, and we use real values for the boundary conditions induced by a piezoelectric actuator, generating displacements like those observed in the literature. We also reproduce the noisy conditions of conventional fluorescence microscopy.
We start with a given modulus μ(x) map ( Fig. 2(c)) and solve the forward problem A = 0 under different boundary conditions g i , i ∈ {1 . . n}, to obtain the different displacement fields u i (Fig. 2(b), n = 1). Then we use u i to warp images I i (x, 0) into images I i (x, 1) ( Fig. 2(a)). Finally, we add noise at the level expected from a microscope and subsample the data to obtain the final pairs of images.
We remark that the magnitude of the displacements quickly decreases as one gets farther from the side where the boundary condition is applied, partially due to the high compressibility of the gel (λ ≈ 2 for many biological applications). In fact, the added noise and the sampling amount to an almost complete loss of information at the opposite side ( Fig. 2(a), white). In addition, the side of the image where the compression is applied is another source of error because the sample changes size slightly.
In the bottom row of figure 2, the results of the minimization (eq. 10) are presented. Figure  2(d) shows the μ(x) solution extracted at each scale of the multi-resolution approach: as more scales are considered, the map is refined. Figures 2(f) and 2(e) show a comparison between the extracted μ(x), u(x) and the original data in respective figures 2(c) and 2(b). While the magnitude and shape of μ(x) are extracted correctly, there are some mismatches caused by the small values (relative to the noise) of the displacement field far from the boundary where the compression is applied, and also by the warping of this side of the image. The reconstruction of u(x) is very precise but is also affected by the same problems (the field is slightly translated upwards with respect to the original).
To alleviate the artifacts introduced by these conditions, we use our framework to combine several measurements of the same sample (Fig. 3). In this case, we work with a more complicated map of the rigidity modulus ( Fig. 3(a)), consisting of three inclusions of different magnitude that occlude each other. By integrating additional data in the form of three different boundary conditions, μ(x) is better recovered (Fig. 3(b)) because the effects from warping the boundary 7 1234567890 ''"" are reduced and because the combination of the displacement fields samples the domain more extensively. In this way, all three inclusions are resolved, localized and ordered in magnitude.
Given the fact that we are using real image data, we are very encouraged by our preliminary results. Our method is able to tackle the compressibility of the material and the experimental noise. We believe it fares well as compared to other approaches found in the literature, whether only theoretical (where the inversion is done directly on the displacement field and under less noise, e.g. [13,29]) or applied (see [6,30]). Finally, the system will be further validated experimentally by using phantoms with inclusions of known elasticity.