Machine Learning in Model-free Mechanical Property Imaging: Novel Integration of Physics With the Constrained Optimization Process

The Autoprogressive Method (AutoP) is a fundamentally different approach to solving the inverse problem in quasi-static ultrasonic elastography (QUSE). By exploiting the nonlinear adaptability of artificial neural networks and physical constraints imposed through finite element analysis, AutoP is able to build patient specific soft-computational material models from a relatively sparse set of force-displacement measurement data. Physics-guided, data-driven models offer a new path to the discovery of mechanical properties most effective for diagnostic imaging. AutoP was originally applied to modeling mechanical properties of materials in geotechnical and civil engineering applications. The method was later adapted to reconstructing maps of linear-elastic material properties for cancer imaging applications. Previous articles describing AutoP focused on high-level concepts to explain the mechanisms driving the training process. In this review, we focus on AutoP as applied to QUSE to present a more thorough explanation of the ways in which the method fundamentally differs from classic model-based and other machine learning approaches. We build intuition for the method through analogy to conventional optimization methods and explore how maps of stresses and strains are extracted from force-displacement measurements in a model-free way. In addition, we discuss a physics-based regularization term unique to AutoP that illuminates the comparison to typical optimization procedures. The insights gained from our hybrid inverse method will hopefully inspire others to explore combinations of rigorous mathematical techniques and conservation principles with the power of machine learning to solve difficult inverse problems.

The Autoprogressive Method (AutoP) is a fundamentally different approach to solving the inverse problem in quasi-static ultrasonic elastography (QUSE). By exploiting the nonlinear adaptability of artificial neural networks and physical constraints imposed through finite element analysis, AutoP is able to build patient specific soft-computational material models from a relatively sparse set of force-displacement measurement data. Physics-guided, data-driven models offer a new path to the discovery of mechanical properties most effective for diagnostic imaging. AutoP was originally applied to modeling mechanical properties of materials in geotechnical and civil engineering applications. The method was later adapted to reconstructing maps of linear-elastic material properties for cancer imaging applications. Previous articles describing AutoP focused on high-level concepts to explain the mechanisms driving the training process. In this review, we focus on AutoP as applied to QUSE to present a more thorough explanation of the ways in which the method fundamentally differs from classic model-based and other machine learning approaches. We build intuition for the method through analogy to conventional optimization methods and explore how maps of stresses and strains are extracted from force-displacement measurements in a model-free way. In addition, we discuss a physicsbased regularization term unique to AutoP that illuminates the comparison to typical optimization procedures. The insights gained from our hybrid inverse method will hopefully inspire others to explore combinations of rigorous mathematical techniques and conservation principles with the power of machine learning to solve difficult inverse problems.
Keywords: Inverse problem, elasticity, ultrasound, finite element analysis, neural network constitutive model

INTRODUCTION
Soft tissues are complex structures that exhibit non-linear, time-dependent elastic properties. Variations in mechanical properties can provide information pertaining to tissue health, including detecting and diagnosing lesions of the breast [6,34], liver [26,47,56], and thyroid [9,10,54], monitoring thermal lesions during ablation therapy [41,53], and characterizing atherosclerotic plaques [2,3,11]. Conventional medical imaging modalities are insensitive to this mechanical contrast. Elasticity imaging methods have been developed to fill this niche and directly evaluate the mechanical properties of soft tissue.
Elasticity imaging can be largely grouped into quasi-static or dynamic methods based on the type of mechanical stimulus applied to induce tissue motion. Comprehensive reviews of the various techniques can be found in several articles [33,35,38,39,46,52]. We focus our discussion on quasi-static ultrasound elastography (QUSE), wherein local tissue displacements are estimated from RF echo frames acquired as an ultrasound probe is slowly pressed into the tissue surface. Quasi-static loading allows time for the force stimulus to propagate through the entire tissue section. The resulting forcedisplacement measurements at points provide an enormous amount of information pertaining to tissue material properties throughout the contiguous volume and at boundaries.
Estimating material properties from force-displacement measurements constitutes the inverse problem in QUSE. The majority of quantitative QUSE approaches are built upon model-based optimization methods. While effective, modelbased methods are prone to modeling errors that can lead to artifacts in reconstructed images of pre-selected material property distributions. More importantly, mathematically defined constitutive models are inherently limited in their ability to capture material behavior. This limitation could preclude discovery of clinically relevant soft tissue material properties.
Data-driven inverse methods circumvent modeling errors by removing the constitutive model assumption and building a softcomputational model of material behavior from measurement data. The Autoprogressive Method (AutoP) is one such method. Unlike other data-driven approaches, AutoP combines physical modeling through finite element analysis (FEA) with artificial neural networks (ANNs) to extract the stress-strain relationship embedded within force-displacement measurements. In the context of elasticity imaging, ANNs characterizing the stress-strain behavior of soft tissues can be interrogated to infer material parameters that best summarize the learned mechanical properties.
Many of the research articles that employed AutoP have described the training process, but none have directly compared and contrasted its operation with the more typical model-based approaches. In this review, we aim to more thoroughly examine the operating principles behind AutoP to better highlight the fundamental differences when compared to model-based inverse methods. Of particular importance is the manner in which physical principles are directly exploited in AutoP to generate increasingly accurate estimates of stress and strain distributions from forcedisplacement measurements under arbitrary geometry and loading. We believe the working principles of AutoP can be applied as a data-driven, physics-guided alternative to other difficult boundary value inverse problems. Section 2 reviews model-based inverse methods in QUSE to establish a basis for comparison. We focus on iterative methods utilizing FEA as the forward solver to keep the discussion concise and comparable to AutoP. Section 3 reviews AutoP and builds intuition of its operation by drawing parallels to typical optimization methods. We also cover a novel neural network architecture developed for QUSE and a regularization term uniquely suited for AutoP. Challenges with non-linear viscoelasticity imaging and potential approaches with alternative ANN structures are briefly discussed in Section 4. Data acquisition and displacement estimation details are omitted from our discussions, but are available in [24,25]. Even though such details are important in practice and will affect the efficacy of QUSE methods, we wish to compare the operation of AutoP and model-based methods at a more conceptual level.

MODEL-BASED INVERSE METHODS
Reconstructing a map of material properties from force-displacement measurements is an ill-posed inverse problem. 1 By combining a prior assumption of the underlying mechanical behavior with a system model incorporating physical principles, model-based inverse methods reduce the solution space and estimate parameter values of a pre-selected constitutive model. The goal of QUSE methods is to find the set of spatially distributed material parameters θ(x) that minimize the difference between displacements U(x) computed in a forward problem and measured displacementsŨ(x), where Ω denotes domain of the scan plane(s). The dependence of variables on position x is assumed and not explicitly stated to simplify notation. The cost function C in Eq. 1 is often defined as the L 2 -norm of the difference between computed and measured displacements, augmented by a regularization term R(θ) acting on the set of parameters, where α is a hyperparameter that controls how strongly the regularization term affects the solution.

Finite Element Analysis
In modern methods, the forward problem is usually solved through finite element analysis (FEA). A brief description of the FEA formulation for solid mechanics follows; more thorough treatments can be found in [12,18]. The strong form of the governing equations for a solid continuum under quasi-static loading is where σ denotes the stress tensor and u are displacements throughout the continuum Ω. Displacement boundary conditions (BCs) are represented byũ whereas surface traction BCs are denoted byp. In the latter two equations, Γ D represents the set of boundary points over which displacements are measured, Γ N is the set of boundary points corresponding to measurements of surface tractions with surface normal n, and Γ D ∩Γ N ∅.
The finite element formulation expresses the boundary-value problem as a system of equations, where N e is the number of elements in the mesh, B e is the straindisplacement matrix for a given element, and the summation symbols represent the assembly operation. The matrix K(θ) is determined directly from the material property matrix D(θ). Any displacement BCs are imposed in global displacement vector U. Likewise, force BCs are included in the global force vector P. 2 We chose to express P in Eq. 8 as a function of stress, instead of as a function of the surface traction BCs, to emphasize the relationship between stress and force. Generally, the solution for a mechanical system is found through FEA by discretizing the domain (i.e., "meshing"), applying load BCs in the global node force and displacement vectors, and solving for the remaining "free" nodal displacements and reaction forces.
FEA solutions of a mechanical system must satisfy both equilibrium and compatibility requirements. The former is specified by Eq. 8 and relates stresses to forces. Compatibility, on the other hand, relates displacements to strains and, in simple terms, precludes elements in the mesh from overlapping or separating during deformation. These physical principles combined with object geometry and a valid constitutive model determine the set of admissible deformations under loading ( Figure 1A). The importance of these principles in solving the inverse problem will be made clear in the following sections.

QUSE as an Optimization Problem
Gradient or Hessian-based optimization techniques are typically applied to solve Eq. 1, wherein the parameter estimates are iteratively updated through the gradient of the error computed in Eq. 2. The solution process for QUSE can be generalized to five steps: The diagram in Figure 1A provides a simplified illustration of the role of FEA within an optimization-based solution method while simultaneously highlighting the fact that the cost function  [1] Apply displacement BCs at the mesh boundary [2] Compute the displacements over the whole scan plane through FEA [3] Compare the computed displacements to the measured values (i.e., evaluate Eq. 2) [4] Update the parameter estimates based on the error [5] If solution converges, exit, otherwise return to [1]. directly affects estimated parameter values. Optimization methods often have intuitive geometric interpretations that aid in understanding the solution procedure. For example, consider the method proposed by Kallel and Bertrand [29]. Define the operator T (E) U as the solution of a FEA where the material is assumed to be linear-elastic with spatial distribution of Young's modulus E(x). The goal is to find E(x) that minimizes the cost function where || · || 2 is the L 2 -norm and R + represents the set of positivevalued real numbers.
Adopting an iterative, Hessian-based solution procedure with k indicating the current iteration, the Young's modulus values can be found, Geometrically, the curvature of the forward operator (the Hessian H), guides the Young's modulus update in Eq. 12. A perhaps more intuitive geometric interpretation is achieved by setting the Hessian equal to the negative identity matrix H −I, changing the update Eq. 11 to In this form, the operator gradient forming the sensitivity matrix J k coupled with the displacement error directs the iterative parameter updates. Iterative QUSE methods typically follow a Hessian or gradient-based optimization procedure. A review of relevant methods can be found in [13,36] Algorithm 1. We will show in Section 3 how displacement errors also influence the material properties learned by ANNs trained in AutoP, albeit indirectly.
QUSE methods utilizing only displacement measurements can provide quantitative parameter estimates up to a multiplicative constant, unless the parameter value(s) is known for some portion of the boundary or interior [4,5]. That is, through compatibility requirements embedded within FEA, at best the distribution of relative parameter values can be inferred from displacement BCs. Incorporating measurements of surface force or boundary stresses can provide the additional information necessary to estimate parameter values exactly [51]. This can be seen directly from Eqs. 5, 8: if no values of reaction forces corresponding to displacement BCs are known, there exists no means to calibrate the magnitude of forces in P. In short, both equilibrium and compatibility requirements affect the solution in model-based QUSE but are not directly exploited.
Force and displacement measurements coupled with FEA are insufficient for overcoming the ill-posed nature of the inverse problem, hence the inclusion of a regularization term in Eq. 2. The choice of regularization depends on some presumed property of the solution; e.g., Tikhonov regularization is founded in energy arguments, L 1 -norm to promote sparse solutions, or total variation to encourage smooth solutions [43]. Although effective and often necessary, the form of regularization will significantly alter the solution in model-based methods and requires skill in choosing appropriately. We will introduce a regularization term in Section 3.1 that is unique to AutoP and arises from physical principles, removing operator judgment from its use.
Model-based QUSE methods are effective when the applied load and true tissue material properties match the selected constitutive model and loading geometry. Choosing the wrong constitutive model can cause artifacts in the resulting elasticity images that corrupts or discards clinically relevant diagnostic information not described by the selected model parameters. The limited flexibility of model-based methods makes them ineffective at discovering material properties. Overcoming limitations of model-based methods and inferring mechanical behavior from measurement data without constraints imposed by constitutive model assumptions may allow investigators to find new clinically relevant material properties.

THE AUTOPROGRESSIVE METHOD
Contrary to model-based inverse methods, the goal of AutoP is to describe the mechanical behavior of a material nonparametrically. By replacing the classic mathematically defined constitutive model with an artificial neural network, any arbitrary stress-strain relationship can be captured. This formulation replaces parameters in Eq. 1 with the connection weights of one or more ANNs 3 . The resulting, trained ANN describing the stress-strain behavior is referred to as a neural network constitutive model (NNCM). One example of a NNCM architecture is shown in Figure 2A (left). The NNCM displayed is a feed-forward, fully connected ANN that accepts a strain vector as input and returns a stress vector at the output. Alternative network architectures have been proposed by other investigators, but this form is the basis for most AutoP applications thus far. Even considering the various NNCM architectures utilized within AutoP, to the best of our knowledge all have comprised at most two hidden layers with up to a few dozen nodes each. Compared to deep-learning methods, NNCMs are very small ANNs.
Existing efforts to build model-free descriptions of material behavior typically rely on experimental test data acquired from a sample in a well-defined loading scenario, such as strip extensiometry, to obtain stress-strain data from forcedisplacement measurements. Recent investigations in datadriven mechanics use this type of measurement data directly (i.e., no NNCM) to solve computational mechanics problems [31,3 Node biases would also comprise the set of parameters if included in the network architecture.
Frontiers in Physics | www.frontiersin.org July 2021 | Volume 9 | Article 600718 32, 37, 50]. Alternatively, stress-strain measurement data can be used to train NNCMs for later use in computational problems (see reviews [14,15,44,55]). More recently, researchers have developed physics-informed deep neural networks (PINNs) that embed physical laws into the cost function used in backpropagation-based training [40]. Some measure of success has been achieved with each of these methods; however, the resulting soft-computational constitutive model, in whichever applicable form, is fully defined by the measurement data acquired from a homogeneous sample. That is, the models are trained once and deployed. When used for inference against input data not within the domain of the training set, the solution becomes unreliable, analogous to modeling errors in model-based inverse methods. Furthermore, these methods rely on the existence of stress-strain training data. Unfortunately, the complexity of soft tissues precludes excision of geometrically precise, homogeneous samples for measurement in well-defined loading scenarios. Nor do ex vivo load tests provide a good measure of soft tissue properties in vivo. In short, current data-driven methods fail to extract relevant information from the available force-displacement measurements acquired in a normal clinical setting.
The Autoprogressive Method [17] circumvents limitations with model-based and other data-driven methods through a tight coupling of NNCMs with FEA and knowledge of both internal and external object geometry. It has been successfully used to build soft-computational mechanical models of many different materials in a multitude of loading scenarios, including soils [19,20], concrete [27,28], steel [30,57], red blood cells [48,49], and the cornea [16]. AutoP has even been applied to thermal constitutive models [1]. More recently, we have been developing an adaptation of AutoP for QUSE using both boundary and interior displacement measurements [22]. The caveat that internal geometry must be known precludes most clinical imaging applications and is addressed by introducing a new type of NNCM, discussed in Section 3.1.
AutoP is a method for generating information from measurements. When applied to inverse problems in continuum mechanics, the goal is to recover the stress and strain fields induced by measured forces and displacements in a model-free way. Replacing the mathematical constitutive model with an ANN is paramount to provide the flexibility inherent to said networks, but alone is not a sufficient strategy. The two key components that make AutoP work are 1) the application of force and displacement BCs in separate FEAs-FEA σ and FEA ε -to obtain physically consistent estimates of stresses and strains through known physical principles and 2) NNCMs participate in the generation of their own training data by acting as the material model in the assembly of the system of equations in Eq. 6. The former point is analogous to model-based inverse methods wherein the parameter estimates in the current iteration affect the displacements computed in the forward problem. For clarity, the FEA equations can be rewritten as We introduce the operator NN(·) to represent inference by the NNCM. The analytical expression for D NN is derived in [21].
Imposing measured forces and displacements in separate FEAs exploits stress equilibrium and compatibility requirements. In FEA σ , measured forces are applied as force BCs (σ · n p) 4 . Stress equilibrium, defined in Eq. 3, relates forces to stresses and we therefore claim the stresses σ σ computed in FEA σ are physically consistent approximations of the true stress field. In a separate FEA ε , corresponding displacement measurements are applied as displacement BCs (u ũ). Because compatibility requirements relate displacements to strain, we claim the strains ε ε computed in FEA ε are physically consistent approximations of the true strain field. The AutoP training process is summarized in Figure 1B and can be broken into the five steps listed in Algorithm 2.
Linear-elastic pretraining is used only to provide NNCMs with a physically consistent starting state. Unlike model-based reconstruction methods that can be sensitive to parameter initialization [4], investigations with AutoP have demonstrated that linear-elastic pretraining does not constrain the mechanical properties learned by the NNCMs.
Convergence in AutoP is determined by comparing the values of a cost function to predefined limits. In most AutoP applications thus far, the cost function has assumed the form [17].
Aside from the use of the L 1 -norm to diminish the effects of outliers, it is important to note that the magnitude of the cost function has no direct effect on the backpropagation-based training of the NNCMs. Even so, training the NNCMs with the stresses and strains computed in FEA σ and FEA ε , respectively, adjust the connection weights so that the material properties described by the NNCMs more accurately resemble the true material properties. We also note that the because the cost function does not directly influence NNCM training, adjustments to Φ(θ) do not significantly affect AutoP whereas altering the cost function in model-based inverse methods can significantly change the resulting parameter estimates.
It is not immediately clear why the combination of the stresses and strains from FEAs solved under separate, but conjugate, boundary conditions should lead the ANNs to the correct material properties. We will demonstrate that it is in fact an implicit displacement error that guides the learning process. Define u σ to be the full set of displacements computed in the solution of FEA σ , including any imposed BCs. Likewise, u ε represents the displacements computed in the solution of FEA ε . By definition of a NNCM, σ NN NN(ε) NN[f u (u)], where f u (·) is a function that receives a displacement vector as input and returns the corresponding strain vector. Let Δu u σ − u ε .
Thus, u σ u ε + Δu and therefore σ σ NN[f u (u ε + Δu)]. What does this mean? Inference by the NNCM using physically consistent strains ε ε as input results in physically consistent stresses σ σ only when the input is augmented by the displacement error.
To better understand this interaction, assume the NNCM describes a linear-elastic material so that σ NN ≈ D NN · ε and deformation is infinitesimal, leading to f u (u + Δu) f u (u) + f u (Δu). We also denote the strains computed in FEA σ as ε σ while σ ε NN(ε ε ) are the stresses output by the NNCM when strains computed in FEA ε are supplied as input to the network. Training a NNCM with a backpropagation-based method consists of minimizing a cost function, in this case where R θm denotes the NNCM connection weights and the overline is used to indicate a constant term. NN is the same NNCM used to solve FEA σ and FEA ε whereas NN(ε ε ; θ) changes each training iteration as the NNCM weights are updated. Applying the aforementioned linear approximations, Eq. 22 can be expressed as If we consider the first iteration of training when D NN D NN (R θ m ), the cost function further reduces to ≈ argmin and reverts back to Eq. 24 after the first training iteration and the NNCM weights R θ m change.
The preceding argument does not imply NNCMs trained in AutoP are limited to learning linear-elastic material properties. Rather, the case of linear materials was chosen as an example to best highlight the effect of displacement error in NNCM training. We emphasize that, similar to model-based methods, the displacement error is a driving force in the solution of the inverse problem as approached through AutoP. However, in the case of AutoP, the effect of displacement error emerges from the inconsistency between stresses and strains computed in FEA σ and FEA ε . It is a characteristic that arises out of the equilibrium and compatibility requirements rather than an explicitly defined error measure to direct the parameter search.
Much like optimization-based methods, the parameter estimates in the current iteration-the weights of the ANNs-affect the solution of the forward problem. Moreover, the error when comparing the forward problem solution to experimental measurements affects the parameter updates. What sets AutoP apart from typical model-based and machine learning techniques is its self-learning paradigm guided by physical principles. When applied to QUSE, the result is a patient-specific soft-computational material model. Figure 2B aims to illustrate the self-learning property of AutoP. Imposing measurement data as BCs in separate FEAs forces the material Algorithm 2 | AutoP Training of NNCMs.
Frontiers in Physics | www.frontiersin.org July 2021 | Volume 9 | Article 600718 properties learned by the ANNs to satisfy physical principles and be self-consistent.

Cartesian NNCMs
The discussion of AutoP thus far has assumed that the finite element mesh conformed to both internal and external object geometry with a different NNCM assigned for each unique material. Eq. 18 implies this fact in that D NN e (θ) is constant over the entire domain of an element. The FEA solution procedure requires the appropriate NNCM to be selected for the element when evaluating the integral. This formulation of AutoP is not suitable for biomedical imaging applications because the internal tissue structure is normally not known a priori nor segmented into discrete regions with homogeneous material properties. As a consequence, the goal of AutoP must be adjusted to infer both material properties and geometry from the set of force-displacement measurements.
Two significant changes were made toward this goal. First, the structure of the NNCM was changed by adding a second network that interacted with the existing NNCM and was responsible for learning spatial information. We deemed these networks Cartesian neural network constitutive models (CaNNCMs) [23][24][25]. The CaNNCM architecture is illustrated in Figure 2A. To identify the individual subnetworks, we refer to the one which learns the "average" stress-strain relationship of the entire object as the material property network (MPN, Figure 2A left) and the network that maps a Cartesian coordinate input to a spatially varying value as the spatial network (SN, Figure 2A right). A thorough description of CaNNCM theory of operation can be found in [23,24]. Unlike NNCMs previously utilized with AutoP, a single CaNNCM is capable of characterizing the heterogeneous material properties of an entire object. The FEA equations are adjusted slightly to accommodate this change: Spatial information is encoded within a set of spatial scaling values S ε x . The SN is responsible for learning the mapping SN : x → S ε x . Given a spatial location and strain vector, the SN provides the corresponding S ε x and the stress vector is then calculated as σ NN (ε, x) NN(ε/S ε x ), where the division is element-wize and the operator NN(·) refers to the MPN. Determining the spatial values is not trivial, although they can be readily computed from the stresses and strains already available in AutoP in conjunction with the MPN, as detailed in [23]. AutoP training must also be adjusted to accommodate CaNNCMs as described in Algorithm 3.
CaNNCMs work well in QUSE because information pertaining to the interior geometry can be gleaned from internal displacement measurements available through US imaging. In unpublished work, we found that CaNNCMs trained in AutoP using only boundary information were unable to recover internal structure away from the surface. Part of the shortcoming is due to the lack of surface force distribution measurements. Rather, the force measured experimentally is the sum of forces over the face of the US transducer that provides information about the mean stress distribution over the surface. Coupling these boundary data with prior knowledge of internal geometry is sufficient for a set of NNCMs to learn the heterogeneous material properties of an object. Therefore, in order for CaNNCMs to learn both material properties and geometry, some measurement data containing spatial information must be provided.
Interestingly, we found that incorporating internal displacements in FEA ε was insufficient for accurately estimating material property distributions, particularly when measurements were acquired by loading from a single direction [25]. The reason can be seen from Eq. 27. By making the stress term explicitly dependent on spatial position, additional force measurements populating P are required to ensure the internal resisting forces computed within the integral are correctly calibrated. Acquiring measurements from multiple sides and loading angles provides an immense amount of information for CaNNCMs at the cost of significantly increased FE modeling complexity and computational load.
An alternative strategy is to add a regularization term to augment the available information. Our solution was to add additional physical constraints during the calculation of S ε x in Step 4a of AutoP. To understand the functionality of the regularization term, we first summarize the detailed methods in [23]. The spatial values are computed by minimizing a cost function very similar to Eq. 21, with the primary difference being the parameters over which the function is minimized changes from the MPN connection weights to the spatial values, The purpose of the spatial values is to further reduce the error in Eq. 21 by storing the information not captured by the spatially independent MPN weights.
Implementing Eq. 28 in AutoP as defined will lead to erroneous material property distribution estimates. Perhaps the most intuitive explanation starts by noting the CaNNCM is pretrained as a homogeneous, linear-elastic material. Then, when computing spatial values in Step 4a, the S ε x found as the solution of Eq. 28 map a heterogeneous strain field to a homogeneous stress field because the surface force BCs applied in FEA σ contain no spatial information. [0] Pretrain CaNNCM as homogeneous, linear-elastic material [1] Apply force BCs in FEA σ , compute σ σ [2] Apply displacement BCs in FEA ε , compute ε ε [3]  We introduced σ-matching regularization to inject spatial information into Eq. 28 by imposing the physical constraint that conjugate force-displacement BCs imposed in separate FEAs should produce the same stress fields, Note that the stress terms are constant, hence the multiplication by sign(σ σ ) to account for the lack of gradient information. We again chose the L 1 -norm to reduce the effects of outliers. If the same linear approximations from Eq. 21 are applied to the σ-matching term, we can show that the displacement error is influencing the result: Adding the regularization term to the cost function for computing S ε x with free parameter β σ to control its influence on the solution, Eq. 28 becomes A significant difference between regularization terms typically used in model-based inverse methods and 33) is the appearance of σ-matching within the L 2 -norm. Because it is constant, simply adding the σ-matching term to the existing cost function would have no effect on the computed spatial values. Adding σ-matching as we have in Eq. 33 can be interpreted as adjusting the (constant) stresses σ σ computed in FEA σ with spatial information encoded within displacement measurements through the displacement error.
We demonstrated in [25] the significant improvement in material property estimates when σ-matching was incorporated into AutoP training, even with a relatively sparse set of forcedisplacement measurements. Contrary to typical forms of regularization, σ-matching is based on a physical principle that will always be true, not a prior assumption about the property the solution should exhibit. Furthermore, it is uniquely suited to AutoP because of the way equilibrium and compatibility requirements are exploited by FEA σ and FEA ε . Other constraints on physical properties will likely be incorporated into AutoP as CaNNCMs are further developed for imaging the non-linear and timedependent properties of soft tissues.

NON-LINEAR ELASTICITY IMAGING WITH AUTOP
The preceding discussion of NNCMS and AutoP was limited to objects under loads inducing small deformation. Mathematically, the strain tensor in the limit of infinitesimal deformation is expressed as ε 1 2 (∇u + ∇u T ). Under this linear kinematic assumption, appropriate when induced strains are less than 5%, elastic solids can be adequately described by a linearelastic constitutive model.
Probing the non-linear mechanical properties of deep structures in QUSE requires application of large loads to sufficiently deform distal tissues. In these scenarios, a large deformation framework must be invoked wherein the strain tensor includes geometrically non-linear terms. As an example, one potential definition of geometrically non-linear strain is ε 1 2 (∇u + ∇u T + ∇u T ∇u). Application of large compressive loads coupled with the non-linear and time-dependent mechanical behavior of soft tissues implies that one or two linear-elastic parameters are ineffective for clinical QUSE. More importantly, the complex non-linear and viscoelastic properties of soft tissues may provide a significant amount of diagnostic information [42], reinforcing the value of reconstructing elasticity images beyond shear modulus.
The problem then arises of selecting an appropriate nonlinear and/or viscoelastic constitutive model that describes tissue behavior and maximizes the diagnostic information in QUSE images. Complexity is further increased when considering the likelihood that different tissues or tissue states (i.e., healthy tissue, benign lesions, and malignant tumors) are best described by different material models. Current model-based approaches assume homogeneity of the underlying constitutive model for the entire imaged tissue section and estimate the spatial variation of the corresponding model parameters. There exists no modelbased QUSE inverse method capable of discovering the set of constitutive models-and estimating values of the corresponding parameters-that adequately describes a heterogeneous group of tissues. The combination of AutoP with CaNNCMs can fill this void.
AutoP has been applied to build soft-computational models of complex materials undergoing small and large deformation when the total geometry was known (e.g. [19,27,28,30,45,58]). To achieve non-linear, time-dependent material property imaging with AutoP, CaNNCMs must be modified to 1) accommodate geometrically non-linear deformation to facilitate non-linear elasticity imaging and 2) incorporate additional stress-strain state information as network input to properly capture a heterogeneous distribution of complex mechanical behaviors. Prior successes with AutoP lead us to believe there exists an adaptation to the CaNNCM architecture and/or AutoP that will accomplish both tasks.
The difficulty of this problem can be understood by comparing the information required for a NNCM to compute a stress vector for materials of varying complexity. For linear-elastic materials, the NNCM output is defined as σ NN NN(ε). Non-linear, viscoelastic materials require time information to fully define the stress-strain response, σ NN NN(ε, t). It is further complicated for objects under large deformation because of the need to account for differences in Lagrangian and Eulerian descriptions of stress and strain. It is clear the NNCM inputs must be changed. Our current work is toward the goal of identifying a NNCM architecture capable of learning the stress-strain response of hyperelastic materials.
Projecting further, adjustments to the CaNNCM architecture will also be necessary to characterize a heterogeneous distribution of non-linear, viscoelastic materials. If the operating principle of linear-elastic CaNNCMs is retained, the stress-strain response would take the form σ NN NN(ε/S ε x , t), where spatial information is again encoded within the spatial scaling values S ε x . Less obvious from this relationship is the necessity of stressstrain state information at the spatial network input to account for the non-linear and time-dependent stress responses. This requirement can be understood by expanding on the explanation provided in [24] for the operation of spatial scaling values. When all materials are linear-elastic, the stiffness matrix D NN (recall its use in Eq. 18) computed from the weights of the MPN is constant. The role of S ε x is to vary the effective "stiffness" described by D NN based on position. Because the materials considered were linearelastic, the modified stiffness matrix is independent of time and state of deformation, meaning S ε x is constant at x. Conversely, D NN computed from a non-linear MPN will not be constant, nor will the rate at which it changes be the same at all points. Therefore, values of S ε x will need to change based not only on position, but also the current (or previous) stress-strain states.

DISCUSSION
Quantitative elasticity imaging is a challenging problem. Much of the difficulty arises from the limited measurement data that can be acquired during a clinical US imaging exam. Model-based inverse methods address the ill-posed inverse problem by constraining the solutions through geometric and constitutive model assumptions, well-defined force loads, and application of regularization to the solution. While effective, deviation of measurement data from the model assumptions can introduce large artifacts in the reconstructed material property images, potentially corrupting diagnostic information, and the regularization method is almost never grounded in physical laws. Furthermore, parameters associated with the pre-selected constitutive model may only capture a portion of the tissue mechanical properties, missing information that could be important to clinical decision making.
AutoP and NNCMs are a fundamentally different approach to the inverse problem in elasticity imaging. Rather than constraining the solution through model assumptions, AutoP extracts material property information from force-displacement measurements through a combination of physical modeling and machine learning. In this review, we presented AutoP in a new way to highlight similarities of the procedure with well-known optimization methods. Specifically, we demonstrated how the principles of stress equilibrium and strain compatibility are exploited through application of measurement data as BCs in separate FEAs, from which naturally arises an implicit displacement error that drives NNCM training.
The primary benefit of AutoP over model-based inverse methods is the potential to discover 3-D mechanical properties of soft tissues in vivo. Full computational modeling of applied force loads means measurements need not be acquired in a constrained manner; e.g., uniaxial compression by a large platen. Nor is 3-D displacement data necessary to reconstruct material properties throughout a volume, as we demonstrated in [25]. However, the use of arbitrary force-displacement measurements as input to AutoP significantly increases the computational load. Furthermore, we expect intelligent sampling strategies will need to be developed in order to maximize the information content of measurement data and reduce the computational burden of AutoP. Good sampling strategies will likely be essential for effectively imaging nonlinear, time-dependent material properties.
One particularly exciting prospect of CaNNCMs and AutoP is the ability to identify mechanical properties from a heterogeneous distribution of material types. Classic modelbased inverse methods must first assume the form of the constitutive model for the entire tissue section. The solution is the spatial variation of corresponding parameters that best fits the measurement data. AutoP reverses this order to first build a non-parametric, soft-computational model that describes the spatially varying material properties encoded in forcedisplacement data. Parameter estimation occurs after relevant information has been extracted from measurements. Successfully developing this capability would represent a paradigm shift in not only QUSE, but inverse identification of material properties in general.
AutoP is at its core a method to solve boundary value inverse problems. It is naturally suited for elasticity imaging because the physical principles of stress equilibrium and strain compatibility are built into the forward model. These principles relate the available force-displacement measurement data to the desired stress and strain distributions. By combining the physical principles through forward modeling with the flexibility of ANNs, virtually any type of continuum material can be described non-parametrically. AutoP could be applicable to other inverse problems with well-defined physical laws relating measurement data to the desired quantities. Nevertheless, AutoP is just one example of integrating physical modeling with machine learning to build powerful methods of solving inverse problems.

AUTHOR CONTRIBUTIONS
The manuscript was written by CH with substantial intellectual input from each author.

ACKNOWLEDGMENTS
Research reported in this publication was supported by NCI and NIBIB of the National Institutes of Health under Award Numbers R01 CA168575 and R21 EB023402. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.