Deformable Models for Surgical Simulation: A Survey

A. Challenge of Realistic Simulation

• Patient-specific tissue properties, incomplete automation of image segmentation and mesh generation, ill-conditioned mechanics, and complex boundary conditions complicate realistic simulation.Anisotropic, heterogeneous, nearly incompressible tissues contribute to numerical expense, inaccuracy, or instability.Source: A. Challenge of Realistic Simulation · p. 1

B. Challenge of Real-Time Simulation

• Visual feedback requires at least 30 Hz and haptic feedback at least 1,000 Hz, while simplifying computation to meet real-time demands can reduce accuracy.The realism and real-time requirements affect one another: improving one is mainly achieved at the detriment of the other.Source: B. Challenge of Real-Time Simulation · p. 1

II. DEFORMABLE MODELS FOR SURGICAL SIMULATION

• The survey covers mechanical-event soft-tissue deformation and divides models into heuristic, continuum-mechanical, and other methodologies.Its scope excludes deformation induced by thermal or electrical events, though these may be simulated alongside mechanical events at added computational cost.Source: B. Challenge of Real-Time Simulation · p. 1, §II · p. 2

A. Heuristic Modeling Methodology

• Geometrically based methods such as FFD and deformable splines enable fast deformation, but their control-point manipulation does not represent underlying tissue mechanics.Their limited physical realism led to their being mainly superseded by physically based models.Source: A. Heuristic Modeling Methodology · p. 2
• MSM models tissue as masses linked by springs, with nodal forces summed from connected springs and positions updated through time-stepping.Its simple implementation and computational efficiency support interactive simulation; spring stiffness, penalty forces, and dampers can represent heterogeneity, near incompressibility, and viscoelasticity, respectively.Source: Figure 1 · p. 2, Figure 2 · p. 3, A. Heuristic Modeling Methodology · p. 3
• MSM is generally non-convergent under mesh refinement, can introduce artificial anisotropy and heterogeneity, and has parameter fits that may fail when topology or boundary conditions change.Discrete spring parameters do not map simply to tissue constitutive laws, making nonlinear stress-strain behavior difficult to reproduce accurately.Source: A. Heuristic Modeling Methodology · pp. 3–4
• ChainMail adjusts voxel-positioned chain elements when neighboring geometric constraints are violated, enabling interactive deformation of medical volumes containing millions of voxels.It can simulate nonlinear force-displacement behavior, hysteresis, and stress relaxation, but constraint parameters are empirical and dynamic behavior is difficult to realize.Source: A. Heuristic Modeling Methodology · pp. 3–4, Figure 3 · p. 4
• Shape matching also manipulates point positions to meet geometric constraints, but realistic deformation depends on choosing cluster stiffness without accounting for tissue material properties.Heuristic models share ambiguity between discrete model parameters and tissue constitutive laws, often requiring optimization whose results may not transfer across changed conditions.Source: A. Heuristic Modeling Methodology · p. 4

1.1) Finite Element Method

• FEM offers physical accuracy by assembling element equations to determine nodal displacements, but its high computational cost challenges interactive deformation.The mesh must conform to the problem domain, and mesh-based solutions can fail when large deformation inverts elements.Source: B. Continuum-Mechanical Methodology · pp. 5, 8

1.2) Simplification of FEM

• FEM simplifications reduce computational complexity through lumped explicit FEM, matrix condensation, precomputation, and surface-only BEM, but may compromise accuracy or restrict material scope.BEM assumes an isotropic, homogeneous interior; strain-limiting FEM reduces equation solving to geometric constraints but confines element deformation range.Source: B. Continuum-Mechanical Methodology · pp. 5–6

1.3) Total Lagrangian Formulation

• The total Lagrangian formulation permits precomputation of spatial derivatives; for 2,535 nodes and 2,200 hexahedrons under ellipsoid indentation, it solves in 2.1 ms versus 10.6 ms for the updated formulation.GPU implementations of TLED are reported at speed gains up to 16.8x and 56.3x over equivalent CPU implementations.Source: B. Continuum-Mechanical Methodology · p. 6

1.4) Model Reduction

• Model reduction projects full system response into a smaller-dimensional subspace to reduce deformation degrees of freedom, with offline response extraction and online reduced-model simulation.The survey describes POD, nonlinear POD with ANM, PGD, and DEIM; PGD obtains responses at kilohertz rates, and DEIM gives a 10x speed gain over classical POD without empirical interpolation.Source: B. Continuum-Mechanical Methodology · p. 7

1.5) Element-Related Issues

• Low-order elements improve computational efficiency but can exhibit numerical deficiencies: under-integrated hexahedrons have a zero-energy mode, and tetrahedrons can lock for nearly incompressible materials.The survey describes hourglass control and improved average nodal pressure (IANP) as approaches to these problems.Source: B. Continuum-Mechanical Methodology · p. 7, Figure 5 · p. 7, Figure 6 · p. 7

2.2) Meshless Total Lagrangian Explicit Dynamics

• Meshless total Lagrangian explicit dynamics (MTLED) uses arbitrarily dispersed particles and background integration cells, running at half the speed of hexahedral TLED but three times faster than similar tetrahedral TLED.The method also faces challenges in numerical integration and boundary conditions; tetrahedral background cells address volume inaccuracy, while displacement corrections can enforce essential boundary conditions.Source: B. Continuum-Mechanical Methodology · p. 8

2.3) Smooth Particle Hydrodynamics and Point-Collocation-based Method of Finite Spheres

• SPH and PCMFS interpolate particle state variables over a smoothing distance using a kernel; SPH variants address incompressibility or model damage and failure.DEM requires calibration of microscopic parameters, while Peridynamics can model soft materials with any Poisson’s ratio.Source: B. Continuum-Mechanical Methodology · p. 8, Figure 8 · p. 9

2) Meshless Approach

• Meshless methods avoid mesh-topology degradation at large deformation, but can struggle in sparsely sampled regions and depend heavily on node placement for physical accuracy.MTLED is also unsuitable when a single node’s reaction force or displacement is needed, despite accuracy in overall reaction forces.Source: B. Continuum-Mechanical Methodology · pp. 8–9

4) Others

• Machine-learning models map inputs to deformation outputs; one reported real-time approach achieves position errors below 0.3 mm, while data-driven stress–strain relationships also offer a route to parameter identification.The RBFN approach trains offline on FEM responses to prescribed displacements and reconstructs deformation online; performance depends on training data and the chosen learning algorithm.Source: C. Other Modeling Methodologies · p. 9
• Costa’s fast deformation model uses Pascal’s principle and volume conservation for soft tissues formed by fibers and fluid, but applies only to fluid-filled objects and has no dynamic behavior.The model represents near incompressibility through volume conservation.Source: C. Other Modeling Methodologies · p. 9

D. Summary of Representative Deformable Models

• Table I frames the comparison by model, year, hardware, speed, and application.The passage introducing the table says it summarizes representative models alongside advent times, hardware settings, and surgical applications.Source: D. Summary of Representative… · p. 9
• The listed TLED entry gives the year 2007 and Intel Pentium hardware.Source: D. Summary of Representative… · p. 10

A. Linear Elasticity

• Most existing surgical-simulation methods use linear elasticity, which assumes linear stress–strain behavior and infinitesimal strains.For isotropic, homogeneous materials, the constant material matrix is defined by the two Lamé coefficients.Source: A. Linear Elasticity · p. 10
• Linear elasticity relates strain to displacement through a shape function, and finite-element assembly yields the global equilibrium system K u = f.The global system uses a 3n × 3n stiffness matrix and 3n × 1 displacement and force vectors, where n is the number of nodes.Source: A. Linear Elasticity · p. 10
• Linear elasticity is efficient because its global stiffness matrix remains constant and can be precomputed, but its solution is valid only for small deformations.Large rotations can produce unrealistic volume growth.Source: A. Linear Elasticity · p. 10
• Geometric nonlinear models can handle large-deformation geometry but retain linear material laws, so they do not capture material nonlinearity.Changing Lamé coefficients to model nonlinear material behavior does not comply with constitutive laws because they are inherent material attributes.Source: A. Linear Elasticity · pp. 10–11

B. Nonlinear Hyperelasticity

• Total Lagrangian formulations represent stress and strain relative to the initial configuration using second Piola-Kirchhoff stress and Green-Saint Venant strain.The deformation gradient maps the initial particle position to its deformed position.Source: B. Nonlinear Hyperelasticity · p. 11
• The strain energy is expressed as a function of Green-Saint Venant strain, and the second Piola-Kirchhoff stress is derived from that energy.For isotropic, homogeneous material, strain energy can be expressed as a function of the right Cauchy-Green tensor invariants.Source: B. Nonlinear Hyperelasticity · p. 11
• Studied total Lagrangian approaches for large soft-tissue deformation include TLED, implicit FEM, MTLED, MJED, and total Lagrangian SPH.The review also points to further details on nonlinear hyperelasticity.Source: B. Nonlinear Hyperelasticity · p. 11

IV. COMPARISON OF INTERNAL FORCES

• The dynamic balance is m ü + f_int = f_ext, and internal forces affect both model accuracy and computational complexity.Table II compares internal forces among different deformable models.Source: §IV · p. 11, Table II · p. 11
• Mass-spring models define forces through springs between neighboring mass points, while linear FEM uses stiffness, displacement, shape-function, and material matrices.The spring description includes spring stiffness and rest length.Source: §IV · p. 12
• Co-rotational FEM updates an element rotation matrix at each time step, while linear and nonlinear TMM use stiffness tensors or parameters.TLED uses a vector form of second Piola-Kirchhoff stress; MJED decomposes strain energy into terms involving the Jacobian and tensor invariants.Source: §IV · p. 12
• MTLED distributes integration points over the domain, normalized total Lagrangian SPH uses particle volumes and kernel gradients, and fibers-fluid modeling includes pressure and fiber forces.The listed formulations specify their integration-point weights, particle reference distances, and fiber-force terms.Source: §IV · p. 12

V. ANISOTROPY, VISCOELASTICITY, AND COMPRESSIBILITY

• Soft tissues’ anisotropy, viscoelasticity, and near-incompressibility are properties that deformable models need to consider for realistic simulation.Fiber microstructure and functional components contribute to directional behavior, while stress relaxation, creep, and hysteresis characterize viscoelasticity.Source: §V · p. 12

A. Anisotropy

• Heuristic models assign different component parameters to represent anisotropy, while continuum methods modify the strain energy density function to represent directional dependence.The survey cautions that discrete heuristic representations produce only arbitrary anisotropic behaviors.Source: A. Anisotropy · p. 12
• For transversely isotropic materials, the initial fiber direction a_0 introduces pseudo-invariants I_4 and I_5 into the strain energy density function.I_1, I_2, and I_3 are invariants of C; I_4 and I_5 arise from anisotropy introduced by a_0.Source: A. Anisotropy · p. 12
• The fiber-induced anisotropy contributes additional I_4- and I_5-dependent terms to the second Piola-Kirchhoff stress S.Source: A. Anisotropy · p. 12
• Orthotropic materials use two initial fiber directions, a_0 and b_0, with a third mutually orthogonal preferred direction, and add pseudo-invariants I_6 and I_7.The strain energy density function includes I_1 through I_7.Source: A. Anisotropy · p. 12

B. Viscoelasticity

• Continuum viscoelastic models use a time-dependent strain energy density function expressed as a convolution integral.The survey also lists Maxwell, Kelvin-Voigt, standard linear solid, generalized Wiechert, and double Maxwell-arm Wiechert models; these comprise linear dashpots and springs.Source: B. Viscoelasticity · p. 13
• The relaxation function is a generalized Maxwell model expressed as a Prony series with positive constants, and its coefficients can be normalized so the function's initial value is one.Source: B. Viscoelasticity · p. 13
• The time-dependent strain energy density modifies the second Piola-Kirchhoff stress into a time-dependent stress tensor, which is discretized into an incremental update for each time step.The update uses constant coefficients A_i and B_i determined from the time increment and relaxation parameters.Source: B. Viscoelasticity · p. 13

C. Compressibility

• For isotropic soft tissues with some compressibility, the strain energy density separates into isochoric and volumetric parts.The isochoric invariants are based on the isochoric right Cauchy-Green tensor C̅ = J^-2/3 C, and J is the volume ratio.Source: C. Compressibility · p. 13

VI. NUMERICAL TIME INTEGRATION

• Soft-tissue dynamics are integrated by converting the second-order equation into two first-order systems using a proxy velocity, then discretizing velocity and acceleration with finite differences.Forward or backward finite-difference estimates yield explicit or implicit integration schemes.Source: §VI · p. 13

A. Explicit Integration

• Explicit integration is computationally efficient and easy to implement because it advances from known current values without inverting the stiffness matrix at each time step.Mass lumping can also split the global system into independent node equations for parallel computation.Source: A. Explicit Integration · p. 14
• Explicit integration is only conditionally stable, so its time step must satisfy the CFL condition and can be restricted by near-incompressibility and viscoelasticity.The resulting small time steps require more iterations per simulation frame; proposed techniques can enlarge the step, but do not eliminate conditional stability.Source: A. Explicit Integration · p. 14

B. Implicit Integration

• Implicit integration is unconditionally stable and can use large time steps for collisions and stiff near-incompressible tissue, but each load step requires solving a nonlinear system.Its time step remains limited by numerical convergence and accuracy, and its computation per iteration is substantially greater than explicit integration.Source: B. Implicit Integration · pp. 14–15
• Explicit integration avoids solving a large system and is efficient, whereas implicit integration has a more global response and better controls equilibrium accuracy in tool-tissue interactions.Explicit integration can accumulate time-discretization error and, for large deformation, distorted elements can reduce its critical time step; dynamic relaxation improves convergence at the cost of accuracy.Source: B. Implicit Integration · p. 15

• No single deformable model addresses both realistic and real-time surgical simulation; different approaches were developed to meet specific needs.Source: §VII · p. 15
• Geometric heuristics favor speed over physical accuracy; MSM optimization can be slow or fail when its assumed conditions change, while ChainMail supports interactive deformation of large medical volumes.Source: §VII · p. 15
• Shape matching offers unconditional stability, but these models struggle to associate parameters with tissue constitutive laws and therefore produce physically plausible simulations.Source: §VII · p. 15
• Continuum methods can improve accuracy at greater computational cost; FEM simplifications and precomputation constrain accuracy or topology changes, while meshless methods avoid element distortion but depend on particle distribution and represent local forces less accurately.Source: §VII · pp. 15–16
• Machine-learning models encode tissue behavior in trained network coefficients, but their simulated results depend strongly on the chosen learning algorithm.Source: §VII · p. 16
• Table III compares the main deformation approaches by physical accuracy and computational performance as a reference for selecting a strategy for different simulation conditions.Source: §VII · p. 16, Table III · p. 16
• Cutting and tearing can require removing or splitting elements, using X-FEM, or modeling meshless discontinuities; resulting stiffness updates may degrade explicit integration's critical time step.Source: §VII · p. 16
• The Truth Cube provides physical displacement data from compression and indentation tests, but its regular geometry and well-characterized conditions do not capture surgical complexity.Source: §VII · p. 16, Figure 9 · p. 16

• The survey classifies deformable models into heuristic, continuum-mechanical, and other methodologies, and reviews linear and nonlinear modeling, internal forces, time integration, and complex biomechanical behaviors.Source: §VIII · p. 17
• Future modeling should link passive mechanical deformation with active electrical activity and incorporate patient-specific tissue-property measurements rather than assuming properties are fixed and known.Source: §VIII · pp. 17–18
• Integrating soft-tissue models into robotic control loops could support automatic planning and precise control of robotic surgical tasks.Source: §VIII · p. 18