BME 552: Computational Biomechanics
Estimated study time: 9 minutes
Table of contents
Sources and References
Primary texts — Zienkiewicz, Taylor, and Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th ed. (Butterworth-Heinemann). Humphrey and O’Rourke, An Introduction to Biomechanics: Solids and Fluids, Analysis and Design, 2nd ed. (Springer).
Supplementary texts — Winter, Biomechanics and Motor Control of Human Movement, 4th ed. (Wiley). Pandy, Computer Modeling and Simulation of Human Movement (Annual Reviews). Bonet and Wood, Nonlinear Continuum Mechanics for Finite Element Analysis, 2nd ed. (Cambridge).
Online resources — OpenSim (Stanford) open-source musculoskeletal modelling with tutorials. FEBio open-source nonlinear FE for biomechanics. SimTK open repository of biomechanical models. MIT OCW 2.080J Structural Mechanics. Mark Panjabi lecture archives on spine biomechanics (open).
Chapter 1: Foundations
1.1 Why Computational
Biological tissues are geometrically complex, mechanically nonlinear, and coupled across length scales. Analytical solutions exist only for idealized geometries and constitutive laws. Computational methods — finite element, multibody dynamics, smoothed-particle hydrodynamics — extend analysis to the realistic geometries and loading conditions needed for device design, surgical planning, and injury prediction.
1.2 Continuum Mechanics Primer
Deformation is described by the gradient \( \mathbf{F} = \partial \mathbf{x}/\partial \mathbf{X} \), with Green–Lagrange strain \( \mathbf{E} = \tfrac{1}{2}(\mathbf{F}^T\mathbf{F} - \mathbf{I}) \). Second Piola–Kirchhoff stress \( \mathbf{S} \) is work-conjugate to \( \mathbf{E} \). Linear elasticity specializes to small strain \( \boldsymbol{\varepsilon} = \tfrac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T) \) with Cauchy stress \( \boldsymbol{\sigma} = \mathbf{C}:\boldsymbol{\varepsilon} \).
Chapter 2: Finite Element Method
2.1 Weak Form and Discretization
Starting from equilibrium \( \nabla\cdot\boldsymbol{\sigma} + \mathbf{b} = 0 \), the principle of virtual work gives the weak form
\[ \int_\Omega \boldsymbol{\sigma}:\delta\boldsymbol{\varepsilon}\, dV = \int_\Omega \mathbf{b}\cdot\delta\mathbf{u}\, dV + \int_{\partial\Omega_t} \mathbf{t}\cdot\delta\mathbf{u}\, dS . \]Discretizing \( \mathbf{u} = \sum_i N_i \mathbf{u}_i \) with shape functions \( N_i \) yields the element equations, assembled into the global system \( \mathbf{K}\mathbf{u} = \mathbf{f} \).
2.2 Elements and Integration
Common elements are tetrahedra (4-node linear, 10-node quadratic) and hexahedra (8-node linear, 20-node quadratic). Gauss quadrature integrates element stiffness. Full integration can lock in incompressible or bending-dominated problems; reduced integration with hourglass control is the pragmatic remedy.
2.3 Nonlinearities
Geometric nonlinearity (large deformation), material nonlinearity (hyperelastic, viscoelastic, plastic), and contact nonlinearity (unknown contact area) dominate biomechanics. Newton–Raphson iteration, arc-length continuation, and augmented-Lagrangian contact enforcement solve the resulting nonlinear systems.
Chapter 3: Soft-Tissue Mechanics
3.1 Hyperelasticity
Incompressible or nearly incompressible soft tissues are modelled by strain-energy density \( W(\mathbf{C}) \) functions. Neo-Hookean, Mooney–Rivlin, Ogden, and Holzapfel–Gasser–Ogden (for collagen-reinforced tissues) provide increasing fidelity:
\[ W_{\text{HGO}} = \frac{\mu}{2}(I_1 - 3) + \frac{k_1}{2k_2}\sum_{i=4,6}\!\left\{\exp\!\left[k_2(\bar{I}_i - 1)^2\right] - 1\right\} . \]Parameters are fit to uniaxial and biaxial test data; extrapolation beyond tested deformation ranges is hazardous.
3.2 Viscoelasticity and Poroelasticity
Generalized Maxwell models represent time-dependent response through Prony series \( G(t) = G_\infty + \sum_i G_i e^{-t/\tau_i} \). Poroelasticity (Biot theory) captures fluid-saturated tissues such as cartilage and intervertebral disc. Biphasic models predict creep and stress relaxation responses relevant to joint mechanics and implant loading.
3.3 Anisotropy
Fibre-reinforced tissues — tendon, ligament, myocardium, arterial wall — exhibit strong anisotropy. Constitutive models include fibre distributions through structure tensors or probability functions. Imaging (diffusion MRI, histology) supplies fibre orientation data that personalize the model.
Chapter 4: Musculoskeletal Multibody Dynamics
4.1 Rigid-Body Multibody
Body segments modelled as rigid bodies connected by joints yield equations of motion
\[ \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) + \mathbf{G}(\mathbf{q}) = \boldsymbol{\tau} + \mathbf{J}^T\mathbf{F}_{\text{ext}} . \]OpenSim assembles such models with joint coordinates, inertial properties from regression or imaging, and generalized forces.
4.2 Muscle-Tendon Actuators
Hill-type muscle models (BME 355) supply forces as functions of activation, length, and velocity. Wrapping surfaces route muscle lines of action around joints and bony landmarks. Moment arms are computed by tendon-excursion or partial derivatives of length with respect to joint angle.
4.3 Static and Dynamic Optimization
Inverse dynamics computes net joint moments from kinematics and external forces. Static optimization distributes moments among redundant muscles by minimizing
\[ J = \sum_i a_i^p, \quad p = 2\text{--}3 , \]subject to physiological force limits. Dynamic optimization extends to time-dependent problems, solving for muscle excitations that produce a target motion.
Chapter 5: Applications
5.1 Injury Simulation
Car-crash and sports-impact simulations use full-body FE models (GHBMC, THUMS) with validated head, neck, thorax, and extremity components. Injury criteria — HIC for head, Nij for neck, chest compression — are computed from model outputs and correlated with injury probability curves.
5.2 Orthopaedic and Cardiovascular Implants
Implant design uses FE to verify strength, fatigue, and biomechanical interaction with surrounding tissue. Stress shielding in total hip replacement is minimized by porous or lattice stems that reduce stiffness mismatch. Stent deployment, coronary hemodynamics, and heart-valve performance likewise benefit from validated models.
5.3 Surgical Planning and Rehabilitation
Patient-specific models built from imaging inform osteotomy planning, tumor resection, and orthognathic surgery. Rehabilitation modelling evaluates exoskeleton assistance, prosthesis alignment, and movement retraining by simulating predicted outcomes before patient trial.
Chapter 6: Verification, Validation, and Uncertainty
6.1 Verification vs Validation
Verification ensures the code solves the chosen mathematical model correctly — mesh convergence studies, manufactured solutions, benchmark problems. Validation ensures the model represents the physical system — comparison with experiment on independent quantities, under the same loading, with quantified error bars. The ASME V&V 40 standard codifies this distinction for medical devices.
6.2 Sensitivity and Uncertainty
Biological inputs — tissue properties, geometry, boundary conditions, loads — are uncertain and patient-variable. Local sensitivity gives gradient information; global sensitivity (Sobol, Morris) reveals interactions. Monte Carlo simulation with sampled inputs produces credible intervals on model predictions, the appropriate currency for clinical use.
6.3 Credibility Evidence
The FDA’s evolving guidance on computational modelling in regulatory submissions expects credibility evidence matched to model risk. The process combines V&V, applicability assessment (is the validated model applicable to the claim being made?), and uncertainty quantification. A computational model without credibility evidence is not usable for regulatory decisions regardless of its apparent sophistication.
6.4 Workflow
A complete project: define clinical question; acquire imaging and loading data; build geometry (segmentation, meshing); select constitutive models with calibrated parameters; impose boundary conditions; solve; post-process to clinically relevant quantities; verify mesh convergence; validate against independent data; propagate uncertainty; report results with explicit domain of applicability. Tools change; the workflow endures.
