FEM model of thorax
A voxel dataset of the chest was segmented from CT slices of a male volunteer for use in the simulation of the electrophysiological activity of the heart [12]. The CT images were acquired using a CT scanner (Toshiba Aquilion; Toshiba Medical Systems, Japan) with tube voltage 120 kVp tube current of 390 mA. The slice thickness is 0.5 mm and pixel spacing is 0.72 mm for 512 × 512 pixels. This dataset was used to perform a 3D reconstruction of the sternum, ribs, vertebrae, intercostal muscles, and diaphragm which was then rendered into the STL format. Because segmenting the intercostal muscles and diaphragm from the CT slices was difficult, 3D models of the intercostal muscles and diaphragm were constructed in the following way: the positions of intercostal muscles and diaphragm attached to the bones were first determined by referencing an anatomy textbook [13]; then the surfaces of intercostal muscles were created by connecting the muscle attachments between upper and lower rib bones; finally, regarding that diaphragm is the boundary of the chest cavity, the surface of the diaphragm was generated as to attach the lung bottom. The created models of intercostal muscles and diaphragm were confirmed by clinicians. The finite element meshes shown in Fig. 1a, b, consisting of 252,795 tetrahedral elements and 63,057 nodes, were then generated from the 3D models by using the mesh generation software ANSYS ICEM CFD (ANSYS, Inc.). By assuming the intercostal muscles and diaphragm to be incompressible material, the mixed u-p formulation was introduced with three displacement nodes and one hydrostatical pressure node.
Because the diaphragm comes into contact with and slides against the chest wall during breathing [14], contact conditions were applied between the diaphragm and chest wall. By defining the distance between the diaphragm and chest wall as gap function g and representing the contact force with η, the contact conditions were introduced as \(g \ge 0\), \(\eta \le 0\), \(g \cdot \eta = 0\) by applying Lagrange multiplier method [15]. The first condition ensures that the diaphragm cannot penetrate the chest wall. The second condition means that the contact force has to be compression force. The third condition implies that the contact force can only be generated at where the contact is occurring, i.e. when g = 0. All of the degrees of freedom at the bottom of the vertebrae were fixed as a displacement boundary condition. Abdominal pressure and pleural pressure were respectively applied to the lower surface of the diaphragm and inner surface of the thorax with linear variations from 0 to 2 kPa [7] and from −0.5 to −0.75 kPa [5] during normal quiet breathing. The fibre direction of the intercostal muscles was assigned based on human cadaver data obtained from Loring (personal communications, 2012), as shown in Fig. 1c. For the diaphragm, we referenced an anatomy textbook [13] to determine the fibre direction, as shown in Fig. 1d. We assumed that the fibre distribution had a radial pattern from the central tendon to the bottom edge.
Transversely isotropic hyperelastic material model of respiratory muscles
The intercostal muscles and diaphragm were considered to be incompressible hyperelastic materials. The function \(U_{V} = { \det }\varvec{C} - 1\) represents the incompressible constraint and is enforced by the Lagrange multiplier \(\xi\), where \({ \det }\varvec{C}\) is the determinant of the right Cauchy-Green tensor \(\varvec{C}\). To exclude the effect of the volume-changing deformation, the strain energy \(U_{R}\) was assumed to depend on the first reduced invariant of \(\varvec{C}\). This was defined as \(\overline{I}_{1}^{C} = J^{ - 2/3} {\text{tr}}\varvec{C}\), where \(J = { \det }\varvec{F}\) denotes the determinant of the deformation gradient \(\varvec{F}\). By denoting t as the external conservative force on the Neumann boundary, the total potential energy \(\Phi\) of the muscles can be formulated as a function of the displacement u and Lagrange multiplier \(\xi\), as shown in Eq. (1):
$$\Phi \left( {\varvec{u},\xi } \right) = \int_{\varOmega } {U_{R} {\text{d}}\varOmega } + \int_{\varOmega } {\xi U_{V} {\text{d}}\varOmega } - \int_{\partial \varOmega } {\varvec{t} \cdot \varvec{u}{\text{dS}}} .$$
(1)
Given the principle of stationary potential energy, the variation in \(\Phi\) is equal to zero:
$$\delta \Phi = \int_{\varOmega } {\frac{{\partial U_{R} }}{{\partial \varvec{E}}}\delta \varvec{E}{\text{d}}\varOmega } + \int_{\varOmega } {\left( {\xi \frac{{\partial U_{V} }}{{\partial \varvec{E}}}\delta \varvec{E} + \delta \xi U_{V} } \right)} {\text{d}}\varOmega - \int_{\partial \varOmega } {\varvec{t} \cdot \delta \varvec{u}{\text{dS}}} = 0,$$
(2)
where \({\mathbf{E}}\) is the Green–Lagrange strain tensor. Thus, the weak form of the governing equation is as follows:
$$\left\{ \begin{array}{l} \int_{\varOmega } {\delta \varvec{E}\left( {\frac{{\partial U_{R} }}{{\partial \varvec{E}}} + \xi \frac{{\partial U_{V} }}{{\partial \varvec{E}}}} \right){\text{d}}\varOmega } - \int_{\partial \varOmega } {\varvec{t} \cdot \delta \varvec{u}{\text{dS}}} = 0 \hfill \\ \int_{\varOmega } {\delta \xi U_{V} {\text{d}}\varOmega } = 0 \hfill \\ \end{array} \right..$$
(3)
The Lagrange multiplier corresponds to the hydrostatic pressure \(p\) by \(p = - 2\xi\) [16]. In this study, the strain energy function \(U_{R}\) proposed by Martins et al. [11] was introduced to model the skeletal muscles. This model extends the traditional Hill’s three-element muscle model to three dimensions. As shown in Fig. 2, the parallel element (PE) and series element (SE) are nonlinear springs which represent passive behaviour. The contractile element (CE) produces a contractile force when the muscle is excited.
The muscle force \(F\) can be assumed to be the sum of the forces in the PE and CE or SE:
$$F = F^{\text{PE}} + F^{\text{SE}} = F^{\text{PE}} + F^{\text{CE}} .$$
(4)
Then, the nominal stress along the fibre direction \(T\) can be written as:
$$T = T^{\text{PE}} + T^{\text{SE}} = T^{\text{PE}} + T^{\text{CE}} .$$
(5)
The stress \(T^{\text{PE}}\) is defined as follows:
$$T^{\text{PE}} = T_{0}^{\text{M}} f_{\text{PE}} \left( {\overline{\lambda }_{f} } \right)$$
(6)
$$f_{\text{PE}} \left( {\overline{\lambda }_{f} } \right) = \left\{ \begin{aligned} 2aA\left( {\overline{\lambda }_{f} - 1} \right)e^{{a\left( {\overline{\lambda }_{f} - 1} \right)^{2} }} ,\quad {\text{for }}\overline{\lambda }_{f} > 1 \hfill \\ 0 ,\quad \quad \quad \quad \quad \quad \quad \quad {\text{otherwise}} \hfill \\ \end{aligned} \right.,$$
(7)
where \(\overline{\lambda }_{f}\) is the stretch ratio in the fibre direction given by \(\overline{\lambda }_{f} = \left[ {J^{ - 2/3} \varvec{C}:\left( {\varvec{N} \otimes \varvec{N}} \right)} \right]^{1/2}\). N is the initial direction of the muscle fibre, \(T_{0}^{\text{M}}\) denotes the maximum muscle nominal stress determining the maximum muscle contraction force, a and A are material constants [11], and \(f_{\text{PE}} \left( {\overline{\lambda }_{f} } \right)\) is the passive force–length relationship of the parallel element PE [11]. The stress \(T^{\text{CE}}\) is given by
$$T^{\text{CE}} = T_{0}^{\text{M}} f_{\text{CE}} \left( {\overline{\lambda }_{f} } \right)\alpha \left( {\text{t}} \right)\gamma ,$$
(8)
where γ represents the muscle’s activation level over a range from 0 to 1, α represents the time dependence of the activation and for which the maximum value can be 1, and \(f_{\text{CE}} \left( {\overline{\lambda }_{f} } \right)\) is the active force–length relationship of the contractile element CE. In this study, this was considered to be the active force–length relationship for the muscle tissue. Figure 3 shows the experimental data for the active force–length relationship of the intercostal muscles and diaphragm as measured from an adult baboon [17].
According to the experimental data, \(f_{\text{CE}} \left( {\overline{\lambda }_{f} } \right)\) of the intercostal muscles can be formulated as
$$f_{\text{CE}}^{Int} \left( {\overline{\lambda }_{f} } \right) = \left\{ \begin{array}{ll} 0.85\overline{\lambda }_{f} - 0.48, &\quad {\text{for }}0.6 \le \overline{\lambda }_{f} < 0.7 \hfill \\ 2.43\overline{\lambda }_{f} - 1.59, &\quad {\text{for }}0.7 \le \overline{\lambda }_{f} < 0.8 \hfill \\ 3.32\overline{\lambda }_{f} - 2.30, &\quad {\text{for }}0.8 \le \overline{\lambda }_{f} < 0.9 \hfill \\ 3.10\overline{\lambda }_{f} - 2.10, &\quad {\text{for }}0.9 \le \overline{\lambda }_{f} < 1.0 \hfill \\ - 2.30\overline{\lambda }_{f} + 3.30, &\quad {\text{for 1}} .0\le \overline{\lambda }_{f} < 1.1 \hfill \\ - 3.20\overline{\lambda }_{f} + 4.29, &\quad {\text{for 1}} .1\le \overline{\lambda }_{f} < 1.2 \hfill \\ 0, &\quad{\text{otherwise}} \hfill \\ \end{array} \right..$$
(9)
Then, \(f_{\text{CE}} \left( {\overline{\lambda }_{f} } \right)\) of the diaphragm can be formulated as
$$f_{{{\text{CE}}}}^{{Dia}} \left( {\overline{\lambda }_{f} } \right) = \left\{ {\begin{array}{*{20}l} {1.80\overline{\lambda }_{f} - 0.80,} & {\quad {\text{for }}0.6 \le \overline{\lambda }_{f} < 0.7} \hfill \\ {2.37\overline{\lambda }_{f} - 1.20,} & {\quad {\text{for }}0.7 \le \overline{\lambda }_{f} < 0.8} \hfill \\ {1.61\overline{\lambda }_{f} - 0.59,} & {\quad {\text{for }}0.8 \le \overline{\lambda }_{f} < 0.9} \hfill \\ {1.42\overline{\lambda }_{f} - 0.42,} & {\quad {\text{for }}0.9 \le \overline{\lambda }_{f} < 1.0} \hfill \\ { - 0.55\overline{\lambda }_{f} + 1.55,} & {\quad {\text{for 1}}.0 \le \overline{\lambda }_{f} < 1.1} \hfill \\ { - 2.11\overline{\lambda }_{f} + 3.27,} & {\quad {\text{for }}1.1 \le \overline{\lambda }_{f} < 1.2} \hfill \\ {0,} & {\quad {\text{otherwise}}} \hfill \\ \end{array} } \right..$$
(10)
The time dependence of the activation \(\alpha \left( {\text{t}} \right)\) of the intercostal muscles and diaphragm was fitted to experimental data obtained from dogs [18], as shown in Fig. 4. Concerning the parameters \(f_{\text{CE}} \left( {\overline{\lambda }_{f} } \right)\) and \(\alpha \left( {\text{t}} \right)\), by considering that the biological characteristics are similar for mammalia and the expected difference between species can be represented by the different magnitudes, in this study, the normalized \(f_{\text{CE}} \left( {\overline{\lambda }_{f} } \right)\) and \(\alpha \left( {\text{t}} \right)\) from animal experiment shown in Figs. 3 and 4 were introduced and their magnitudes were then adjusted to be able to reproduce the behaviour of human.
The muscles were modelled as incompressible transversely isotropic hyperelastic materials for which the strain energy function can be written as follows:
$$U_{R} = U_{I} \left( {\overline{I}_{1}^{C} } \right) + U_{f} \left( {\overline{\lambda }_{f} ,\alpha ,\gamma } \right)$$
(11)
$$U_{I} = ce^{{b\left( {\overline{I}_{1}^{C} - 3} \right) - 1}}$$
(12)
$$U_{f} \left( {\overline{\lambda }_{f} ,\alpha ,\gamma } \right) = U_{\text{PE}} \left( {\overline{\lambda }_{f} } \right) + U_{\text{CE}} \left( {\overline{\lambda }_{f} ,\alpha ,\gamma } \right)$$
(13)
$$U_{\text{PE}} \left( {\overline{\lambda }_{f} } \right) = T_{0}^{\text{M}} \int_{1}^{{\overline{\lambda }_{f} }} {f_{\text{PE}} \left( \lambda \right){\text{d}}\lambda }$$
(14)
$$U_{\text{CE}} \left( {\overline{\lambda }_{f} ,\alpha ,\gamma } \right) = T_{0}^{\text{M}} \int_{1}^{{\overline{\lambda }_{f} }} {f_{\text{CE}} \left( \lambda \right)\alpha \left( {\text{t}} \right)\gamma {\text{d}}\lambda } ,$$
(15)
where \(U_{I}\) and \(U_{f}\) are the strain energies stored in the isotropic hyperelastic matrix and muscle fibre, respectively, and \(c\) and \(b\) are material constants. For an incompressible hyperelastic material, by denoting \(p\) as the hydrostatic pressure, the second Piola–Kirchhoff stress tensor S can be derived as follows:
$$\begin{aligned} \varvec{S} & = - p\varvec{C}^{ - 1} + \frac{{\partial U_{R} }}{\partial E} = - p\varvec{C}^{ - 1} + \frac{{\partial U_{I} }}{{\partial \overline{I}_{1}^{C} }}\frac{{\partial \overline{I}_{1}^{C} }}{{\partial \varvec{E}}} + \frac{{\partial U_{f} }}{{\partial \overline{\lambda }_{f} }}\frac{{\partial \overline{\lambda }_{f} }}{{\partial \varvec{E}}} \\ & = - p\varvec{C}^{ - 1} + \left\{ {bce^{{b\left( {\overline{I}_{1}^{C} - 3} \right)}} } \right\}\left( {2\varvec{I} - \frac{2}{3}\overline{I}_{1}^{C} \varvec{C}^{ - 1} } \right) \\ &\quad + \left\{ {T_{0}^{M} \left( {f_{PE} \left( {\overline{\lambda }_{f} } \right) + f_{CE} \left( {\overline{\lambda }_{f} } \right)\alpha \left( t \right)\gamma } \right)} \right\}\left( {\overline{\lambda }_{f}^{ - 1} \left( {\varvec{N} \otimes \varvec{N}} \right) - \frac{1}{3}\overline{\lambda }_{f} \varvec{C}^{ - 1} } \right). \\ \end{aligned}$$
(16)
The material parameters used for the intercostal muscles and the diaphragm were \(T_{0}^{\text{M}}\) = 0.3 MPa [17], c = 0.009348 MPa, and b = 1.4939 MPa [19]. The Lamé parameters μ and λ of the isotropic elastic model were 569.23 and 853.85 MPa, respectively for bone [6]; 65.38 and 98.08 MPa, respectively, for cartilage [6]; and 5.18 and 10.05 MPa, respectively, for tendons [20].
Simulation procedure
In this study, by omitting the inertial force, the static analyses were performed to examine and validate the introduced Hill-type transversely isotropic hyperelastic continuum skeletal muscle model for thorax deformation. The muscle contraction was activated by assigning parameter γ in Eq. (8) to be greater than zero. First, we individually activated the intercostal muscles and diaphragm to separately investigate the effect of each muscle on the chest deformation, and also to confirm the representative motions of the ribs and the diaphragm. Then, we simultaneously activated the intercostal muscles and diaphragm in order to simulate the thorax deformation. Finally, the resultant thorax deformation was compared with 4D-CT data during normal quiet breathing. The 4D-CT images of a healthy adult male were acquired using a 320-MDCT scanner (Aquilion One, Toshiba Medical Systems, Japan) for three phases during inspiration and expiration, respectively. To reduce radiation exposure, the detector row was reduced to 160, the tube voltage was reduced to 80 kVp and the tube current was reduced to 30 mA. The field of view was set to 16 cm to cover the diaphragm. The data were collected at Tokuyama Central Hospital, Japan Community Healthcare Organization, based on the Declaration of Helsinki. The procedure was approved by the Ethics Committee of the hospital.