Biophysical Modeling of Respiratory Organ Motion

) image sequences, often without concerning knowledge about anatomy and physiology in detail. In contrast, biophysical approaches aim at identification of anatomical and physiological aspects of breathing dynamics that are to be modeled. In the context of radiation therapy, biophysical modeling of respiratory organ motion commonly refers to the framework of continuum mechanics and elasticity theory, respectively. Underlying ideas and corresponding boundary value problems of those approaches are described in this chapter, along with a brief comparison to image registration-based motion field estimation.

4.1 Introduction

Approaches for estimation of respiratory organ motion in 4D image sequences can be divided into two main types: image registration and biophysical modeling [32, 56, 62]. Image registration aims at estimating motion fields by directly extracting them from the 4D image sequences, commonly by maximizing similarity measures between the individual 3D frames of the image sequences; knowledge about breathing physiology and physical properties of organs are (usually) not taken into account [25, 32]. In contrast, biophysical modeling means to identify aspects of the physiological processes and to describe them within a physics-based mathematical formulation. These approaches are also often referred to as biomechanical or biophysics-based (image) registration [4, 11, 44] and so contrasted to other types of image registration such as intensity- or feature-based. In addition, biophysical models are often directly associated with Finite Element Methods (FEM) because FEM are frequently applied to solve the problem formulation derived during the modeling process (usually partial differential equations, PDEs), especially for approaches within the framework of continuum mechanics. Strictly speaking, it should be noted that synonymous use of the terms FEM-based and biophysical modeling is misleading because FEM represent only a technique, amongst others, to numerically solve PDE systems. Therefore, subsequently only the term of biophysical motion modeling is used.

Biophysical modeling of breathing mechanics initiated long before the introduction of modern 4D imaging devices (see [29, 37] and references therein as examples). A diversity of approaches can be found in the literature and the spatial scales range from microscopic (e.g. analysis and simulation of alveolar structures in the lungs) to global ranges. Multi-scale models are also being developed, cf. the lung physiome project [16], but in the context of radiation therapy, only large scale models are commonly applied. These are models that assess macroscopic structure behavior and interactions [9]. Underlying techniques cover, for example, mass-spring models, in which organs are described by a set of distributed masses connected by springs [9, 30], or more generally particle systems, where the organ is represented by a collection of particles [18, 28]. For mass-spring models, the system behavior is mainly characterized by the spring parameters and applied external forces; using particle systems, organ deformation is controlled by specified physical particle properties, assumed interaction potentials and the external forces. Although these techniques tend to be intuitive and computationally efficient, they are not necessarily accurate. The behavior of mass-spring models, for instance, is described to be dependent on the topology of the springs and masses [20]. Further, proper values for the spring constants/particle properties are not always easy to derive and say little about the material being modeled [20, 40]. Therefore, the most prominent model formulations refer directly to the framework of continuum mechanics, and predominantly elasticity theory. The main ideas of these approaches are described and exemplified in this chapter.

The chapter is therefore organized as follows: First, models that focus on motion of a single organ are described, with the lungs and the liver as examples (Sect. 4.2). Extensions of the single organ models to provide motion estimation for several anatomical structures are presented in Sect. 4.3. In the last section, comparison studies of biophysical modeling approaches and image registration-based motion field estimation are briefly reflected. We close with conclusions on potential benefits by combining the two worlds.

Table 4.1

Notation of continuum mechanics/elasticity theory used in this chapter, mainly following the notation introduced in [10]

$\varGamma _i$	Undeformed/initial organ geometry for breathing phase
$\partial \varGamma _i$	Surface of $\varGamma _i$
$\varvec{n}$	Outer-pointing normal to the surface $\partial \varGamma _i$
$\varphi$	Injective and orientation preserving transformation $\varphi =\varvec{1}+\varvec{u}$ of organ $\varGamma _i$
$\varvec{1}$	Identity map
$\varvec{u}$	Displacement field parameterizing the transformation $\varphi$ . Typically, $\varvec{u}=\left( u_x,u_y,u_z\right) ^T$ represents the sought motion field estimation
$\hat{\varvec{u}}$	Prescribed displacements for $\varvec{x}\in \partial \varGamma _i$ (Dirichlet boundary conditions)
$\varvec{\varepsilon }$	Linear or infinitesimal strain tensor
$\varvec{E}$	Green-St.Venant (Green-Lagrange, Lagrange) strain tensor
$\varvec{e}$	Eulerian-Almansi strain tensor
$\varvec{\sigma }$	Cauchy or true stress tensor
$\varvec{S}$	First Piola-Kirchhoff stress tensor
$\varvec{\varSigma }$	Second Piola-Kirchhoff stress tensor
	Strain energy function. Used to define the stress–strain relation for hyperelastic materials
$\lambda$	First Lamé parameter; no physical interpretation
$\mu$	Second Lamé parameter, shear modulus; measure for the resistance to volume preserving shear deformations, ratio of shear stress to shear strain
	Young’s modulus; relative material stiffness within the elastic range, defined as ratio of uniaxial (tensile) stress to uniaxial (tensile) strain
$\nu$	Poisson’s ratio; measure of the compressibility, defined as ratio of relative transverse contraction to relative longitudinal stretching
$\varvec{f}$	Density of external body forces (e.g. gravity)
$\tilde{(\cdot )}$	$\tilde{(\cdot )}$ indicates that the term refers to the deformed organ geometry. Example: $\varGamma _i$ represents the undeformed organ geometry for breathing phase , whilst $\tilde{\varGamma }_i=\varphi \left( \varGamma _i\right)$ describes its deformed state

To describe the modeling ideas and corresponding boundary value problems, it is referred to the notation of elasticity theory. The main terms are summarized in Table 4.1, but please refer to corresponding textbooks like [10, 63] for detailed explanations and derivations.

4.2 Single Organ Motion as a Problem of Elasticity Theory

Modeling respiratory organ motion using elasticity theory is motivated by the view that organ movements can be described as elastic body deformations. Current literature on single organ modeling focuses primarily on two structures: the lungs as the central organ for breathing, and the liver as the largest abdominal organ.

4.2.1 Modeling Macroscopic Lung Motion

Modeling and simulation studies regarding lung elasticity can be found already in the 1970s (and most likely before) [38, 59]. The pioneer studies were based on only simplified lung shape-like geometry definitions due to the lack of modern (4D) tomographic imaging devices and limited computational power. The ideas can often be considered to be the ground of current biophysical breathing motion modeling attempts, but the models were rarely intended or suited to be used in radiation therapy context. This section therefore focuses on state-of-the-art approaches motivated by that application.

These approaches usually aim at explicit modeling of the process of lung ventilation. Conditions of lung deflation tend to be more sophisticated from the perspective of modeling (key words: passive process, retraction forces), and so the majority of the studies focuses rather on lung inflation than deflation. In the next paragraph, the physiological principles of inflation are briefly described first. Following, the modeling ideas are explained and the corresponding boundary value problems defined.

4.2.1.1 Background: Thoracic Anatomy and the Physiology of Lung Inflation

Anatomically, each lung is surrounded by a pleura sac, which is made up of two layers: parietal pleura and visceral pleura (Fig. 4.1). The parietal pleura is adherent to the internal surface of the thoracic cavity, the diaphragm, and the mediastinum. The visceral pleura covers the outer surface of the lungs. Both layers join at the root of the lung, which is the point of entry of the pulmonary vasculature, nerves, and airways into the lung. The space enclosed by the pleurae is the pleural cavity. It is filled with an incompressible fluid and subject to the so-called (intra-)pleural pressure, which is generally negative (relative to the atmospheric pressure) [27].

Fig. 4.1

Illustration of the anatomical terminology

In lung inflation, the breathing muscles (diaphragm, outer intercostal muscles) contract and the thoracic cavity expands. As a consequence, the intrapleural pressure takes larger negative values. This, in turn, forces parietal and visceral pleura to keep contact and finally the lungs to follow the expansion of the thoracic cavity. Driven by diaphragmatic breathing, the lung deformation occurs predominantly in superior-inferior direction and so—due to the liquid in the pleural gap—the visceral pleura slides down the parietal pleura. This behavior is the starting point of most current modeling approaches.

4.2.1.2 Lung Inflation as Elastic Boundary Value Problem

Despite detailed anatomical knowledge of the lungs: The majority of current studies on modeling macroscopic lung behavior approximate lung tissue to be homogeneous and isotropic [1, 6, 52, 55, 56, 62]. Inner lung structures like the bronchial tree and pulmonary vessels are ignored, each of which having in principle different elastic properties [3, 56]. Then, given the lung at its undeformed state (subsequently assumed to be end expiration, EE) and described by a closed and connected subset $\varGamma _{\mathrm{EE }}$ of $\mathbb R ^{3}$ , the sought injective and orientation preserving deformation is denoted by $\varphi :\varGamma _{\mathrm{EE }}\rightarrow \mathbb R ^{3}$ , which is parameterized by a displacement field (= motion field estimation) $\varvec{u}:\varGamma _{\mathrm{EE }}\rightarrow \mathbb R ^{3}$ and $\varphi =\mathbf{1 }+\varvec{u}$ .

Further, inertia effects are also usually neglected and a quasi-static formulation is employed [25, 51]. This can be seen as being motivated by assuming breathing frequencies to be only low and, consequently, the lungs to remain in quasi-static equilibrium—which, in turn, allows for an easier model formulation because the notations of elastostatics can be applied.

Starting with these approximations and simplifying assumptions, elastic boundary value problems (BVPs) are defined to determine the transformation or the displacement field within elasticity theory, including the specification of the governing equations described in the following: kinematic equations (strain-displacement relationship), constitutive equations (material properties), the equations of equilibrium, and appropriate boundary conditions.

Kinematic Equations (Strain-Displacement Relationship)

The kinematic equations describe the relationship between the displacements, interpreted as movements of the continuum or organ particles, and the resulting changes of the particle configuration, characterized by the regional distribution of strain tensors. The principle concept of strain is to measure how much a given displacement $\varvec{u}$ differs locally from a rigid body displacement. Strain tensor formulations are usually expressed in terms of the displacement gradient $\mathbf{\nabla }\varvec{u}$ , defined by

$\begin{aligned} \mathbf{\nabla }\varvec{u}=\left( \begin{array}{lll} \frac{\partial u_x}{\partial x} &{} \frac{\partial u_x}{\partial y} &{} \frac{\partial u_x}{\partial z}\\ \frac{\partial u_y}{\partial x} &{} \frac{\partial u_y}{\partial y} &{} \frac{\partial u_y}{\partial z}\\ \frac{\partial u_z}{\partial x} &{} \frac{\partial u_z}{\partial y} &{} \frac{\partial u_z}{\partial z}\\ \end{array} \right) \end{aligned}$

(4.1)

with

and

being the components of $\varvec{u}$ . In classical mechanics, different strain tensor formulations are in use, with the following being most common in the given context.¹

Infinitesimal Strain Tensor. The (Cauchy’s) infinitesimal or linear strain tensor is composed of the linear (normal) strain along , , and axes,

$\begin{aligned} \varepsilon _x=\frac{\partial u_x}{\partial x}, \quad \varepsilon _x=\frac{\partial u_y}{\partial y}, \quad \varepsilon _x=\frac{\partial u_z}{\partial z}, \end{aligned}$
and the (engineering) shear strain components,

$\begin{aligned} \gamma _{xy}=\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}, \ \gamma _{yz}=\frac{\partial u_y}{\partial z}+\frac{\partial u_z}{\partial y}, \ \gamma _{xz}=\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x}, \ \end{aligned}$
which describe the angular change at any point between two lines crossing this point before and after deformation. They are combined, resulting in the symmetric tensor

$\begin{aligned} \varvec{\varepsilon }=\left( \begin{array}{ccc} \varepsilon _x &{} \frac{\gamma _{xy}}{2} &{} \frac{\gamma _{xz}}{2}\\ \frac{\gamma _{yx}}{2} &{} \varepsilon _y &{} \frac{\gamma _{yz}}{2}\\ \frac{\gamma _{zx}}{2} &{} \frac{\gamma _{zy}}{2} &{} \varepsilon _z\\ \end{array} \right) =\frac{1}{2}\left( \mathbf{\nabla }\varvec{u}^T+\mathbf{\nabla }\varvec{u}\right) . \end{aligned}$

(4.2)
Note that the infinitesimal strain tensor is an approximate measure, assuming the first derivatives of the components being so small that higher order terms (squares, products of partial derivatives) are negligible [8, 10]. If larger shape changes are to be expected (i.e. undeformed and deformed configurations are assumed to be substantially different), strain tensors of finite strain theory (also called large deformation or large displacement theory) are applied.
Green-St.Venant Strain Tensor. Within finite strain theory, the Green-St.Venant Strain Tensor (also known as Green or Lagrangian finite strain tensor) measures the strain with respect to the undeformed geometry. It is defined by

$\begin{aligned} \varvec{E}= \frac{1}{2} \left( {\varvec{\nabla }}\varvec{u}^T+{\varvec{\nabla }}\varvec{u}+{\varvec{\nabla }}\varvec{u}^T{\varvec{\nabla }}\varvec{u}\right) . \end{aligned}$

(4.3)
and therefore introduces a non-linear strain-displacement relationship. As before, the off-diagonal elements of $\varvec{E}$ can be referred to as shear strains, the diagonal elements as normal strains, and it becomes obvious that Eq. 4.2 represents a linearization of Eq. 4.3; for further detailed descriptions of the underlying physical motivation and corresponding derivations please refer to, e.g., [8, 10].
Eulerian-Almansi Strain Tensor. The counterpart of the Green-St.Venant Strain Tensor within finite strain problem settings, but defined with reference to the deformed configuration, is given by the Almansi or Eulerian (finite) strain tensor,

$\begin{aligned} \varvec{e}=\frac{1}{2}\left( {\varvec{\nabla }}\tilde{\varvec{u}}^T+{\varvec{\nabla }}\tilde{\varvec{u}}-{\varvec{\nabla }}\tilde{\varvec{u}}^T{\varvec{\nabla }}\tilde{\varvec{u}}\right) . \end{aligned}$

(4.4)
Thereby and hereinafter, $\tilde{(\cdot )}$ indicates that the corresponding term (here: $\tilde{\varvec{u}}$ ) refers to the deformed geometry.

Focusing on respiratory lung motion, occurring displacements are generally considered to be too large to assume that changes in geometry do not influence the (local) lung behavior [52]. Thus, most modeling approaches assume geometrical non-linearity/a non-linear strain-displacement relation; they make use of finite strain theory and corresponding strain tensors.

Constitutive Equations (Material Properties)

The constitutive equations relate stresses to the imposed strains and therefore determine the actual behavior of the modeled media. Stresses can be interpreted as internal forces (force exerted per unit area) arising under the impact of the applied external forces. They are measured using stress tensors. Similar to the above paragraph different variants exist and are used in the given context:

Cauchy Stress Tensor. The Cauchy or true stress tensor $\varvec{\sigma }:\tilde{\varGamma }_{\mathrm{EE }}\rightarrow \mathbb R ^{3\times 3}$ is defined with respect to the deformed geometry $\tilde{\varGamma }_{\mathrm{EE }}=\varphi (\varGamma _{\mathrm{EE }})$ , i.e. it relates forces to areas in the deformed configuration. To define $\varvec{\sigma }$ , consider a small cube in the object of interest with the normals of the cube faces being oriented along the canonical unit vectors $\varvec{e}_x$ , $\varvec{e}_y$ , $\varvec{e}_z$ of $\mathbb R ^{3}$ . Let further $\varvec{T}^{(\varvec{e}_i)}$ be the stress vector acting on the face with its normal being oriented along $\varvec{e}_i$ . Then, the Cauchy stress tensor is given by

$\begin{aligned} \varvec{\sigma }= \left( \varvec{T}^{(\varvec{e}_x)}, \varvec{T}^{(\varvec{e}_y)}, \varvec{T}^{(\varvec{e}_z)}\right) ^T =\left( \begin{array}{lll} \sigma _{xx} &{} \sigma _{xy} &{} \sigma _{xz}\\ \sigma _{yx} &{} \sigma _{yy} &{} \sigma _{yz}\\ \sigma _{zx} &{} \sigma _{zy} &{} \sigma _{zz}\\ \end{array} \right) . \end{aligned}$

(4.5)
Similar to the strain tensors, the diagonal elements are referred to as normal stresses and the off-diagonal elements as shear stresses. The stress components are positive, if they act in positive direction of the coordinate axes.
First and Second Piola-Kirchhoff Stress Tensors. The first and second Piola-Kirchhoff stress tensors $\varvec{S}:\varGamma _{\mathrm{EE }}\rightarrow \mathbb R ^{3\times 3}$ and $\varvec{\varSigma }:\varGamma _{\mathrm{EE }}\rightarrow \mathbb R ^{3\times 3}$ are defined with respect to the undeformed geometry. The first Piola-Kirchhoff or Lagrangian stress tensor links forces in the deformed and areas in the undeformed configuration. Its relation to the Cauchy stress tensor is given by

$\begin{aligned} \varvec{S}=\det (\varvec{\nabla }\varphi )\varvec{\sigma }\varvec{\nabla }\varphi ^{-T} \end{aligned}$

(4.6)
with $\varvec{\nabla }\varphi =\mathbf{1 }+\varvec{\nabla }\varvec{u}$ being the Jacobian matrix of the transformation $\varphi$ .

The second Piola-Kirchhoff stress tensor relates forces and areas in the initial configuration. Its relationship to $\varvec{S}$ is given by

$\begin{aligned} \varvec{\varSigma }=\varvec{\nabla }\varphi ^{-1}\varvec{S}. \end{aligned}$

(4.7)

Note that for infinitesimal deformations, $\varvec{\sigma }$ , $\varvec{S}$ and $\varvec{\varSigma }$ are identical; the specification whether the stress is measured with respect to the undeformed or deformed geometry is only necessary in finite strain theory [8]. However, as mentioned before, the use of finite strain theory is currently the standard approach in modeling respiratory lung motion and constitutive equations are usually formulated using the second Piola-Kirchhoff stress tensor [26, 35, 55, 56], which features a computationally favorable symmetric structure.

No consensus exists, however, about the correct formulation of the constitutive equations for lung motion modeling in detail. For the majority of material property parameters it is at least impractical to experimentally measure the values, and so only a limited number of studies can be found about the measurement of mechanical properties of lung tissue [1]. And even for existing experiments and measurements, it remains widely unclear how far reported parameter values are suitable for dimensioning elastic constants during modeling. While experiments are mainly conducted in vitro, the actual in vivo behavior may, e.g., be influenced by other processes than contraction of the breathing muscles, for instance by tissue perfusion. As aggravating aspects, lung tissue properties are also affected by various factors such as age [34], lung parenchyma distortion [21], or tissue location [61]. This makes the formulation of constitutive equations difficult—especially when aiming at precise and patient-specific motion modeling. Currently, usually hyperelastic material models are chosen. Hyperelastic materials can be defined as represented by a stress-strain relation that is determined by a strain energy density function $$W$$

so that $\varvec{\varSigma }(\varvec{E})=\frac{\partial W}{\partial \varvec{E}}$ (cf. [8, 63]). Thus, it is guaranteed that the stress state depends only on the instantaneous strain state and is independent of history or rate of loading [10]. Doing so, most modeling approaches can be classified into two groups:

St. Venant-Kirchhoff material model. The St. Venant-Kirchhoff material model is the probably most common of the hyperelastic models. It holds for isotropic materials and features a linear stress-strain relationship:

$\begin{aligned} \varvec{\varSigma }\left( \varvec{E}\right) =\lambda \ {\mathrm{tr}}\left( \varvec{E}\right) \mathbf{1 }+2\mu \varvec{E}. \end{aligned}$

(4.8)
$$$\lambda >0$$” src=”/wp-content/uploads/2016/07/A217865_1_En_4_Chapter_IEq70.gif”></SPAN> and <SPAN id=IEq71 class=InlineEquation><IMG alt=$ 0$$” src=”/wp-content/uploads/2016/07/A217865_1_En_4_Chapter_IEq71.gif”> are called the Lamé constants. Other pairs of elastic constants could be used instead, e.g. Young’s modulus ² . The constants can be converted by

$\begin{aligned} Y=\frac{\mu \left( 3\lambda +2\mu \right) }{\lambda +\mu } \quad {\text{ and }} \quad \nu =\frac{\lambda }{2\left( \lambda +\mu \right) }. \end{aligned}$

(4.9)
Experimentally motivated material models. Zeng et al. conducted an experimental investigation of human lung material properties based on ex vivo lung tissue specimen (specimen tested within 48 h after death) [61]. The measurements yield a non-linear stress-strain relationship, which could be fitted by a theoretically derived strain potential per unit volume (cf. [27, 61] for details):

$\begin{aligned} W&=\frac{1}{2}c \left( \ \exp \left( a_1E_{11}^2+a_2E_{22}^2+2a_4E_{11}E_{22}\right) \right. \nonumber \\&\quad +\exp \left( a_1E_{11}^2+a_2E_{33}^2+2a_4E_{11}E_{33}\right) \nonumber \\&\quad +\left. \exp \left( a_1E_{33}^2+a_2E_{22}^2+2a_4E_{33}E_{22}\right) \right) . \end{aligned}$

(4.10)
The mean values of , and were found to be 11.8 g/cm, 0.43, 0.56, and 0.32. $E_{ij}$ refers to the components of the Green-St. Venant strain tensor.