A Numerical PDEs Approach to Evolution Equations in Shape Analysis Based on Regularized Morphoelasticity

Abstract: This work studies a variational formulation and numerical solution of a regularized morphoelasticity problem of shape evolution. The foundation of our analysis is based on the governing equations of linear elasticity, extended to account for volumetric growth. In the morphoelastic framework, the total deformation is decomposed into an elastic component and a growth component, represented by a growth tensor $G$. While the forward one-step problem -- computing displacement given a growth tensor -- is well-established, a more challenging and relevant question in biological modeling is the inverse problem in a continuous sense. While this problem is fundamentally ill-posed without additional constraints, we will explore parametrized growth models inscribed within an optimal control problem inspired by the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework. By treating the growth process as a path within a shape space, we can define a physically meaningful metric and seek the most plausible, energy-efficient trajectory between configurations. In the construction, a high-order regularization term is introduced. This elevates the governing equations to a high-order elliptic system, ensuring the existence of a smooth solution. This dissertation focuses on the issue of solving this equation efficiently, as this is a key requirement for the feasibility of the overall approach. This will be achieved with the help of finite element solvers, notably from the FEniCSx library in Python. Also, we implement a Mixed Finite Element Method, which decomposes the problem into a system of coupled second-order equations as a treatment of these high-order systems that have significant computational challenges.

• The paper introduces a novel variational PDE framework that combines linear elasticity with additive growth to model shape evolution.
• It employs a mixed finite element method alongside high-order regularization to ensure numerical stability for complex deformations.
• The study tackles both forward and inverse problems, paving the way for physics-informed neural networks in rapid morphoelastic inference.

A Numerical PDEs Approach to Evolution Equations in Shape Analysis Based on Regularized Morphoelasticity

Introduction

This work develops a rigorous variational and numerical framework for evolution equations in geometric shape analysis derived from biomechanical morphoelasticity. Rather than depending on image-based or Eulerian grid-based approaches, the methodology formulates shape evolution via continuum mechanics, primarily exploiting Lagrangian representations. The fundamental modeling strategy is to merge linear elasticity with volumetric growth, articulated through an additive/symmetric growth tensor $\mathbf{G}$, and to embed this into a high-order regularized dynamical system for applications in computational anatomy and shape statistics.

The forward and inverse evolution problems are considered. For forward modeling, given a growth distribution, the work delivers a variational PDE (augmented by appropriate regularization) and addresses computational challenges using a Mixed Finite Element Method (MFEM) with the FEniCSx platform. From the inverse perspective, the work exposes the ill-posedness of the general growth identification problem and situates a solution in a physically motivated optimal control/geodesic framework, leveraging the Large Deformation Diffeomorphic Metric Mapping (LDDMM) theory to define efficient Riemannian metrics on non-linear shape spaces.

Physical and Mathematical Foundations

Linear Elasticity and Lagrangian Formulation

The modeling base is the linear elasticity PDE with isotropic, homogeneous constitutive response, small strain, and small displacement assumptions. The variational (weak) form is derived directly from the principle of minimum potential energy to guarantee coercivity and positive definiteness, properties that are essential for establishing metric structure on the shape manifold. These concepts, well-established in elasticity theory [gould2018introduction, sadd2020elasticity], ensure the well-posedness of the static problem and form the bedrock for subsequent extensions.

Deformation is visualized via direct mapping of the domain $\Omega$ by the displacement field $\mathbf{u}$, with the configurations displayed as pre- and post-deformation meshes.

Figure 2: Deformed configuration ($\Omega_1$) of the cantilever beam under gravitational loading, highlighting significant bending in the $-z$ direction.

Morphoelasticity and Additive Growth Modeling

To model biological growth or resorption, the elastic strain is additively decomposed into an elastic and growth component as $$\varepsilon_g = \varepsilon(\mathbf{u}) - \mathbf{G}$$. The stress responds only to the elastic fraction, and the growth tensor $\mathbf{G}$ is often specified as an isotropic Gaussian spatial field centered in the reference domain.

The new variational problem possesses identical bilinear (elastic) structure but incorporates an additional linear term due to growth. Importantly, for biologically sensible and numerically tractable models, the growth field must be parametrized and tied to Lagrangian material points, not just Eulerian coordinates.

Figure 4: Reference configuration ($\Omega_0$) of the sphere for Example 2.

Figure 1: Deformed configuration ($\Omega_1$) of the sphere under a surface traction.

Simulation Scenarios and Growth Dynamics

The paper provides multifaceted simulation results to demonstrate modeling capabilities:

• Cantilever beams and spheres with/without gravity and surface pressure.
• Volumetric growth: Both positive (expansion) and negative (atrophy) cases demonstrate the flexibility of the morphoelastic framework.
• Competition between growth and load: By varying the growth amplitude $A_g$, the transition from external force-dominated to growth-dominated deformation is quantified.

  Figure 7: Reference configuration ($\Omega$0) for Model 3: rectangular cantilever beam.

  

  Figure 3: Deformed configuration ($\Omega$1) of the beam, exhibiting gravity-driven bending plus localized internal growth expansion.

  

  Figure 5: Deformed configuration of the sphere with low growth, exhibiting compression dominated by external pressure.

  

  Figure 11: Intermediate growth amplitude counteracts pressure, inducing visible upward displacement.

  

  Figure 6: High growth amplitude overcomes external load, leading to volumetric expansion.

These results affirm that the model qualitatively captures phenomena such as threshold behavior and competition between internal/external drivers.

The Inverse Problem and Shape Space Geometry

Ill-posedness and Regularization

Given a displacement field $\Omega$2 (typical in biological imaging), the task of inferring the growth tensor $\Omega$3 is generally non-identifiable: the null space of the divergence of the stress operator is rich, and infinitely many distinct growth fields may generate the same observed deformation. Only by restricting $\Omega$4 (e.g., to isotropic, parametric families) can the inverse mapping be rendered unique.

Geometric Framework and Riemannian Metric

The solution is to move from a static, one-shot mapping to a time-dependent, optimal control path within the infinite-dimensional shape space. Charon and Younes [charon2023shape] formalized this by using elastic energy dissipation as a Riemannian metric on shape space manifolds, yielding well-defined geodesic distances. The physical cost of a trajectory (interpreted as minimum total potential energy dissipated via growth/remodeling) now rigorously encodes anatomical evolution.

The model is discretized for computation as a time-stepping sequence—each step involves updating the growth field in the Lagrangian frame, solving a regularized equilibrium, and evolving the reference domain.

High-Order Regularization and Mixed Finite Element Implementation

LDDMM Regularization

To ensure solution smoothness, avoid mesh entanglement, and guarantee diffeomorphic evolution, a high-order Sobolev penalization term is introduced. The regularization operator $\Omega$5 (with $\Omega$6 in 3D) draws directly from LDDMM theory [beg2005computing], enforcing that the displacement field possesses sufficient regularity (in $\Omega$7, so $\Omega$8 or better).

This changes the governing PDE system from second to higher-order (up to eighth-order in practical examples). Standard finite elements cannot handle this directly, so a mixed formulation is used wherein a sequence of auxiliary variables decouple the operator into solvable second-order blocks.

Computational Results

The solver is demonstrated across a spectrum of shapes: convex/concave, with varied boundary constraints and both positive and negative growth.

Figure 8: Initial shape for simulation with localized growth.

Figure 15: Initial ellipsoid geometry, facilitating asymmetric growth-induced bulges.

Figure 17: Kidney-shaped (non-convex, smooth) domain for morphoelastic simulation.

Figure 9: Initial sphere geometry for eccentric/negative growth cases.

Figure 10: Initial configuration exhibiting non-trivial geometric features.

The solver successfully produces plausible deformations without mesh self-intersection or collapse, and the computational cost is demonstrated to be practical (∼ seconds for moderate mesh sizes and 4 steps).

Limitations and Future Developments

The computational complexity increases significantly for the inverse problem, particularly with high-order regularization and large deformation. Traditional PDE-constrained optimization/adjoint methods are not tractable for these scenarios. To overcome this, the paper advocates for Physics-Informed Neural Networks (PINNs) and Neural Operator architectures (e.g., DeepONet [Lu2021]) to bypass mesh generation and directly encode PDE constraints via automatic differentiation. This reformulates both forward and inverse morphoelastic problems as mesh-free, functional regressions, enabling substantially faster inference and extending to real-time applications.

The extension to operator learning frameworks portends rapid surrogates for geodesic shooting in shape space, dramatically lowering the cost of statistical analysis and registration tasks in high-throughput biomedical settings.

Conclusion

This paper presents an advanced computational pipeline for modeling and simulating evolution equations in shape analysis based on regularized morphoelasticity. By marrying explicit physical modeling, smoothness-enforcing high-order regularization, and rigorous variational numerics, the work builds a practical yet mathematically principled bridge between continuum biomechanics and geometric data representation. The demonstration on varied shapes and load/growth scenarios underscores the scope and flexibility of the approach. Looking forward, incorporating data-driven solvers and neural operator models promises to enable efficient solutions to ill-posed inverse problems and support large-scale applications in computational anatomy and shape statistics.