NM, JP-P, MCe, and MG: contributed conception and design of the study; NM, JP-P, and MCo: developed the uncertainty quantification study and its implementation in the HTC environment; NM: performed the computational studies, processed the data and obtained results; NM, MCe, PM, and MG: interpreted the results and discussed the resulting conclusions; NM: wrote the first draft of the manuscript and all authors contributed to the manuscript revision, read and approved the submitted version.

#### ACKNOWLEDGMENTS

This work was partly supported by the Spanish Ministry of Economy and Competitiveness under the Maria de Maeztu Units of Excellence Program (MDM-2015-0502), by the AGAUR grant 2016-PROD-00047, the European Union Seventh Framework Program (FP7/2007-2013), Grant agreement 304857, HEAR-EU project and the QUAES Foundation Chair for Computational Technologies for Healthcare.

International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (Lake Buena Vista; Orlando, FL).

**Conflict of Interest Statement:** Author PM was employed by company MED-EL, Austria.

The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Mangado, Pons-Prats, Coma, Mistrík, Piella, Ceresa and González Ballester. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Using Non-linear Homogenization to Improve the Performance of Macroscopic Damage Models of Trabecular Bone

#### Francesc Levrero-Florencio1,2 and Pankaj Pankaj 2 \*

*1 Computational Cardiovascular Science, Department of Computer Science, University of Oxford, Oxford, United Kingdom, 2 Institute for Bioengineering, School of Engineering, The University of Edinburgh, Edinburgh, United Kingdom*

Realistic macro-level finite element simulations of the mechanical behavior of trabecular bone, a cellular anisotropic material, require a suitable constitutive model; a model that incorporates the mechanical response of bone for complex loading scenarios and includes post-elastic phenomena, such as plasticity (permanent deformations) and damage (permanent stiffness reduction), which bone is likely to experience. Some such models have been developed by conducting homogenization-based multiscale finite element simulations on bone micro-structure. While homogenization has been fairly successful in the elastic regime and, to some extent, in modeling the macroscopic plastic response, it has remained a challenge with respect to modeling damage. This study uses a homogenization scheme to upscale the damage behavior from the tissue level (microscale) to the organ level (macroscale) and assesses the suitability of different damage constitutive laws. Ten cubic specimens were each subjected to 21 strain-controlled load cases for a small range of macroscopic post-elastic strains. Isotropic and anisotropic criteria were considered, density and fabric relationships were used in the formulation of the damage law, and a combined isotropic/anisotropic law with tension/compression asymmetry was formulated, based on the homogenized results, as a possible alternative to the currently used single scalar damage criterion. This computational study enhances the current knowledge on the macroscopic damage behavior of trabecular bone. By developing relationships of damage progression with bone's micro-architectural indices (density and fabric) the study also provides an aid for the creation of more precise macroscale continuum models, which are likely to improve clinical predictions.

Keywords: trabecular bone, multiscale modeling, parameter estimation, continuum damage, finite element method, homogenization, biomechanics, high performance computing

### 1. INTRODUCTION

The growth of older population around the world in the last few decades has caused an increase in problems which can be associated to deteriorated mechanical properties of bone; osteoporosis is the clearest example of one such condition.

Computer models have been extensively employed to evaluate the mechanical response of bone and bone-implant systems under a range of loading scenarios (Pankaj, 2013). Previous

#### Edited by:

*Alfons Hoekstra, University of Amsterdam, Netherlands*

#### Reviewed by:

*Bradley John Roth, Oakland University, United States Kumari Sonal Choudhary, University of California, San Diego, United States*

> \*Correspondence: *Pankaj Pankaj<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology*

Received: *19 January 2018* Accepted: *27 April 2018* Published: *17 May 2018*

#### Citation:

*Levrero-Florencio F and Pankaj P (2018) Using Non-linear Homogenization to Improve the Performance of Macroscopic Damage Models of Trabecular Bone. Front. Physiol. 9:545. doi: 10.3389/fphys.2018.00545*

**210**

studies have assumed bone to be homogeneous (Completo et al., 2009; Conlisk et al., 2015), i.e., its properties do not vary from point to point in space or heterogeneous (Helgason et al., 2008; Schileo et al., 2008; Tassani et al., 2011), i.e., its properties vary with location (these are typically assigned on the basis of greyscale values observed in micro-computed tomography scans). However, in the large majority of studies, bone is assumed to be linear elastic and isotropic, i.e. its properties at a certain point in space are the same in all directions. It is well- recognized that the cellular microstructure of trabecular bone renders it anisotropic (Turner et al., 1990; Odgaard et al., 1997), i.e., properties at a point in space vary in different directions. Finite element (FE) analysis of the bone microstructure, in which the solid and pore phases are explicitly modeled, has been used to evaluate the homogenized anisotropic linear elastic properties of bone in the past two decades. Morphology-elasticity relationships that use bone density and fabric have also been established, with fabric typically measured through the mean intercept length (MIL) fabric tensor (Harrigan and Mann, 1984). These relationships establish links between density, fabric, and the components of the stiffness tensor (Zysset, 2003). More recently, some studies have attempted the evaluation of homogenized yield behaviour (Cowin, 1986; Wolfram et al., 2012; Levrero-Florencio et al., 2016).

Homogenized FE models of the whole bone can include microstructural information at the continuum (macroscopic) level and can thus improve the assessment of the behavior of bone and bone-implant systems in clinical scenarios. Homogenization relies on averaging the strains and stresses over a representative volume element (RVE) of the considered material; it is the most widely used multiscale approach to study the macroscopic behavior of trabecular bone. Homogenization of an RVE in the post-elastic regime requires examining its response to a wide range of loading scenarios (Bayraktar et al., 2004; Levrero-Florencio et al., 2016, 2017a). It is important to note that, in experiments, it is not possible to test multiple load cases after a certain load threshold has been surpassed because permanent deformations and/or damage caused during the first loading case will affect the behavior in subsequent loading cases. Therefore, computational means provide an attractive alternative. Nonetheless, the need for fine resolution to recreate a biofidelic geometry of the bone microstructure leads to micro-FE (µFE) systems of several tens of millions of degrees of freedom. The need to undertake multiple load cases each in non-linear regime requires the usage of high performance computing (HPC) platforms and software which can take advantage of them.

Although the damage behavior of bone has been considered in a few studies (Keaveny et al., 1994a; Garcia et al., 2009; Shi et al., 2010; Schwiedrzik and Zysset, 2013; Lambers et al., 2014), there are apparent limitations to most of the employed mathematical formulations. For example, most macroscopic damage models of trabecular bone employ an isotropic damage evolution, i.e., a "basic," or single scalar isotropic formulation, as mentioned in Carol et al. (2002), and do not take into account that the development of damage may be related to the load case being considered (Levrero-Florencio et al., 2017a). The authors have previously conducted a series of uniaxial simulations which show that damage develops differently in tension−compression, and in normal−shear (Levrero-Florencio et al., 2017a).

This study has a number of aims. Firstly, it extends the study performed in Levrero-Florencio et al. (2017a) by adding 12 biaxial macroscopic cases in the normal strain space. The second aim is to examine the suitability of certain damage mechanisms by fitting different damage laws to the damaged macroscopic stiffness tensors. The study then investigates the possible relationships between the macroscopic damage behavior of trabecular bone and its density and fabric description, by including these micro-architectural indices as additional data in the fitting procedure. The data for these formulations is obtained computationally through homogenization-based multiscale simulations run on a HPC platform with an in-house developed parallel implicit FE code.

### 2. NOTATION

The mathematical operators defined in this section largely follow the notation used in Wu and Li (2008), Schwiedrzik et al. (2013), and Levrero-Florencio et al. (2016). Compact tensor notation is used throughout this study, with indicial notation within brackets being used in this section to clarify certain tensorial operations, or in specific sections where further clarification might be required.

As a general rule, scalars are denoted with Greek or Latin italic characters (e.g., λ or a, respectively); vectors, or first-order tensors, are denoted by Latin bold lower-case characters (e.g., **a**); second-order tensors are denoted with Greek or Latin bold upper-case characters (e.g., σ or **A**, respectively); and fourthorder tensors are denoted by Latin double-barred upper-case characters (e.g., A).

Tensorial operations are denoted as follows. Single contraction of tensorial entities may appear as **a** · **b** (aibi), **a** · **B** (aiBij), **Ab** (Aijbj), or **AB** (AikBkj), note that the scalar product symbol (·) only appears when the first entity to be contracted is a first-order tensor; double contraction of tensorial entities may appear as **A** : **B** (AijBij), A : **B** (AijklBkl), **A** :B (AijBijkl), or A : B (AijmnBmnkl). Different tensor products have been defined, which include **a** ⊗ **b** (aibj), **A** ⊗ **B** (AijBkl), **A**⊗**B** (AikBjl ), **A**⊗**B** (AilBjk ), or **A**⊗**B** = 1 2 (**A**⊗**B** + **A**⊗**B**) ( 1 2 [AikBjl + AilBjk]).

Curly brackets {·} are used to represent vector projections of second-order tensors, such as

$$\begin{Bmatrix} \mathbf{A} \end{Bmatrix} = \begin{Bmatrix} A\_{11} & A\_{22} & A\_{33} & A\_{12} & A\_{13} & A\_{23} \end{Bmatrix}^{\mathrm{T}}.\tag{1}$$

Square brackets [·] are used, in conjunction with parentheses (·), to indicate priority in the order of mathematical operations; an important exception occurs when square brackets are used to represent the matrix projection of a fourth-order tensor, such as

$$[\mathbb{A}] = \begin{bmatrix} A\_{1111} & A\_{1122} & A\_{1133} & A\_{1112} & A\_{1113} & A\_{1123} \\ A\_{2211} & A\_{2222} & A\_{2233} & A\_{2212} & A\_{2213} & A\_{2223} \\ A\_{3311} & A\_{3322} & A\_{3333} & A\_{3312} & A\_{3313} & A\_{3323} \\ A\_{1211} & A\_{1222} & A\_{1233} & A\_{1212} & A\_{1213} & A\_{1223} \\ A\_{1311} & A\_{1322} & A\_{1333} & A\_{1312} & A\_{1313} & A\_{1323} \\ A\_{2311} & A\_{2322} & A\_{2333} & A\_{2312} & A\_{2313} & A\_{2323} \end{bmatrix} . \tag{2}$$

Double vertical bars k(·)k are used to represent the Frobenius norm of the matrix (·), such as the Frobenius norm of the following 3 × 3 symmetric matrix,

$$\|\mathbf{[A]}\| = \sqrt{A\_{11}^2 + A\_{22}^2 + A\_{33}^2 + 2A\_{11}^2 + 2A\_{13}^2 + 2A\_{23}^2}.\tag{3}$$

### 3. MATERIALS AND METHODS

#### 3.1. Computational Methods

This section follows the "Materials and Methods" section in Levrero-Florencio et al. (2017a). The authors used µCT images of trabecular bone samples to create detailed FE models, which ranged from 10 to 30 million elements, representing the solid phase of bone for a cubic trabecular bone samples (which includes both solid phase and pores) of size 5 mm. In the study conducted by Levrero-Florencio et al. (2017a), plasticity and damage were considered for the solid phase post-elastic properties and nine uniaxial strain cases were investigated (load cases 1 to 9 of **Table 1**) representing: three tensile cases (+ε11, +ε22, and +ε33), three compressive cases (−ε11, −ε22, and −ε33), and three shear cases (ε12, ε13, and ε23). The macroscopic damage behavior was studied by using an appropriate homogenizationbased multiscale technique, which is explained later.

Trabecular bone is an anisotropic material; its anisotropy may be quantified with a fabric tensor, which indicates how directionally distributed a material is. The Mean Intercept Length (MIL) fabric tensor is used in this study because it is widely used in trabecular bone studies, and it performs slightly better than other fabric measures (Kabel et al., 1999; Zysset, 2003). The magnitude of an eigenvalue of the MIL fabric tensor denotes the proportion of material which is aligned in the direction expressed in the correspondent eigenvector. The fabric tensors are normalized by a trace equal to 3 (Zysset, 2003).

In this study, 10 out of the 12 samples employed in Levrero-Florencio et al. (2017a) were subjected to 12 additional biaxial strain cases in the normal strain space (**Table 1**, cases 10–21). Kinematic uniform boundary conditions (i.e., conditions in which displacements, or macroscopic strains, are controlled) were used for all analyses; these are known for providing an upper bound for the macroscopic stiffness tensor and macroscopic yield surface of trabecular bone (Wang et al., 2009; Panyasantisuk et al., 2015). An example of how boundary conditions are implemented can be seen in **Figure 1**, which corresponds to load case 4 in **Table 1**. The morphological indices of these samples are shown in **Table 2**. BV/TV stands for bone volume over total volume and it is a surrogate for density, DOA stands for degree of anisotropy and it is the ratio of the highest to the lowest eigenvalues of the MIL fabric tensor, and SMI stands for structure model index and it ranges from rod- (SMI = 3) to plate-shaped (SMI = 0) microstructure.

The 10 samples were aligned with the MIL fabric tensor eigenvectors, with the eigenvalues sorted in descendent order (m1 > m2 > m3). The samples were then meshed with trilinear hexahedra and subjected to the aforementioned straincontrolled load cases; the largest mesh consisted of ∼27M degrees of freedom, leading to square sparse stiffness matrices of up to 27M×27M elements. The considered constitutive law at the

tissue level was isotropic with coupled plasticity and damage (the former captures irrecoverable deformations while the latter takes accounts for stiffness reduction), meaning that damage and plasticity interact with each other and evolve at the same time; the considered yield surface was Drucker-Prager (Tai et al., 2006; Carnelli et al., 2010; Panyasantisuk et al., 2015) with yield values corresponding to 0.41% strain in tension and 0.83% strain in

image of one of the used trabecular bone specimens; this particular sample led to a FE mesh of ∼21M degrees of freedom.

TABLE 2 | Morphological indices of the 10 used specimens.

compression (Bayraktar and Keaveny, 2004). Linear isotropic hardening corresponding to 5% of the undamaged elastic slope (Wolfram et al., 2012; Sanyal et al., 2015) was used. At the tissue level, damage evolution was assumed to be isotropic and it was obtained from Schwiedrzik and Zysset (2013, 2015). The maximum damage was capped at 0.9 (90% isotropic stiffness reduction) to avoid numerical difficulties related to the loss of positive-definiteness of the stiffness matrix; this was performed by using

$$D(\varepsilon^{\mathcal{P}}) = D\_{\max} \left( 1 - e^{-k\_{\mathcal{P}} \varepsilon^{\mathcal{P}}} \right) \tag{4}$$

where ε p = kε pk is the accumulated plastic strain, Dmax is the maximum damage, and kp is a parameter obtained from Schwiedrzik and Zysset (2015).

The µFE simulations were run on a Cray XC30 supercomputer hosted by ARCHER (UK National Supercomputing Service), with an in-house version of ParaFEM (Smith et al., 2013; Levrero-Florencio et al., 2017b) which solves implicit quasi-static finite strain elastoplasticity problems in a highly scalable message passing interface-based (MPI) parallel fashion. Each simulation took from 40 to 120 min when using 1,920 cores, depending on the considered load case, with biaxial compression-compression load cases taking the longest. In order to improve the convergence aspect of the local (constitutive level, i.e., at each integration point) Newton–closest-point projection method (Newton-CPPM), two additional schemes were implemented: (a) a line search as in the primal-CPPM scheme described in Pérez-Foguet and Armero (2002) and (b) an improved trial predictor (Bicani ´ c´ and Pearce, 1996; de Souza Neto et al., 2008). In the latter scheme, if the first Newton-CPPM fails to converge, it is restarted but this time with the initial guess for stress as σ proj, which is the stress returned to the frozen yield surface, i.e., no hardening or damage evolution. If these two mechanisms do not work, to ensure that a possible local lack of convergence does not influence the results of the µFE simulations, lack of convergence of the CPPM scheme is broadcasted to all MPI processes in order to cut down the time increment to half of its value. The initial, and maximum, step size corresponded to 0.1% macroscopic strain Frobenius norm and was allowed to decrease to a minimum of 0.001%, if global (structural level, i.e., the global stiffness matrix) or local convergence was not achieved. The global solution scheme employed was Newton-Raphson, and a Jacobi, or diagonally, preconditioned conjugate gradients method was used as the linear algebraic solver.

The macroscopic elastic stiffness tensor was calculated at each time increment by using the homogenization procedure described by van Rietbergen et al. (1995, 1996), in which the macroscopic elastic stiffness tensor E is

$$\mathbb{E} = \frac{1}{V} \int\_{\Omega} (1 - D\_{\mu}) \mathbb{E}\_{\mu} : \mathbb{M} \, \mathrm{d}V,\tag{5}$$

which, in a FE setting, is equivalent to

$$\mathcal{E} = \frac{1}{V} \sum\_{i=1}^{n\_{\text{cls}}} \sum\_{j=1}^{n\_{\text{ips}}} (1 - D\_{\mu \text{ ij}}) \mathbb{E}\_{\mu \text{ ij}} : \mathbb{M}\_{ij} \det(\mathbf{J}\_{ij}) \mathbf{w}\_j,\tag{6}$$

and where V is the volume of the cubic region (5×5×5 = 125 mm3 ), Dµ is the damage at the solid phase, Eµ is the solid phase undamaged stiffness tensor, nels is the total number of elements in the considered mesh, nips is the number of integration points in a trilinear hexahedron, det(**J**ij) is the determinant of the Jacobian of the transformation from normal to natural coordinates, wj is the weight of the corresponding integration point, and M is the local structure tensor, which relates the solid phase strain εµ to the average strain tensor ε, such that

$$\mathfrak{s}\_{\mu} = \mathbb{M} : \mathfrak{s}.\tag{7}$$

This tensor M was determined by solving six completely linear FE systems for six macroscopic uniaxial strain cases (three tensile or compressive and three shear). For each of these cases, the tissue strains calculated represent one of the six columns of the matrix projection of M (Hollister and Kikuchi, 1992). The assumption made was that the samples are aligned in their orthotropic axes as they were aligned with the MIL fabric tensor eigenvectors (Odgaard et al., 1997). Macroscopic strain points were defined by using the 0.2% strain criterion (Wolfram et al., 2012; Levrero-Florencio et al., 2017a), and it was extended to define further 0.3, 0.4, and 0.5% strain levels. The corresponding damaged slope to calculate these strain points is determined at each time step, depending on the load case. The following is an example for the biaxial tensile case ε11 = ε22 > 0 (load case 10 in **Table 1**). Since the macroscopic strains are small, the assumption of linear kinematics can be considered at the macroscale; thus, the homogenized infinitesimal stress can be obtained through the macroscopic infinitesimal strain and the macroscopic stiffness tensor, such as

$$
\begin{aligned}
\begin{Bmatrix}\sigma\_{\text{hom},11} \\ \sigma\_{\text{hom},21} \\ \sigma\_{\text{hom},3} \\ \sigma\_{\text{hom},12} \\ \sigma\_{\text{hom},13} \\ \sigma\_{\text{hom},23} \\ \sigma\_{\text{hom},23} \end{Bmatrix} &= \begin{Bmatrix}E\_{1111} & E\_{1122} & E\_{1133} & E\_{1112} & E\_{1113} & E\_{1123} \\ E\_{2211} & E\_{222} & E\_{2233} & E\_{2212} & E\_{2312} & E\_{2313} \\ E\_{3111} & E\_{3222} & E\_{1233} & E\_{1212} & E\_{1213} & E\_{1223} \\ E\_{1311} & E\_{3122} & E\_{1333} & E\_{1312} & E\_{1313} & E\_{1323} \\ E\_{2311} & E\_{2322} & E\_{2333} & E\_{2312} & E\_{2313} & E\_{2313} \end{Bmatrix} \begin{Bmatrix} \varepsilon\_{\text{l}} \\ \varepsilon\_{\text{l}} \\ 0 \\ 0 \\ 0 \\ 0 \end{Bmatrix}, \\ &= \begin{Bmatrix} E\_{1111}\varepsilon\_{\text{l}} + E\_{1122}\varepsilon\_{\text{l}} \\ 0 \\ 0 \\ 0 \\ 0 \end{Bmatrix}, \end{aligned} \tag{8}$$

where σ hom is the homogenized stress tensor, leading to

$$\|\sigma\_{\text{hom}}\| = \sqrt{E\_{1111}^2 \varepsilon\_{11}^2 + E\_{1122}^2 \varepsilon\_{22}^2 + E\_{2211}^2 \varepsilon\_{11}^2 + E\_{2222}^2 \varepsilon\_{22}^2}, \quad \text{(9)}$$

with the damaged slope being (note that in the considered biaxial cases |εii| = εjj )

$$K\_{\rm dam} = \sqrt{E\_{1111}^2 + E\_{1122}^2 + E\_{2211}^2 + E\_{2222}^2}.\tag{10}$$

#### 3.2. Theoretical Framework of Damage

The previously described µFE simulations, together with the homogenization-based multiscale procedure, were used to derive the damaged macroscopic stiffness tensors of the considered samples, for different load scenarios (**Table 1**) and load levels (0.2, 0.3, 0.4, and 0.5% strain norm). These stiffness tensors were used as data points for a minimization procedure (described in the following subsections), which was used to fit the macroscopic damage behavior to several theoretical damage models: single scalar isotropic formulation, three scalars anisotropic formulation, and isotropic/anisotropic combined formulation with tension/compression asymmetry.

Coupled damage and plasticity were considered for the µFE simulations. However, the focus of this study is on the macroscopic damage behavior of trabecular bone and therefore no plasticity is assumed at the macroscale. This is why, in the following, the total strain ε is used instead of the elastic strain ε e .

#### 3.2.1. Basic Concepts and Description of the Baseline Model

Let us consider the theoretical framework of elastic degradation by using state variables, from which the different damage constitutive models are derived (Carol et al., 1994, 2002; Murakami, 2012). The starting point of the theoretical framework is the assumption of a Helmholtz free energy potential per unit reference volume ψ of the considered material, from which the state equations are derived. The free energy potential may be expressed as

$$\begin{aligned} \psi(\mathfrak{s}, D\_k, R\_k) &= \psi^\mathfrak{e}(\mathfrak{s}, D\_k) + \psi^D(R\_k) \\ &= \frac{1}{2} \mathfrak{s} : \mathbb{E}(\mathbb{E}\_0, D\_k) : \mathfrak{s} + \frac{1}{2} \sum\_{k=1}^l K\_k R\_k^2, \end{aligned} \tag{11}$$

where ε is the infinitesimal strain tensor, E and E0 are, respectively, the damaged and undamaged stiffness tensors, Dk are a set of l scalar damage variables; Rk and Kk are, respectively, a set of l variables and l parameters controlling the size and hardening of the (damage) dissipation potential functions Fk (Equation 16).

Time derivative of Equation (11) yields

$$\boldsymbol{\psi} = \frac{\partial \boldsymbol{\psi}}{\partial \mathbf{s}} : \dot{\mathbf{s}} + \sum\_{k=1}^{l} \frac{\partial \boldsymbol{\psi}}{\partial D\_{k}} \dot{\mathbf{D}}\_{k} + \sum\_{k=1}^{l} \frac{\partial \boldsymbol{\psi}}{\partial R\_{k}} \dot{\mathbf{R}}\_{k},\tag{12}$$

which, when used in the Clausius-Duhem inequality for isothermal processes

$$\boldsymbol{\sigma} : \dot{\mathbf{s}} - \boldsymbol{\rho} \,\, \dot{\boldsymbol{\psi}} \succeq \mathbf{0},\tag{13}$$

gives rise to the dissipation inequality

$$\phi = \left(\sigma - \rho \frac{\partial \psi}{\partial \mathbf{z}}\right) \cdot \dot{\mathbf{s}} - \sum\_{k=1}^{l} \rho \frac{\partial \psi}{\partial D\_k} \dot{D}\_k - \sum\_{k=1}^{l} \rho \frac{\partial \psi}{\partial R\_k} \dot{R}\_k$$

$$= \sum\_{k=1}^{l} Y\_k \dot{D}\_k + \sum\_{k=1}^{l} B\_k \dot{R}\_k \ge 0,\tag{14}$$

where ρ is the density of the considered material, σ = ρ ∂ψ ∂ε , Yk = −ρ ∂ψ ∂Dk , and Bk = −ρ ∂ψ ∂Rk = KkRk .

The evolution equations of Dk and Rk are derived from the corresponding dissipation potential functions Fk , leading to

$$
\dot{D}\_k = \dot{\eta}\_k \frac{\partial F\_k}{\partial Y\_k}; \quad \dot{R}\_k = \dot{\eta}\_k \frac{\partial F\_k}{\partial B\_k}, \tag{15}
$$

where γ˙k are indeterminate multipliers. Since Fk also delimit the undamaged region of the considered material, the nonnegativeness of Equation (14) is assured (Murakami, 2012). Linear, a priori uncoupled, criteria for Fk are considered in this study (each Dk is related to a single Fk ), such that

$$F\_k(Y\_k, B\_k) = Y\_k - (B\_k + B\_{k,0}) = Y\_k - (K\_k R\_k + B\_{k,0}) \le 0,\text{ (16)}$$

where Bk,0 are the initial sizes of Fk , i.e., when Rk = 0. These linear functions are considered for the sake of simplicity and also because data on additional strain points is needed so that more complex, non-linear, evolution expressions of the dissipation potentials may be taken into account.

Energy equivalence is adopted here since it automatically induces major symmetry in the stiffness and compliance tensors. This leads to

$$\begin{split} \psi^{\varepsilon}(\mathfrak{s}, D\_{k}) &= \frac{1}{2} \mathfrak{s} : \mathbb{E}(\mathbb{E}\_{0}, D\_{k}) : \mathfrak{s} = \frac{1}{2} \mathfrak{s} : \mathbb{M}^{\mathrm{T}}(D\_{k}) : \mathbb{E}\_{0} : \mathbb{M}(D\_{k}) : \mathfrak{s} \\ &= \frac{1}{2} \mathfrak{s}\_{\mathrm{eff}}(\mathfrak{s}, D\_{k}) : \mathbb{E}\_{0} : \mathfrak{s}\_{\mathrm{eff}}(\mathfrak{s}, D\_{k}), \end{split} \tag{17}$$

where M is the fourth-order damage effect tensor which depends on the considered damage formulation, and A T is defined so that A T ≡ A T ijkl = Aklij.

#### 3.2.2. Numerical Solution of the Damage Models

Equations (15, 16) are integrated with Backward Euler. Residual equations for each of the variables to be sought can be formulated, with a format similar to that of CPPM equations of computational plasticity (Armero and Pérez-Foguet, 2002; Pérez-Foguet and Armero, 2002), so that

$$\begin{Bmatrix} R\_{D,k} \\ R\_{R,k} \\ F\_k \end{Bmatrix} = \begin{Bmatrix} D\_{k,n+1} - D\_{k,n} - \Delta \gamma\_{k,n+1} \frac{\partial F\_k}{\partial Y\_k} \Big|\_{n+1} \\ R\_{k,n+1} - R\_{k,n} - \Delta \gamma\_{k,n+1} \frac{\partial F\_k}{\partial B\_k} \Big|\_{n+1} \\ Y\_{k,n+1} - (K\_k R\_{k,n+1} + B\_{k,0}) \end{Bmatrix} \tag{18}$$

where n stands for the nth time increment, and the vertical bar means "evaluated at".

The resulting set of non-linear equations (Equation 18) can be solved with a numerical scheme, for instance a Newton-Raphson approach. The first step is to calculate the Jacobian of the system, and therefore the residuals (Equation 18) are linearized, leading to (time subscripts are dropped for convenience from now onwards)

$$\begin{Bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{0} \end{Bmatrix} = \begin{Bmatrix} \mathbf{d}D\_{j} \Big( \delta\_{jk} - \Delta \boldsymbol{\eta}\_{k} \frac{\partial}{\partial D\_{j}} \frac{\partial F\_{k}}{\partial Y\_{k}} \Big) - \mathbf{d}\Delta \boldsymbol{\eta}\_{k} \frac{\partial F\_{k}}{\partial Y\_{k}} \\ \mathbf{d}R\_{k} - \mathbf{d}\Delta \boldsymbol{\eta}\_{k} \frac{\partial F\_{k}}{\partial B\_{k}} \\ \mathbf{d}D\_{j} \frac{\partial F\_{k}}{\partial D\_{j}} + \mathbf{d}R\_{k} \frac{\partial F\_{k}}{\partial R\_{k}} \end{Bmatrix}.\tag{19}$$

where δij = ( 0 if i 6= j 1 if i = j is the Kronecker delta. The specific

expressions for the derivatives of the Jacobian are presented for each of the considered damage models in the following sections.

The resulting Newton-Raphson scheme to solve for Dk , Rk , and 1γk is

$$\begin{Bmatrix} D\_k \\ R\_k \\ \Delta \boldsymbol{\chi}\_k \end{Bmatrix}\_{m+1} = \begin{Bmatrix} D\_k \\ R\_k \\ \Delta \boldsymbol{\chi}\_k \end{Bmatrix}\_m - \begin{bmatrix} \delta\_{jk} - \Delta \boldsymbol{\chi}\_k \frac{\partial}{\partial D\_j} \frac{\partial F\_k}{\partial \boldsymbol{Y}\_k} & 0 & -\frac{\partial F\_k}{\partial \boldsymbol{Y}\_k} \\ 0 & 1 & -\frac{\partial F\_k}{\partial B\_k} \\ \frac{\partial F\_k}{\partial D\_j} & \frac{\partial F\_k}{\partial R\_k} & 0 \end{Bmatrix}^{-1} \begin{Bmatrix} R\_{D,k} \\ R\_{R,k} \\ F\_k \end{Bmatrix}\_m \tag{20}$$

where m stands for the mth iteration of the Newton-Raphson scheme.

#### 3.2.3. Damage Models

This section describes the three main models, and their variants, used in this study. The first two models, single scalar isotropic model (section <IP_ADDRESS>) and three scalars anisotropic model (section <IP_ADDRESS>) are mainly used to assess the BV/TV and fabric eigenvalue dependencies of macroscopic damage models of trabecular bone. The proposed model (section <IP_ADDRESS>), we believe, is a considerable improvement upon the usually employed single scalar isotropic formulation.

#### **<IP_ADDRESS>. Single scalar isotropic formulation**

In this simple damage formulation a single scalar damage variable D equally affects all the components of the stiffness tensor, i.e., all directions are equally affected by damage. The damage effect tensor is

$$\mathbb{M} = (1 - D)\mathbb{I}\_{\text{sym}},\tag{21}$$

where Isym = **I** ⊗ **I**.

The Helmholtz free energy potential for this model is

$$\psi(\mathfrak{s}, D, R) = \frac{1}{2}\mathfrak{s} : (1 - D)^2 \mathbb{E}\_0 : \mathfrak{s} + \frac{1}{2}K\mathbb{R}^2,\tag{22}$$

which leads to the following expressions for the conjugate thermodynamic associated variables

$$\begin{aligned} \mathfrak{\mathfrak{a}} &= (1 - D)^2 \mathbb{E}\_0 : \mathfrak{a} \\ Y &= -\frac{1}{2} \mathfrak{a} : \frac{\partial \mathbb{E}}{\partial D} : \mathfrak{a} = \mathfrak{s} : (1 - D) \mathbb{E}\_0 : \mathfrak{s} \\ B &= KR \end{aligned} \tag{23}$$

and to the following expressions for the derivatives in Equation (20)

$$\begin{aligned} \frac{\partial}{\partial D} \frac{\partial F}{\partial Y} &= 0\\ \frac{\partial F}{\partial Y} &= 1\\ \frac{\partial F}{\partial D} = \frac{\partial Y}{\partial D} &= -\frac{1}{2}\mathbf{s} : \frac{\partial^2 \mathbb{E}}{\partial D^2} : \mathbf{s} = -\mathbf{s} : \mathbb{E}\_0 : \mathbf{s} \\ \frac{\partial F}{\partial R} &= -K. \end{aligned} \tag{24}$$

BV/TV dependence is included in this model by defining K = K0,iso ρ o and B = B0,iso ρ p , where ρ is the BV/TV of the considered sample, and o and p are the exponents expressing the BV/TV dependency.

#### **<IP_ADDRESS>. Three scalars anisotropic formulation**

In the anisotropic damage formulation a damage scalar for each principal direction of the sample is considered (D1, D2, and D3), meaning that each of these three orthogonal directions has a different damage behavior (as previously stated, these orthogonal directions are parallel to the axes of the cubic sample). Since the range of post-elastic strains applied to the sample is relatively small, it is assumed that no rotation of the orthotropic axes occurs. The damage effect tensor is

$$\mathbb{M} = (\mathbb{I}\_{\text{sym}} - \mathbb{D}),\tag{25}$$

where

∂D ∂D1 = 1 α α 0 0 0 α 0 0 0 0 0 α 0 0 0 0 0 0 0 0 β 0 0 0 0 0 0 β 0 0 0 0 0 0 0 ; ∂D ∂D2 = 0 α 0 0 0 0 α 1 α 0 0 0 0 α 0 0 0 0 0 0 0 β 0 0 0 0 0 0 0 0 0 0 0 0 0 β ; ∂D ∂D3 = 0 0 α 0 0 0 0 0 α 0 0 0 α α 1 0 0 0 0 0 0 0 0 0 0 0 0 0 β 0 0 0 0 0 0 β . (26)

in which α and β are parameters which determine how the components of the stiffness tensor are affected by the different damage scalars.

The Helmholtz free energy potential is

$$\psi(\mathbf{s}, D\_k, R\_k) = \frac{1}{2}\mathbf{s} : \mathbb{E}(\mathbb{E}\_0, D\_k) : \mathbf{s} + \frac{1}{2} \sum\_{k=1}^3 K\_k R\_k^2,\tag{27}$$

which leads to the following expressions for the conjugate thermodynamic associated variables

$$\begin{aligned} \boldsymbol{\sigma} &= \mathbb{E} : \mathbf{s} = [(\mathbb{I}\_{\text{sym}} - \mathbb{D}) : \mathbb{E}\_0 : (\mathbb{I}\_{\text{sym}} - \mathbb{D})] : \mathbf{s} \\ Y\_k &= -\frac{1}{2} \mathbf{s} : \frac{\partial \mathbb{E}}{\partial D\_k} : \mathbf{s} \\ B\_k &= K\_k R\_k \end{aligned} \tag{28}$$

and to the following expressions for the derivatives in Equation (20)

$$\begin{aligned} \frac{\partial}{\partial D\_j} \frac{\partial F\_k}{\partial Y\_k} &= 0\\ \frac{\partial F\_k}{\partial Y\_k} &= 1 \end{aligned}$$

$$\begin{aligned} \frac{\partial F\_k}{\partial D\_j} &= \frac{\partial Y\_k}{\partial D\_j} = -\frac{1}{2} \mathfrak{s} : \frac{\partial}{\partial D\_j} \frac{\partial \mathbb{E}}{\partial D\_k} : \mathfrak{e} \\\ \frac{\partial F\_k}{\partial R\_k} &= -K\_k \\\ \frac{\partial \mathbb{E}}{\partial D\_k} &= -\left[\frac{\partial \mathbb{D}}{\partial D\_k} : \mathbb{E}\_{\mathbb{D}} : \mathbb{M} + \mathbb{M} : \mathbb{E}\_{\mathbb{D}} : \frac{\partial \mathbb{D}}{\partial D\_k}\right] \\\ \frac{\partial}{\partial D\_j} \frac{\partial \mathbb{E}}{\partial D\_k} &= \frac{\partial \mathbb{D}}{\partial D\_k} : \mathbb{E}\_{\mathbb{D}} : \frac{\partial \mathbb{D}}{\partial D\_j} + \frac{\partial \mathbb{D}}{\partial D\_j} : \mathbb{E}\_{\mathbb{D}} : \frac{\partial \mathbb{D}}{\partial D\_k} \end{aligned}$$

Fabric eigenvalue dependencies are included in this model by defining Kk = K0,aniso m q k and Bk = B0,aniso mr k , where mk is the MIL fabric eigenvalue corresponding to the kth orthotropic direction of the sample; and q and r are the exponents expressing the fabric eigenvalue dependency.

#### **<IP_ADDRESS>. Combined formulation with tension/compression asymmetry**

We propose a combined isotropic/anisotropic damage formulation, which consists of four damage scalars: a single scalar defines the isotropic part of the model (Diso); and three scalars define the anisotropic part of the model, one for each of the three orthotropic directions (D1, D2, and D3). As in the previous cases, the isotropic damage scalar equally affects all directions, while each of the three orthotropic damage scalars only affect their corresponding orthogonal direction. It is assumed that there is no rotation of the orthotropic axes. The tension/compression asymmetry is included in the damage effect tensor, such that

$$\mathbf{M} = \mathbb{I}\_{\text{sym}} - \mathbb{D}\_{\text{iso}} - \sum\_{i=1}^{3} [1 + \eta \,\mathrm{H}(-\mathbf{m}\_{i} \cdot \mathbf{s} \,\mathrm{m}\_{i}) \mathbb{D}\_{\text{aniso},i}],\tag{30}$$

where

$$\mathbb{D}\_{\text{iso}} = (1 - D)\mathbb{I}\_{\text{sym}},\tag{31}$$

$$\mathbb{D}\_{\text{aniso},i} = \frac{\partial \mathbb{D}}{\partial D\_i} D\_i \tag{32}$$

with ∂D ∂Di being defined in Equation (26), η is the parameter governing the tension/compression asymmetry, **m**i is the i th fabric tensor eigenvector, and H(·) is the Heaviside function defined as

$$\mathbf{H}(\cdot) = \begin{cases} 1 & \text{if } (\cdot) > 0 \\ 0 & \text{if } (\cdot) \le 0 \end{cases}. \tag{33}$$

The Helmholtz free energy potential for this model is

$$\psi(\mathfrak{s}, D\_k, R\_k) = \frac{1}{2}\mathfrak{s} : \mathbb{E}(\mathbb{E}\_0, D\_k) : \mathfrak{s} + \frac{1}{2} \sum\_{k=1}^4 K\_k R\_k^2,\tag{34}$$

BV/TV and fabric eigenvalue dependencies are included in this model by defining Kiso = K0,iso ρ o ; Kk,aniso = K0,aniso ρ tm q k , k ∈ {1, 2, 3}; Biso = B0,iso ρ p ; and Bk,aniso = B0,aniso ρ umr k , k ∈ {1, 2, 3}, where o and p are the exponents expressing BV/TV dependency of the isotropic part of the model; and t, u, q, and r are, respectively the exponents expressing BV/TV and fabric eigenvalue dependencies of the anisotropic part of the model. The rest of expressions in the model are the same to those in section <IP_ADDRESS>.

#### 3.3. Fitting of the Different Damage Laws

The different damage constitutive models described in the previous section are fitted to the macroscopic damage response obtained from the homogenization-based multiscale µFE simulations. The constitutive laws were fitted by using a particle swarm optimization scheme (particleswarm, MATLAB R2017b, MathWorks Inc.), followed by a gradient-based scheme (fmincon, MATLAB R2017b, MathWorks Inc.) to enhance the final tuning of the parameters, as it is assumed that when particleswarm finishes, the solution is already within the proximity of a minimum. The minimization problem is thus defined as

$$\min \sum\_{i=1}^{n} \left( \| \left[ \mathbb{E}\_{\text{pred}}(\theta\_{i}) - \mathbb{E}\_{\mu \text{FE}} \right] \|\_{i} \right)^{2},\tag{35}$$

where n is the number of samples×load cases×strain levels, which means that the damage results for each sample, each considered load case, and each considered strain level (i.e., 0.2, 0.3, 0.4, and 0.5%) are used in the parameter fitting procedure; k[Epred]k is the Frobenius norm of the matrix projection of the damaged stiffness tensor predicted by the considered theoretical damage model, k[EµFE]k is the Frobenius norm of the matrix projection of the damaged stiffness tensor calculated through homogenization, and θ are the s different parameters of the considered damage model.

This minimization problem (Equation 35) involves the fitting of parameters which govern the size of the damage dissipation potentials (i.e., the surface containing the elastic regime, in which damage does not develop; it is the damage analog to the yield surface in plasticity), and therefore the solution of the CPPM scheme may involve negative 1γk , which are not physical solutions. The CPPM scheme is used in computational plasticity and/or damage contexts to solve the corresponding non-linear equations (Equation 20). If the loading state of a sample is found within the elastic regime (i.e., inside of the yield surface in a plasticity context, or inside the damage dissipation potential in a damage context), no equations need to be solved as plasticity and/or damage related quantities would not further develop. Thus, these undesired values of 1γk will arise only if the loading state of the considered sample is not outside of the damage dissipation potential. In order to avoid these, the minimization problem is modified with a penalty term to avoid such unwanted situations, such that

$$\min \sum\_{i=1}^{n} \left[ \left( \| \left[ \mathbb{E}\_{\text{pred}} (\theta\_s) - \mathbb{E}\_{\mu \text{FE}} \right] \| \right)\_i \right]^2 + \sum\_{k=1}^{l} \mathbb{H} (-\Delta \chi\_{i,k}) K\_{\text{perm}} (\mathbf{e}^{\left| \Delta \chi\_{i,k} \right|} - 1) \right],\tag{36}$$

where Kpen is a large (penalty) constant.

The initial choice of a solver not based on gradients is because the addition of this penalty term breaks the C 1 continuity of the functional to be minimized, and its global non-convexity is assumed a priori. The specific choice of particle swarm optimization over other methods not based on gradients, such as genetic algorithm, is established on the superior computational efficiency of particle swarm optimization over the genetic algorithm (Panda and Padhy, 2008).

The goodness of the fitting procedure was analyzed with the standard error of the estimate (SEE). This is calculated as

$$\text{SEE}(\%) = 100 \frac{\sqrt{\sum\_{i=1}^{n} (\|\left[\mathbb{E}\_{\text{pred}} - \mathbb{E}\_{\mu \text{FE}}\right]\|\_{i})^2}}{\sqrt{\sum\_{i=1}^{n} (\|\left[\mathbb{E}\_{\text{pred}} - \mathbb{E}\_{0}\right]\|\_{i})^2}}. \tag{37}$$

#### 4. RESULTS

#### 4.1. Evaluation of the µFE Results

For all load cases in **Table 1**, the considered samples were subjected to several strain levels, leading to different damage levels. The resulting macroscopic damaged stiffness tensors and the macroscopic strain Frobenius norms were measured at 0.2, 0.3, 0.4, and 0.5% strain levels by using the 0.2% strain criterion (Wolfram et al., 2012). This theoretically leads to damage and macroscopic strain Frobenius norms being evaluated, respectively, at 0–0.3% (with 0% being considered as macroscopic yield) macroscopic plastic strain Frobenius norms. The macroscopic strain Frobenius norms at 0.5% strain level for each load case are shown in **Figure 2** in the form of boxplots. It can be seen from this figure that within each group (T, C, S, or MA), higher macroscopic strain Frobenius norms correspond to compression-dominated load cases (load cases 4–6, 13, 17, and 21 in **Figure 2**).

Damage is evaluated by subtracting the damaged stiffness tensor from the undamaged stiffness tensor and calculating the Frobenius norm of its matrix projection (k[E0 − Edam]k). The values of these norms for each of the considered load cases are shown in **Figure 3**; the damage shown corresponds to the 0.5% strain level. Due to the alignment of the samples and ordering of their fabric eigenvalues (m1 > m2 > m3), it can be seen from this figure that within each group (T, C, S, or MA), higher damage values are seen where the fabric tensor eigenvalues are the largest (i.e., load cases 1, 4, 7, and 10–13 in **Figure 3**). Moreover, higher damage values are also seen in load cases that are compression-dominated (load cases 4–6, 13, 17 and 21). These higher damage values in uniaxial compression, or in compressive-dominated multi-axial load cases, compared to tension load cases indicate a possible

tension/compression asymmetry in the damage behavior at the macroscopic level. It is important to mention that, although damage values were measured at the same strain levels according to the 0.2% strain criterion, the macroscopic strain Frobenius norms (**Figure 2**) were considerably larger in compression than in tension.

Multi-linear regressions in log-log space were performed to establish possible relationships between damage and the microarchitectural indices of the considered samples. These regressions were between k[E0 − Edam]k at 0.5% strain level, BV/TV, fabric eigenvalues and macroscopic strain Frobenius norms, such as

$$\begin{aligned} \log(\|\left[\mathbb{E}\_0 - \mathbb{E}\_{\text{dam}}\right]\|) &= A + B \log(\text{BV/TV}) + C \log(m\_1) \\ &+ D \log(m\_2) + E \log(\|\mathfrak{s}\_0\|) \end{aligned} \tag{38}$$

where m1 and m2 are the fabric eigenvalues corresponding to directions 1 and 2 (only shear and multi-axial load cases have two directions); A, B, C, D, and E are the constants in the regression. These regressions were performed separately for the following

TABLE 3 | Results from the multi-linear regressions between k[E0 − Edam]k at 0.5% strain level, BV/TV, fabric eigenvalues, and macroscopic strain Frobenius norms, in log-log space.

*Regressions were performed for uniaxial tension (T), uniaxial compression (C), combined uniaxial tension and uniaxial compression (T* ∪ *C), shear (S), and multi-axial (MA) in normal strain space.*

sets of load cases: uniaxial tension, uniaxial compression, combined uniaxial tension and uniaxial compression, shear, and multi-axial load cases in normal strain space. The results from these regressions can be seen in **Table 3**. **Table 3** shows that both BV/TV and fabric eigenvalues have a significant effect (p ≤ 0.05), and that damage expressed as per Equation (38) is directly proportional to the micro-architectural indices, with the slopes for BV/TV being substantially larger than those for the fabric eigenvalues. The coefficients of determination (R 2 ) show that only the multi-linear model of the multi-axial load cases in normal strain space behaves poorly in comparison to the rest.

The component-wise fraction between the matrix projection of E0 − Edam at 0.5% strain level and the matrix projection of E0 (i.e., the i-th and j-th component of E0 − Edam is divided by the i-th and j-th component of E0) leads to the 6 × 6 matrix with components

$$[\mathbf{D}]\_{ij} = \frac{[\mathbb{E}\_0 - \mathbb{E}\_{\text{dam}}]\_{ij}}{[\mathbb{E}\_0]\_{ij}}.\tag{39}$$

This matrix depicts the component-wise ratio of the damaged and undamaged coefficients for each sample and load case. The component-wise mean of [**D**]ij over all the considered samples was calculated and then normalized from 0 to 1 for each of the considered load cases, forming another 6 × 6 matrix (e.g., the new matrix i-th and j-th component is the mean of the **D**ij components of all the samples); the components in E0 which are zero are ignored and not considered in the normalization, i.e., the non-orthotropic coefficients. The resulting 21 normalized matrices are shown in **Figure 4**. These plots suggest that macroscopic damage in trabecular bone is actually anisotropic and dependent on the considered load case. In uniaxial tensile and compressive load cases, it can be observed that the normal components of the stiffness tensor which are related to the considered load case are the most affected ones (e.g., in the load case ε11 > 0, components E1111, E1122, E1133, and the corresponding symmetric counterparts are more affected than the rest). In shear load cases, the corresponding shear component is the most affected one. Considering multiaxial load cases in normal strain space we find that in tensiontension and compression-compression load cases, the most affected components are in the off-diagonals of the matrix the components that are related to the plane which is being loaded (e.g., in the load case ε11 = ε22 > 0, components E1122 and E2211 are more affected than the rest); in tensioncompression/compression-tension load cases, the most affected components are in the matrix diagonal - the components that are related to the plane which is being loaded (e.g., in the load case ε11 = −ε22 > 0, components E1111 and E2222 are more affected than the rest).

### 4.2. Effect of BV/TV and MIL Fabric Tensor on the Damage Behavior

The effect of BV/TV and fabric on the macroscopic damage behavior of trabecular bone was assessed by (1) considering the single scalar isotropic damage model in section <IP_ADDRESS> with and without considering the effect of BV/TV and then comparing the respective values of SEE; and (2) considering the anisotropic damage model in section <IP_ADDRESS> with and without considering the effects of BV/TV and fabric eigenvalues and then comparing the respective values of SEE. In the anisotropic scenario, in the case in which fabric eigenvalues were not included, the order of fabric eigenvalues was randomized to maximise the effect of including fabric in the comparison (the ordering no longer corresponds to m1 > m2 > m3; the corresponding stiffness and strain tensors were reordered accordingly). The minimization scheme was run for five times to ensure that a suboptimal solution was not chosen. This comparison is shown in **Table 4**.

Note that the values of SEE of the anisotropic cases are not considerable lower than those of the isotropic cases. This is because even if the damage is higher in the components related to the considered load case, all the components of the stiffness tensor are damaged, and k[E0 − Edam]k takes into account the reduction of all the components of the stiffness tensor. The exponents that express BV/TV dependency are considerably larger than those expressing fabric eigenvalue dependency.

### 4.3. Macroscopic Damage Model for Trabecular Bone

A damage model which incorporates both isotropic/anisotropic damage progression and tension/compression asymmetry was implemented and its efficacy in evaluating the macroscopic damage behavior of trabecular bone was assessed. BV/TV and fabric eigenvalue dependencies were considered; BV/TV dependency was included in the isotropic part of the model while both BV/TV and fabric eigenvalue dependencies were included in the anisotropic part. Tension/compression asymmetry was included as shown in section <IP_ADDRESS>. The SEE and the value of the parameters of the model are shown in **Table 5**.

This considered model reduces the SEE in more than 15% with respect to the single scalar isotropic model (SEE = 37.03%). Despite the 13 parameters, a considerably larger number in comparison with the two parameters of the isotropic model, the values of some of these parameters suggest that not all of them need to be considered. For

uniaxial tension, (D–F) correspond to uniaxial compression, (G–I) correspond to shear, and (J–U) correspond to multi-axial in normal strain space.

instance, the value of B0,aniso is very small, which means that these parameters, together with the corresponding exponents expressing BV/TV and fabric eigenvalue dependencies (u and r) could be ignored, reducing the number of parameters to 10. It is important to point out the negative values of η and B0,iso.

### 5. DISCUSSION

The macroscopic damage behavior of trabecular bone has been researched in a few studies, but these are usually restricted to uniaxial load scenarios which only permit the assessment of stiffness reduction in the direction of loading (Keaveny et al., 1994b; Zioupos et al., 2008; Garcia et al., 2009). Consequently these studies are unable to provide a comprehensive constitutive model that can be included in whole-bone simulations. This study investigated the possible relationship between damage at the tissue level and the macroscopic multi-axial damage behavior, by employing a homogenization-based multiscale approach to samples with a relatively wide range of BV/TV and fabric tensor eigenvalues, subjected to multiple loading scenarios. The macroscopic damage behavior of trabecular TABLE 4 | SEEs, BV/TV, and fabric eigenvalue exponents for the isotropic and anisotropic models.

*Models 1 and 2 are isotropic with and without BV/TV dependency, respectively; models 1, 2, and 3 are anisotropic and: (1) without BV/TV and fabric eigenvalue dependencies, (2) with BV/TV dependency only, and (3) with fabric eigenvalue dependency only.*

TABLE 5 | Value of the parameters and SEE of the combined isotropic/anisotropic model with tension/compression asymmetry.

bone was approximated via different continuum damage models: isotropic and anisotropic; with and without BV/TV and fabric eigenvalue dependencies; and with and without tension-compression asymmetry. From the results, it can be concluded that the macroscopic damage behavior of trabecular bone has the following features: BV/TV and fabric eigenvalue dependencies; tension/compression asymmetry; a combined isotropic/anisotropic behaviour. The first two of these features are not unexpected as they play a key role in the evaluation of elastic stiffness (Odgaard et al., 1997; Zysset, 2003), however, the previously unexplored, last feature indicates that damage in trabecular bone is best represented by using both isotropic and anisotropic damage variables. This is likely to be true for most cellular materials.

This study assumed an isotropic model with coupled damage and plastic behavior at the tissue level, which was deemed appropriate as the isotropy assumption at this level is known to result in little to no error in macroscopic results (Cowin, 1997). Isotropic damage at the solid phase level leads to an anisotropic macroscopic damage response with a dependency on the considered load case (Levrero-Florencio et al., 2017a). The variation in the components of the stiffness tensor shows anisotropic damage which depends on the considered load case (**Figure 4**). Shi et al. (2010) suggested that there is a larger proportion of damaged tissue in the longitudinal trabeculae (direction of loading) for uniaxial load cases, which is in agreement with the results presented here, as the most damaged components of the macroscopic stiffness tensor are always the on-axis components. An issue which may make validation of these results very challenging is the use of kinematic uniform boundary conditions; these boundary conditions are extremely difficult, not to say impossible, to reproduce experimentally, especially for the more complex load cases. Most previous studies involving damage in trabecular bone have used isotropic models (Garcia et al., 2009; Schwiedrzik and Zysset, 2013), which may be acceptable for proportional loading scenarios, but not for changing loads or cyclic loading scenarios, such as those arising during physiological activities.

The results show that the macroscopic strain Frobenius norms were considerably larger in macroscopic compression than in macroscopic tension. This is important in the considered context of damage modeling as the thermodynamic stress-like variables governing damage evolution (Yk ) directly depend on the macroscopic strain values, which could explain the higher damage values in compression without the explicit need of modeling tension/compression asymmetry. However, this asymmetry is taken into account because it still leads to a better fit of the damage model and it only consists of one additional parameter. The fact that damage values are higher in compression-dominated load cases compared to tension load cases could be related to the more heterogeneous stress distributions at the solid phase level occurring during macroscopic compression, which includes tensile stresses at the tissue level due to bending and buckling of trabeculae (Stölken and Kinney, 2003). Another important factor to take into account is that the considered model at the tissue level is ductile (i.e., fracture is not incorporated). If fracture was considered at a critical damage threshold, the tension/compression asymmetry would probably be different as tissue damage is more diffused in compression than in tension (Lambers et al., 2014), and therefore a significant decrease of load carrying capacity would occur in tension.

The variation in the components of the stiffness tensor shows anisotropic damage which depends on the considered load case (**Figure 4**). Shi et al. (2010) suggested that there is a larger proportion of damaged tissue in the longitudinal trabeculae (direction of loading) for uniaxial load cases, which is in agreement with the results presented here, as the most damaged components of the macroscopic stiffness tensor are always the on-axis components. An issue which may make validation of these results very challenging is the use of kinematic uniform boundary conditions; these boundary conditions are extremely difficult, not to say impossible, to reproduce experimentally, especially for the more complex load cases. Most previous studies involving damage in trabecular bone have used isotropic models (Garcia et al., 2009; Schwiedrzik and Zysset, 2013), which may be acceptable for proportional loading scenarios, but not for changing loads or cyclic loading scenarios, such as those arising during physiological activities.

Multi-linear regressions between k[E0 − Edam]k, BV/TV, fabric eigenvalues and macroscopic strain Frobenius norms (from **Table 3**). It shows that both BV/TV and fabric eigenvalues are statistically significant. The coefficients of determination suggest that only the regression of k[E0 − Edam]k of the multi-axial load cases in normal strain space behaved poorly in comparison to the others. The slopes of BV/TV are significantly higher than those of fabric eigenvalues, suggesting that BV/TV plays a more important role in these regressions; they also suggest that the higher the BV/TV and fabric eigenvalues, the higher the damage is. Results in Levrero-Florencio et al. (2017a) showed that the damage in the orthotropic coefficients of the macroscopic stiffness tensors do not have significant dependencies on BV/TV or fabric, for each of the considered load cases. In this study the Frobenius norm k[E0 − Edam]k is used instead, which takes into account the damage of all the components of the macroscopic stiffness tensor. Therefore, the slopes and p-values in **Table 3** suggest that lower BV/TV samples have a more anisotropic damage behaviour in the sense that the longitudinal trabeculae are more damaged than the oblique, and that higher BV/TV samples have a more isotropic behavior, or are more damaged in general. Even if fabric eigenvalues have a significant effect on k[E0 −Edam]k, the considerably lower slopes suggest that their relevance is significantly lower than that of BV/TV.

The standard errors of the estimate (SEE) and the exponents with respect to BV/TV and fabric eigenvalues of five different damage models indicate that the SEEs are not substantially different in all these considered models, this is because, despite the anisotropic damage behavior, all the components of the stiffness tensor are damaged (Levrero-Florencio et al., 2017a), suggesting that while a combined isotropic and anisotropic model is most suitable for simulating the macroscopic damage behavior of trabecular bone, an isotropic model is not necessarily poor. The SEEs of the models with dependencies are not substantially lower to those without the dependencies, suggesting that the considered BV/TV and fabric eigenvalue dependencies may not be needed. Nonetheless, the results of the multi-linear regressions (**Table 3**) show significance of BV/TV and fabric eigenvalues when modeling damage. Furthermore, since these five assessed damage formulations only partially model some of the features of the macroscopic damage behavior of trabecular bone mentioned earlier, the dependencies are maintained in the combined isotropic/anisotropic model with tension/compression asymmetry.

It is apparent that the model with a combined isotropic/anisotropic behavior and tension/compression asymmetry is a substantial improvement over the single scalar damage formulation since the SEE is reduced by more than 15% (**Table 5**). Nonetheless, it is important to mention that this model has 13 parameters instead of 2, though the value of the parameter B0,aniso indicates that this parameter and the associated exponents expressing BV/TV and fabric eigenvalue dependencies can be ignored. The negative value of η suggests that if tension-dominated cases had similar strains to those in compression-dominated cases, the damage values would be higher in tension, as a negative value of η implies crack-closure, which is expected as bone could be considered a quasi-brittle material (Hambli, 2013; Mayya et al., 2016). The negative value of B0,iso suggest that, when modeling the damage progression with a linear model, there is an initial presence of damage, which has been previously observed in Levrero-Florencio et al. (2017a) (the intercepts of the y-axis of the damage-accumulated plastic strain plots are not zero).

This study has a number of limitations. As previously mentioned, bone at the solid phase level is assumed to be ductile, i.e., while reduction in stiffness due to damage is included, fracture is not. This is perhaps appropriate for the considered level of loading, but it is indeed not applicable if large strains are applied, as complete fracture of trabeculae can occur. Nawathe et al. (2013) shows that ductile tissue behavior overestimates the experimental yield properties. Another limitation, previously stated in Levrero-Florencio et al. (2017a), is that although there is plenty of experimental data on uniaxial load cases (Keaveny et al., 1997; Bayraktar and Keaveny, 2004; Sanyal et al., 2012; Manda et al., 2016), these physical experiments do not allow evaluation of stiffness for samples subjected to different load cases and the effect of loading in one direction on the behavior in the others. Therefore, a study completely based on numerical simulations is the only alternative even though the results cannot be currently validated experimentally. The use of kinematic uniform boundary conditions in the µFE analyses could also be considered a limitation, as they are known for providing an upper bound of the stiffness tensor (Pahr and Zysset, 2008; Wang et al., 2009) or macroscopic yield (Panyasantisuk et al., 2015), and may also affect the damage morphology when compared to the in situ case (Daszkiewicz et al., 2017). We also assume that the orthotropic directions do not rotate during loading, which may be a valid assumption for the considered range of strains.

Use of a large number of load cases (21) and samples (10) shows that the evolution of the damaged macroscopic stiffness tensor is based on the loading history. By examining relationships between bone microstructural indices (such as BV/TV and fabric) with macroscopic damage constitutive laws, we show that the proposed combined isotropic/anisotropic damage law with tension/compression asymmetry is a viable superior alternative to the widely used single scalar isotropic damage formulation as it reduces the fitting error from 37 to 22%; it does, however, require specification of a larger number of material parameters. The relationships of damage progression with bone's microarchitectural indices (density and fabric) developed in this study provide an approach for the creation of macroscale continuum models that incorporate damage and will, therefore, improve clinical predictions of the behavior of bone and bone-implant systems.

### AUTHOR CONTRIBUTIONS

FL-F designed the study, performed the FE simulations and parameter fittings, and analyzed the data; PP contributed to the design of the study. Both authors contributed to the critical writing and revision of the manuscript.

#### FUNDING

This work was supported by funding from the Engineering and Physical Sciences Research Council grant EP/K036939/1. The authors gratefully acknowledge ARCHER, UK National Supercomputing Service, for access to their Cray XC30

#### REFERENCES

supercomputer under the project "Modelling the nonlinear micromechanical behaviour of bone".

#### ACKNOWLEDGMENTS

We gratefully acknowledge Dr. Lee Margetts, from The University of Manchester, for his assistance with the implementation and development of the used version of ParaFEM.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Levrero-Florencio and Pankaj. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Left Ventricular Trabeculations Decrease the Wall Shear Stress and Increase the Intra-Ventricular Pressure Drop in CFD Simulations

Federica Sacco1,2 \*, Bruno Paun2 , Oriol Lehmkuhl 1 , Tinen L. Iles 3 , Paul A. Iaizzo3 , Guillaume Houzeaux 1 , Mariano Vázquez 1,4, Constantine Butakoff 2 and Jazmin Aguado-Sierra1 \*

1 Barcelona Supercomputing Center (BSC), Barcelona, Spain, 2 PhySense, ETIC, Universitat Pompeu Fabra, Barcelona, Spain, 3 Visible Heart Laboratory, Department of Surgery, University of Minnesota, Minneapolis, MN, United States, 4 IIIA - CSIC, Bellaterra, Spain

#### Edited by:

Timothy W. Secomb, University of Arizona, United States

#### Reviewed by:

Patrick Segers, Ghent University, Belgium Gernot Plank, Medizinische Universität Graz, Austria

#### \*Correspondence:

Federica Sacco<EMAIL_ADDRESS>Jazmin Aguado-Sierra<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology

Received: 27 November 2017 Accepted: 13 April 2018 Published: 30 April 2018

#### Citation:

Sacco F, Paun B, Lehmkuhl O, Iles TL, Iaizzo PA, Houzeaux G, Vázquez M, Butakoff C and Aguado-Sierra J (2018) Left Ventricular Trabeculations Decrease the Wall Shear Stress and Increase the Intra-Ventricular Pressure Drop in CFD Simulations. Front. Physiol. 9:458. doi: 10.3389/fphys.2018.00458

The aim of the present study is to characterize the hemodynamics of left ventricular (LV) geometries to examine the impact of trabeculae and papillary muscles (PMs) on blood flow using high performance computing (HPC). Five pairs of detailed and smoothed LV endocardium models were reconstructed from high-resolution magnetic resonance images (MRI) of ex-vivo human hearts. The detailed model of one LV pair is characterized only by the PMs and few big trabeculae, to represent state of art level of endocardial detail. The other four detailed models obtained include instead endocardial structures measuring ≥ 1 mm2 in cross-sectional area. The geometrical characterizations were done using computational fluid dynamics (CFD) simulations with rigid walls and both constant and transient flow inputs on the detailed and smoothed models for comparison. These simulations do not represent a clinical or physiological scenario, but a characterization of the interaction of endocardial structures with blood flow. Steady flow simulations were employed to quantify the pressure drop between the inlet and the outlet of the LVs and the wall shear stress (WSS). Coherent structures were analyzed using the Q-criterion for both constant and transient flow inputs. Our results show that trabeculae and PMs increase the intra-ventricular pressure drop, reduce the WSS and disrupt the dominant single vortex, usually present in the smoothed-endocardium models, generating secondary small vortices. Given that obtaining high resolution anatomical detail is challenging in-vivo, we propose that the effect of trabeculations can be incorporated into smoothed ventricular geometries by adding a porous layer along the LV endocardial wall. Results show that a porous layer of a thickness of 1.2 · 10−2 m with a porosity of 20 kg/m2 on the smoothed-endocardium ventricle models approximates the pressure drops, vorticities and WSS observed in the detailed models.

Keywords: trabeculae, papillary muscles, left ventricular modeling, left ventricular hemodynamics, porosity

### 1. INTRODUCTION

Computational cardiac modeling has become important as a non-invasive modality to study the overall cardiac function (Trayanova, 2011; Taylor et al., 2013). Recently, regulatory bodies are encouraging and supporting the use of in-silico modeling to reduce animal experimentation. Within this context, models of cardiac hemodynamics have yet to be improved. The majority of the hemodynamic cardiac computational simulations consider simplified geometries with smoothed endocardial surfaces (Doost et al., 2015; Khalafvand et al., 2015; Imanparast et al., 2016), mostly due to a lack of high-resolution, fast and safe in-vivo imaging techniques. It is also true that solving highly detailed models require computationally expensive simulations that can only be carried out using HPC.

In reality the heart anatomy is complex and all individuals have their own unique anatomies. The interior of the cardiac chambers is not smooth: it is populated by PMs, trabeculae of different sizes and false tendons (Gao et al., 2014). In the LV, PMs are the muscles responsible for properly positioning the chordae tendinae during systole to optimize mitral valve leaflet coaptation. Trabeculae are complex muscular structures that are unique to a given human heart, mostly consisting of myocytes, that protrude from the endocardial wall into the interior of the ventricle and present a sponge-like structure. The primary role of the trabeculae in the overall cardiac function remains unknown, but they are often associated with the Purkinje network.

State of the art LV CFD simulations, employing detailed endocardial structure models, have been created from either MRI or computed tomography (CT) in-vivo modalities (Chnafa et al., 2016; Lantz et al., 2016; Vedula et al., 2016), but they only incorporated PMs and a few large trabeculae. In previously reported studies, Chnafa et al. (2016) used 4D MR images to reconstruct the LV geometry, characterized only by PMs, and prescribe physiological deformations using numerical treatments. In this way the author could study blood flow instabilities within the ventricular cavity and found out that high-frequency flow fluctuations can be common in normal LVs. Both Vedula et al. (2016) and Lantz et al. (2016) added also few big trabeculae together with the PMs in their LV geometries and studied the impact of these endocardial structures on the blood flow by comparing simulations results with smoothedendocardium ventricles. Vedula et al. (2016) reconstructed the LV geometry from high-resolution 4D CT scans and applied prescribed mesh deformation based on immersed boundary method. The authors observed a "scrambling" of blood flow vortices produced by PMs and trabeculae, which caused the generation of deeper and more complex vortices that were not present in the smoothed model. In this way, trabeculations help diastolic filling and, during systole, they help ventricular washout by wringing out the blood flow out from the apex. Lantz et al. (2016) extracted the LV endocardial surface and wall motion over time from 4D CT data. In contrast with Vedula et al. (2016) and Lantz et al. (2016) did not observe any deep penetration of the mitral inflow jet toward the apex: the jet strongly interacted with the PMs and was diverted toward the outflow tract. However, both papers have demonstrated that the detailed anatomies of LV endocardium have important influences on blood flow dynamics: in particular, particle tracking used by Lantz et al. (2016) demonstrated that blood flow interacted with trabeculae and PMs, creating vortices around the endocardial spaces between the trabeculations. More vortices appeared during diastole in the detailed LV as compared to the smoothed one, and PMs redirected blood flow and generated a large vortex, which was not present in the smoothed model. Finally, it was shown that the presence of trabeculations created a region where the flow appeared to be stagnant during five cardiac cycles, which is impossible to reproduce with smoothed endocardium models. While Vedula et al. (2016) and Lantz et al. (2016) considered the PMs and few trabeculae, the level of detail and the amount of trabecular structures in LV geometry reconstructions was not as high as in Kulp et al. (2011), who segmented detailed endocardial structures from high-resolution 4D CT data. The authors studied the interaction between trabeculations and the blood flow by deforming the initial 3D mesh in each following frame. Results showed how the complex endocardial surface caused the blood to move through the empty spaces between the trabeculations and fill these cavities during diastole.

In this paper we used highly detailed anatomical LV endocardium models to characterize the effects of trabeculae and PMs on the blood flows using CFD simulations. Four detailed LV geometries were reconstructed from high resolution imaging data of perfusion fixed human hearts (2 male and 2 female), which were obtained at the Visible Heart R Laboratory (Atlas of Human Cardiac Anatomy, RRID:SCR\_015734). Detailed and smoothed endocardial models were reconstructed for each of the four hearts to quantify the differences between these two cases and thus characterize the impacts of PMs and trabeculae on ventricular hemodynamics. The level of detail in these reconstructions of the endocardial structures was, to the best of our knowledge, the highest ever achieved for this kind of study: the average size of the smallest structures reconstructed measures about 1 mm2 in cross-sectional area.

A fifth (male) LV geometry, named control LV, was reconstructed from the human hearts high resolution images dataset, together with its smoothed equivalent. This model was only characterized by PMs and few large trabeculae: in this way we could compare simulation results obtained from the highly detailed models described previously to the ones from an LV geometry which is similar in detail to those present in literature (Kulp et al., 2011; Chnafa et al., 2016; Lantz et al., 2016; Vedula et al., 2016).

Through CFD simulations we aim to characterize the hemodynamics inside detailed vs. smoothed human ventricular anatomies by quantifying the trabecular volume, intraventricular pressure drop, WSS and vorticity within the LV cavities. Furthermore, we propose that a porous layer can be added to the LV endocardium to compensate for the absence of

**Abbreviations:** LV, Left Ventricle; PMs, Papillary Muscles; HPC, High Performance Computing; MRI, Magnetic Resonance Imaging; CFD, Computational Fluid Dynamics; WSS, Wall Shear Stress; CT, Computed Tomography; BSC, Barcelona Supercomputing Center; PBS, Phosphate Buffered Saline; FE, Finite Element; FSI, Fluid Structure Interaction.

trabeculae within smoothed ventricular hemodynamic models. Our main findings show that the presence of trabeculae alters significantly the blood flow by increasing intra-ventricular pressure drop, reducing the shear stress at the ventricular walls and generating multiple secondary vortices, absent in smoothwalled ventricle simulations. Furthermore, our results show that indeed a porous layer can compensate for the absence of trabeculations inside simplified ventricular models by increasing the intra-ventricular pressure drop, by reducing the wall shear stress at the interface of the porous layer and by increasing the amount of vortex structures within the LV.

#### 2. MATERIALS AND METHODS

#### 2.1. Left Ventricular Models

The five LV models used in this work were reconstructed from high-resolution MR images obtained from in vitro perfusion of fixed human hearts. The research uses of these heart specimens have received appropriate approval from both the University of Minnesota's Institution Review Board and LifeSource Research Committee (Minnesota's non-profit procurement donation organization). The hearts were recovered from organ donors whose hearts were not viable for transplantation. Written and informed consents were obtained from the donors families which follow the wishes of the donor. The database is open to public access.

DICOM data sets were acquired utilizing a 3T Siemens scanner with 0.44 × 0.44 mm in-plane resolution and slice thickness of 1 to 1.7 mm. The hearts were fixed with 10% formalin in phosphate buffered saline (PBS) solution for at least 24 h under 40–50 mmHg of pressure and then stored in 10% formalin. The five hearts DICOM datasets are shown in **Figure 1**.

Image segmentation was carried out with Fiji software (Fiji, RRID:SCR\_002285), using the maximum entropy-based thresholding algorithm (Qi, 2014), followed by endocardial surface reconstruction using marching cubes algorithm in Seg3D (Seg3D, RRID:SCR\_002552). The relative high contrast of the images guaranteed that the thresholding produced detailed endocardial models. The smoothed models were generated from the detailed geometries by manually deleting trabeculae and PMs and closing holes on associated surfaces using ReMESH software (ReMESH, RRID:SCR\_015735). Autodesk Meshmixer (Autodesk Meshmixer, RRID:SCR\_015736) sculpting software was then used to adjust the smoothed endocardial surface as to maintain the same outline for both the smoothed and detailed geometries.

The control LV, as a representative of a state of the art anatomical model, was reconstructed using a regularized region growing algorithm of ITK-SNAP (ITK-SNAP Medical Image Segmentation Tool, RRID:SCR\_002010) to get only large scale anatomical detail from the images. The algorithm allowed controlling the smoothness of the extracted contour making it easier to obtain the smoothed surface with just the PMs and a few large trabeculae (approximately 5 mm2 in cross-sectional area). The obtained level of detail for the control LV was similar to the reported models used in recent publications on blood flow analysis in LV such as Vedula et al. (2016) and Chnafa et al. (2016).

In order to let the flow develop, a 50 mm long tube was attached at the inlet (corresponding to the mitral valve orifice) and a 70 mm tube at the outlet (corresponding to the aortic valve orifice). Each tube base matched exactly the corresponding valvular ring plane. Tubes were created for every given LV model using ParaView (ParaView, RRID:SCR\_002516).

The resulting surface meshes were uniformly remeshed using Remesh and then volumetric tetrahedral meshes were generated using an isosurface-stuffing-based algorithm (Labelle and Shewchuk, 2007) with an in-house mesher developed at the Barcelona Supercomputing Center (BSC). The volumetric meshes had adaptive element size, with volumes varying from 10−7 mm3 to 1.9 · 10−2 mm3 , with an average size of 5.7 · 10−5 mm3 . Wireframe zoomed images of the tetrahedral meshes can be found in Figure S1.

The five LV anatomies, both smoothed and detailed, are shown in **Figure 2**. A more detailed view is reported in Figure S2. The medical histories and related information can be found in Table 1 of the Supplementary Materials. The ventricular volume of each mesh is reported in Table 2 of the Supplementary Materials.

#### 2.2. Hemodynamic Simulations

To carry out CFD simulations the walls were defined as rigid, noslip boundary conditions. For the outlet, a stabilizing boundary condition employed a baseline pressure of 10.7 kPa (80 mmHg, a normal end-diastolic arterial pressure) plus an outflow resistance (Bazilevs et al., 2009). Blood viscosity was set to 0.0035 kg/(m · s) and density to 1,060 kg/m3 .

Hemodynamic simulations solving continuity and Navier-Stokes equations for incompressible flows were run on the MareNostrum 4 supercomputer (MareNostrum, RRID:SCR\_015737) and on Archer (ARCHER, RRID:SCR\_015854), UK supercomputer, using Alya, the BSC's in-house, parallel multi-physics, HPC solver (Houzeaux et al., 2009; Vazquez et al., 2016).

Simulations were carried out using a low dissipation finite element (FE) strategy described below. The Navier-Stokes equations for a fluid domain bounded by Ŵ = ∂ within the time interval (t0, tf ) reside in calculating a velocity **u** and a kinematic pressure p so that Equations (1, 2) are satisfied; where ν is the kinematic viscosity, **f** is the vector of external body forces and **S**(**u**) is the rate-of-strain tensor.

$$
\partial\_t \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} - 2\nu \nabla \cdot \mathbf{S}(\mathbf{u}) + \nabla p - \mathbf{f} = \mathbf{0} \qquad \text{in } \Omega \times (t\_0, t\_f)
$$

$$\{\}\tag{1}$$

$$\nabla \cdot \mathbf{u} = 0 \qquad \text{in } \Omega \times (t\_0, t\_f) \tag{2}$$

To obtain weak or variational formulation of the Navier-Stokes equations Equations (1, 2), we introduced the spaces of vector functions **V**D= **H**1 D (), **V**0= **H**1 0 () and Q=L 2 () /ℜ. L 2 () is the space of square-integrable functions, H1 () is a subspace of L 2 () formed by functions whose derivatives belong also to

L 2 (), H1 D () is a subspace of H1 () which satisfies Dirichlet boundary condition on Ŵ. H1 0 () is a subspace of H1 () whose functions are zero on Ŵ; and **H**1 D () and **H**1 0 () are their vector counterparts in a two- or three-dimensional space. (·, ·) determines the standard L 2 inner product. For the evolutionary case, **V**t≡L 2 t0, tf ; **V**D and Qt≡D′ t0, tf ; Q were introduced, where L p t0, tf ; X is the space of time dependent functions in a normed space X so that R tf t0 f p X dt < ∞, 1 ≤ p < ∞ and Qt consists of mappings whose Q-norm is a distribution in time. The weak form of problem (Equations 1, 2) with the boundary conditions is then: Find **u** ∈ **V**t , p ∈ Qt such that Equation (3) is satisfied for every **v**, q ∈ **V**0 × Q.

$$\begin{aligned} \left(\partial\_l \mathbf{u}, \mathbf{v}\right) &+ \left(\mathbf{u} \cdot \nabla \mathbf{u}, \mathbf{v}\right) + 2\nu \langle \mathbf{S} \mathbf{u}, \nabla \mathbf{v} \rangle - \left(p, \nabla \cdot \mathbf{v}\right) + \left(q, \nabla \cdot \mathbf{u}\right) \\ &- \left(\mathbf{f}, \mathbf{v}\right) = \mathbf{0}, \end{aligned} \tag{3}$$

Moreover, in the previous equations the non-linear term convective form reported in Equation (4) was used, which is the most frequent choice in computational practice. Using Equation (2), other non-linear term forms can be derived, which are the same at the continuous level but do have different properties at the discrete level. In Equation (5) we consider the energy, momentum and angular momentum conserving form recently proposed in Charnyi et al. (2017). A non-incremental fractional-step method was used for pressure stabilization. This allows the use of finite element pairs which do not satisfy the inf-sup conditions, like the equal order interpolation for the velocity and pressure used in this work. An energy conserving Runge-Kutta explicit method lately proposed by Capuano et al. (2017) along with an eigenvalue based time-step estimator (Trias and Lehmkuhl, 2011) were used in order to time integrate the set of equations. This methodology, recently proposed by

Lehmkuhl et al. (2017), follows the principles of Verstappen and Veldman (2003), generalized for unstructured finite volumes by Jofre et al. (2013) and Trias et al. (2014) but in a FEM framework. The presented methodology has been successfully validated and benchmarked vs. other popular CFD approaches and experimental data in the bioengineering flows environment in Koullapis et al. (2017).

$$NL\_{conv}\left(\mathbf{u}\right) = \mathbf{u} \cdot \nabla \mathbf{u} \tag{4}$$

$$NL\_{emac}\left(\mathbf{u}\right) = 2\mathbf{S}\left(\mathbf{u}\right)\mathbf{u} + \left(\nabla \cdot \mathbf{u}\right)\mathbf{u} \tag{5}$$

We performed the following simulations:

$$\nu\left(t\right) = \begin{cases} \frac{A\_E}{2} \left(1 + \cos\left(\frac{2\pi\left(t - t\_{\bar{p},E}\right)}{t\_{1,E} - t\_{0,E}}\right)\right) & t\_{0,E} \le t \le t\_{1,E} \\ 0 & t\_{1,E} < t < t\_{0,A} \\ \frac{A\_A}{2} \left(1 + \cos\left(\frac{2\pi\left(t - t\_{\bar{p},A}\right)}{t\_{1,A} - t\_{0,A}}\right)\right) & t\_{0,A} \le t \le t\_{1,A} \end{cases} \tag{6}$$

$$
\mu = -\frac{1}{\mu} k(\nabla p) \tag{7}
$$

#### 2.3. HPC Characteristics

HPC characteristics for both constant and transient inflow simulations in terms of cores, total simulation time and time step are reported in **Tables 1**, **2**. Every simulation, with both constant and transient inflow, was run up to 800 ms. Information on the scalability of the incompressible flow module within Alya multiphysics solver can be found in the works of Houzeaux et al. (2009) and Vazquez et al. (2016). The elements-per-core ratio that was used to run these hemodynamic simulations was about 25,000.

#### 2.4. Geometric Markers

The following geometrical markers were used:

Trabecular volume was calculated as a difference between the volume of the convex hull of the detailed-endocardium LV model and the volume of the model itself.

The angle between inlet and outlet was the angle between the vectors normal to the two valvular planes (mitral and aortic). An illustration can be seen in **Figure 5**.

The distance between inlet and outlet was the distance between the mitral and aortic valves centers. An illustration can be seen in **Figure 5**.

#### 2.5. Hemodynamic Analysis

#### 2.5.1. Intra-ventricular Pressure Drop

The pressure distributions were analyzed within 15 mm long volumes of both inlet and outlet tubes. The sections were chosen right at the inlet of the mitral and at the outlet of the aortic valves. The histograms were normalized to unit area under the curve and bin width was calculated using the Freedman-Diaconis rule. As the pressure distributions were non-Gaussian (**Figure 6**), the intra-ventricular pressure drop was calculated as the difference between the inlet and outlet pressure mode. The pressure difference was then averaged over the ten last time frames (approximately 50 ms of simulations), during steady flow. From these results, we calculated the intra-ventricular pressure drop difference (1Pdiff ) as the difference of the detailed and smoothed pressure drops for every studied LV.

#### 2.5.2. WSS on the Ventricular Walls

WSS histograms were computed for every LV cavity (**Figure 7**) and normalized using the Freedman-Diaconis rule to choose the width of the bins. Given that these distributions were upward skewed, the median was used to analyze the WSS of each model. The total magnitude range and the mode are also reported for each case. In the case of the porous layer simulations, the median and mode WSS were calculated on the interface between the porous layer and the blood flow.

#### 2.5.3. Vorticity

Coherent structures were analyzed in both steady and transient inflow simulations applying the Q-criterion method (Hunt et al., 1988; Chakraborty et al., 2005). The applied thresholds for vortex visualization were 5,000 s−2 and 1,000 s−2 for steady and

TABLE 1 |Constant inflow simulations results and HPC information on both constant and transient inflow

simulations.

 1Pdiff=1Pd−1Ps, differences between detailed and smoothed pressure

 aIn the case of thecontrol LV we calculated the PMs volume, since this LV geometry was only characterized by PMs.

drops.

 bC.I., Constant Inflow simulations.cT.I.,TransientInflowsimulations.

transient inflow simulations respectively. Vortex quantification was done in Paraview, by integrating the contours of the vortices to estimate the total surface area.

#### 3. RESULTS

#### 3.1. Constant Inflow

The results of constant inflow quantification for the five LVs are shown in **Table 1**. In all LVs, the intra-ventricular pressure drops (1Pdiff ) increased in the detailed geometries by an average 0.2 kPa, except for subject D, that elicited the highest 1Pdiff of 2.7 kPa. Also the detailed E LV, exhibited a similar pressure drop to models A, B, and C, regardless of the low percentage of trabecular volume within its geometry. The geometrical markers do not indicate a direct correlation to the pressure drop (see **Table 1**). The pressure drop in models A, B, C, and E correlated best with the distance between their inlet and outlet, however, when model D is included, any good correlation disappears. The pressure drop is not correlated either to the Reynolds number at the inlet of each model.

The magnitudes of the WSS for each case is shown in **Figure 7**. The WSS histograms are shown in **Figure 8**; they were cropped at 3 Pa for visualization purposes, but the maximum values are included in each plot. The median WSS decreased in the detailed geometries on all subjects, except for subject D, as shown in **Table 1**. Notice that in model D, in **Figure 7**, high WSS regions are markedly localized on the PMs and the outflow region. In model D the WSS median values remain relatively similar between the smoothed and detailed models, even though the mode values are significantly lower for the detailed geometry. It is important to point out that the peak WSS was higher on the detailed models in comparison with the smoothed geometries, except for model E.

The vortical structures shown in **Figure 9** are thresholded at 5,000 s−2 and color coded according to the velocity magnitude [m/s]. The smoothed geometries generated fewer and larger vortices, while the detailed LVs showed the disruption of larger structures breaking down into a multitude of small scale vortices. To quantify them, the vortex contour surface area was calculated (shown in **Table 3**). The total surface area of the vortices was smaller in the smoothed (with a mean of 31.6 ± 8.9 · 10−3 m2 ) compared to the detailed LVs (with a mean of 43.5 ± 10.8 · 10−3 m2 ).

#### 3.2. E-A Wave Mitral Valve Inflow Results

From the E-A wave transmitral inflow function 6 time instants were selected, as highlighted in **Figure 4**: early E wave (1), E wave peak (2), late E wave (3), early A wave (4), A wave peak (5) and late A wave (6). **Figure 10** shows the vortices in model D. The Q-criterion values were thresholded at 1,000 s−2 and colored according to the velocity magnitude [m/s]. **Figure 10** and **Table 3** show that the presence of trabeculae created secondary vortices at the early E wave (time instants 1), increasing the total vortex surface from 9.4 · 10−3 m2 in the smoothed to 14.8 · 10−3 m2 in the detailed geometry. The secondary vortices in the detailed LV penetrated deeper between the trabeculations during the late E wave (time instant 3). During the early and peak A wave, a second weaker vortex ring was formed and it mixed with the vortices TABLE 2 | Constant inflow simulations results for A-D LVs with the porous layer of thickness 1.2 · 10−2 m and porosity 20 kg/m2 , along with the corresponding HPC information.

P, mode of the pressure distribution at inlet and outlet.

1P, intra-ventricular pressure drop from inlet to outlet.

1Pdiff = 1Pd − 1Ps, differences between detailed and smoothed pressure drops.

generated during the early filling (time instants 4–5). Here, for both smoothed and detailed geometries, the total surface of the vortices increased due to the mixing vortices but the amount was still higher in the detailed compared to the smoothed LV (25.5 − 22.5 · 10−3 m2 vs. 29.2 − 25.9 · 10−3 m2 respectively).

#### 3.3. Constant Inflow With a Porous Layer

The proposed porous layer produces energy dissipation, which increases the intra-ventricular pressure drop and adds complexity to the blood flow. The results of the sensitivity analysis of the effect on the intra-ventricular pressure drops for each thickness-porosity combination are reported in **Table 4**. By adding a layer of 1.2 · 10−2 m thickness and a porosity of

20 kg/m2 to the smoothed subject A, we obtained the intraventricular pressure drop equal to the one of its detailed case (1.5 kPa). Additionally, the vortex visualization using the Qcriterion thresholded at 5,000 s−2 demonstrates that the amount of vortices in the smoothed LV with the porous layer are similar to the ones in the detailed case (see **Figure 3**), presenting a total vortex surface area of 36.8 · 10−3 m2 . This may indicate that the roughness interacts with the boundary layer as a sandgrain roughness does, without prioritizing any particular flow direction, being then the proposed porous layer model very effective. The presence of the porous layer could also approximate the WSS calculated for the detailed geometry A as can be observed in **Table 2**, providing a relative error of only 0.076.

The same layer with a thickness of 1.2 · 10−2 m and a porosity of 20 kg/m2 was then applied to all the other cases (subjects B, C, D). In subject B and C, the presence of the porous layer increased the intra-ventricular pressure drop to values similar to the ones obtained within the detailed cases (see **Table 2**) providing a relative error of just 0.2 and 0.02 respectively. In subject D, the intra-ventricular pressure drop increased slightly with the presence of the porous layer, but in this case the thickness and porosity values of the porous layer were not able to reproduce the high values of pressure drop obtained inside the detailed model (relative error is 0.46). However, in all the models with the porous layer the WSS was reduced as shown in **Table 2** with the biggest relative error being 0.23 in model D. This table also shows how in all the cases the total vortex surface increased with the presence of the porous layer by providing values slightly higher than those of the detailed geometries in all cases, with the maximum relative error being −0.12 for model D.

#### 4. DISCUSSION

LV endocardiums of humans present a highly trabeculated appearance, which is often ignored in ventricular hemodynamic studies. Even though a few studies have been done to analyze the effect of papillaries and trabeculae on blood flow, no study to our knowledge has ever included small trabeculae of cross-sectional area of 1 mm2 .

In this study we focused on characterizing solely the effects of the geometries on the hemodynamics using CFD simulations. The fact that the walls are rigid makes it impossible to extapolate any of the findings to a clinical or physiologically relevant scenario. This study provides an engineering-like approach to a very complex biological system, by quantifying and characterizing the interaction between endocardial structures and blood flow and providing a potential model to include the effect of the complex structures within the heart without the need of segmenting an extremely complex structure and running large simulations every time. Therefore the absence of the mitral valve may constitute a limitation if we were drawing physiological conclusions, however it is not the case for the kind of study presented in this manuscript.

Given the different metrics analyzed, there appears to be no direct correlation between the volume of trabeculae and the intra-ventricular pressure drop. For the geometrical markers, subject D, for example, presented the highest 1Pdiff and is characterized by the smallest angle between the valvular planes (43.3◦ ). We hypothesize that the 1Pdiff obtained in model D is a result of the location of the PMs, which were positioned right below the inlet, disturbing the blood flow at the inlet, which led to a higher energy dissipation, not observable in the other cases. There is no direct correlation between the angle or distance of inlet and outlet and the 1Pdiff . It is clear that the existence of rugosities along the endocardial walls alter the hemodynamics by creating flow recirculation regions, vortex disruption into secondary vortices, which increase the energy dissipation, hence increasing the intra-ventricular pressure drop in complex ways. A key observation is that the location and orientation of the mitral and aortic valvular rings influence the direction of the flow, and hence, the high WSS visible either on trabeculated regions or on the PMs. We hypothesize that this observation is responsible for the high maximum WSS magnitude on the detailed geometries. Given the surface area of trabeculae or papillaries, WSS tends to be concentrated in small regions (**Figure 7**), increasing its maximum magnitude. This fact may have high implications on local tissue remodeling. However, regardless of the range of WSS, in **Figure 8** it can be observed that in the smoothed geometries the mode WSS was of 0.5–1 Pa, while in detailed meshes the mode drops to approximately 0 Pa. PMs and trabeculae reduced the WSS in LV about 23.5–66.7%. Wall shear stress is an important parameter in biology in general. Mechano-transduction is an important mechanism in biology. Even though it is practically impossible to measure it in-vivo in the LV, and it has never been reported before, the fact that WSS is reduced in the presence of trabeculae may provide an insight of the reason of why such endocardial structures exist. The presence of trabeculae and PMs generated a multitude of secondary vortices that were not present in the smoothed geometries, as shown in **Figure 9**. The overall vortex area decreased in the smoothed LVs (**Table 1**), with a mean of 31.6 ± 8.9 · 10−3 m2 . The presence of detailed endocardial structures increased the amount of vortices with a mean area of 43.5 ± 10.8 · 10−3 m2 . Subject D was intriguing. The effect of trabeculae and PMs led to a higher WSS median. This is because this LV is characterized by large PMs (noticeable in **Figures 1**, **8**), which led to higher WSS concentrated on the PMs. A 4.8% higher WSS was indeed observed in the detailed D case. We hypothesize that the high pressure drop in case D was due to the prominent PMs, which disturbed flow markedly, increasing the energy dissipation within the intra-ventricular volume, and thus, generating a high intra-ventricular pressure drop (see Figure S3). The analysis of more geometries is required to further understand and characterize the effect of endocardial structures on hemodynamics.

The results from the control subject E demonstrated that the presence of only the PMs led to a low 1Pdiff (0.2 kPa), however the main impact appears to be the WSS distributions in comparison with the more detailed geometries. In other words, having only PMs and a few big trabeculae does not significantly modify the WSS.

Using a transient inflows (E-A wave) allowed the study of vortex formations following physiological inputs. In the smoothed LV vortices were nominal and the generation of the vortex rings was clearly visible during the E waves (see the example in **Figure 10**). On the other hand, in the trabeculated ventricles, the vortex rings were disrupted, generating a multitude of secondary vortices. The vortex surface areas are provided in **Table 3**. It can be observed that the total surface areas of the vortices were larger in the trabeculated, in comparison to the smoothed geometry. Furthermore, in the anatomically detailed LV the secondary vortices penetrated deeper between the trabeculations during the late E wave (time instant 3). Notice that the A wave seems to produce higher vorticity in this model. We hypothesize that the reason for this is that we are starting our simulation with an organized zero flow all throughout the model. Recirculation within the cavity at zero flow (diastasis) would create higher vorticity in the second inflow wave. In the smoothed case vortices were more compact, while in the anatomically detailed LV there were multiple vortices that tended to be pushed toward the apex during the late A wave (time instant 6). This finding is in accordance with the results from the work of Vedula et al. (2016), in which the authors suggests that the observed behavior may help increasing LV washout. On the other side, the main vortex ring disruption and secondary vortices formation due to LV trabeculations was not seen in the work of Lantz et al. (2016), where they noticed that the presence of PMs and trabeculae generates a large vortex in the middle of the LV cavity.

Preliminary results from replacing the trabeculae with a porous layer show that it is possible to obtain an intra-ventricular pressure drop similar to those generated by the detailed endocardial models. Only for subject D, the intra-ventricular

TABLE 3 | Total vortex surface [m2 ] for the six instants (1−6) of the synthetic E-A wave in subject D.

pressure drop was not as high as the one observed in the detailed model. This is due to the presence of big PMs right below the mitral valve inlet, which, as explained previously, highly disturbed the flow and increased the energy dissipation.

Moreover, the addition of a porous layer on the smoothed geometries helped to reduce the WSS median in all the models providing small relative errors 0.07, 0.05, 0.23, and 0.16 respectively. Again, for model D, the median WSS had a bigger relative error mostly due to the large impinging of flow on the PMs, which increase the WSS median for that specific case. Finally, the presence of the porous layer disrupted the main vortices into smaller ones, increasing the total vortex surface to values slightly higher than the observed in the detailed cases, however, the relative errors range from −0.1 to −0.12, therefore reproducing the hemodynamic behavior in terms of vorticity.

#### 4.1. Limitations

The main limitation of this study is the use of CFD, without fluid-structure interaction (FSI), valves or moving walls. The lack of motion prevents us from comparing any measured value to in-vivo heart function, however, as was mentioned before, this study attempts to characterize solely the geometry effects on hemodynamics. Future work involves the use of FSI to compare our findings to in-vivo measurements. Another limitation is the small sample size, which limits the generalization and statistical significance of our results.

The generation of smoothed ventricles from the trabeculated ones provides a degree of variability in the geometries created. We removed the trabeculae keeping the overall shape of the ventricle and the volume unchanged. This however is observer dependent and requires a fair amount of user interaction. In models A–D, the smoothing procedure led to slightly smaller estimated LV volumes due to the elimination of the trabeculae. In the case of model E, the ventricular surface was primarily characterized by PMs, presenting just a few trabeculae, which led to a larger volume in the smoothed heart.

An important limitation is that valves were not considered in our simulations. The presence of the mitral valve will direct the blood flow jet to create impinging and this will have some influence on the interaction between the flow and the detailed endocardial structures. The valves have also been reported to create a vortex ring right below its leaflets, as observed in previous studies (Töger et al., 2012; Vedula et al., 2016), which is impossible to capture in this study.

Constant inflow simulations.

### 5. CONCLUSIONS

The highly anatomically detailed LV models developed in this study present a level of geometric information that was never achieved before. The simulations performed highlight the differences between blood flow CFD simulations in detailed vs. smoothed human ventricular models. The presence of detailed structures increase the intra-ventricular pressure drop, create multiple secondary vortices and decrease the WSS within the LV cavity. The amount of trabeculations have no direct correlation with the 1Pdiff , which was noted highest in the female LV D case. To the best of our knowledge, our study analyzed for the first time intra-ventricular pressure drops to investigate the effects of trabeculae and PMs on LV hemodynamic modeling. LV hemodynamics in detailed geometries are more complex than we anticipated, hence, a detailed study with more subjects is necessary and ongoing. Furthermore, our results confirm that neglecting detailed endocardial structures prevent computational models from recreating the complex blood flow behavior within the ventricles. Given that HPC simulations and high resolution MRI data are not always accessible, we propose that a simulated porous layer on the endocardial wall of smoothed LV models can potentially substitute the highly detailed geometries. Finally, we demonstrated that by adding a layer of 1.2 · 10−2 m thickness and 20 kg/m2 porosity to the smoothed cases we obtained pressure drops and WSS similar to the ones in the detailed LVs. The porous layer also increased the amount of secondary vortices, close to the amount observed inside the trabeculated models.

#### 6. RESOURCE IDENTIFICATION INITIATIVE

Atlas of Human Cardiac Anatomy, RRID:SCR\_015734.

MareNostrum Supercomputer, BSC, Barcelona, RRID:SCR\_015737. ARCHER, UK National Supercomputing Service, RRID:SRC\_015854. Fiji, RRID:SCR\_002285. Seg3D, RRID:SCR\_002552.

#### REFERENCES

ReMESH, RRID:SCR\_015735. ITK-SNAP Medical Image Segmentation Tool, RRID:SCR\_002010. Autodesk Meshmixer, RRID:SCR\_015736. Paraview, RRID:SCR\_002516.

### AUTHOR CONTRIBUTIONS

FS: Conception, drafting, data analysis and interpretation of data. BP: Data preprocessing, data analysis. TI, PI: Acquisition of data and critical revision. CB, JA-S: Conception, drafting, analysis, critical revision. MV, GH, OL: Implementation of the solvers, critical revision.

### FUNDING

This paper has been partially funded by CompBioMed project, under H2020-EU.<IP_ADDRESS> European Union's Horizon 2020 research and innovation programme, grant agreement n◦ 675451. FS is supported by a grant from Severo Ochoa (n◦ SEV-2015- 0493-16-4), Spain. CB is supported by a grant from the Fundació La Marató de TV3 (n◦ 20154031), Spain. TI and PI are supported by the Institute of Engineering in Medicine, USA, and the Lillehei Heart Institute, USA.

### ACKNOWLEDGMENTS

The DICOM datasets were provided by the Visible Heart R Laboratory, obtained by MRI scanning perfusion fixed hearts that were graciously donated by the organ donors and their families through LifeSource. Part of the simulation hours were provided by the CompBioMed project in the Archer supercomputer, EPCC, UK.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2018.00458/full#supplementary-material

could be the implications? Ann. Biomed. Eng. 44, 3346–3358. doi: 10.1007/s10439-016-1614-6

Research, Proceedings of the Summer Program, (Stanford, CA), 193–208.

fractional flow reserve: scientific basis. J. Am. College Cardiol. 61, 2233–2241. doi: 10.1016/j.jacc.2012.11.083

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Sacco, Paun, Lehmkuhl, Iles, Iaizzo, Houzeaux, Vázquez, Butakoff and Aguado-Sierra. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Scalable and Accurate ECG Simulation for Reaction-Diffusion Models of the Human Heart

#### Mark Potse1,2,3 \*

1 CARMEN Research Team, Inria Bordeaux Sud-Ouest, Talence, France, 2 Institut de Mathématiques de Bordeaux , UMR 5251, Université de Bordeaux, Talence, France, 3 IHU Liryc, Electrophysiology and Heart Modeling Institute, Foundation Bordeaux Université, Pessac-Bordeaux, France

Realistic electrocardiogram (ECG) simulation with numerical models is important for research linking cellular and molecular physiology to clinically observable signals, and crucial for patient tailoring of numerical heart models. However, ECG simulation with a realistic torso model is computationally much harder than simulation of cardiac activity itself, so that many studies with sophisticated heart models have resorted to crude approximations of the ECG. This paper shows how the classical concept of electrocardiographic lead fields can be used for an ECG simulation method that matches the realism of modern heart models. The accuracy and resource requirements were compared to those of a full-torso solution for the potential and scaling was tested up to 14,336 cores with a heart model consisting of 11 million nodes. Reference ECGs were computed on a 3.3 billion-node heart-torso mesh at 0.2 mm resolution. The results show that the lead-field method is more efficient than a full-torso solution when the number of simulated samples is larger than the number of computed ECG leads. While the initial computation of the lead fields remains a hard and poorly scalable problem, the ECG computation itself scales almost perfectly and, even for several hundreds of ECG leads, takes much less time than the underlying simulation of cardiac activity.

#### Edited by:

Mariano Vázquez, Barcelona Supercomputing Center, Spain

#### Reviewed by:

Arun V. Holden, University of Leeds, United Kingdom Mohammad Hasan Imam, American International University-Bangladesh, Bangladesh

#### \*Correspondence:

Mark Potse<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology

Received: 14 January 2018 Accepted: 27 March 2018 Published: 20 April 2018

#### Citation:

Potse M (2018) Scalable and Accurate ECG Simulation for Reaction-Diffusion Models of the Human Heart. Front. Physiol. 9:370. doi: 10.3389/fphys.2018.00370 Keywords: numerical modeling, electrocardiogram, high-performance computing, reaction-diffusion model, bidomain model, lead fields

### 1. INTRODUCTION

The electrocardiogram (ECG) is one of the most common tools in present-day medicine, yet its relation with the molecular biology of the heart is still poorly understood. The ECG witnesses the collective activity of about a million current-generating transmembrane proteins in each of the heart's muscle cells (Hille, 2001). Many of these proteins have been identified and their actions have been captured in mathematical models that predict their collective behavior on the scale of a cell (Noble and Rudy, 2001). By coupling millions of these membrane models one can create a model of whole-heart electrophysiology. Such models generate crucial insights in the functional effects of molecular-level changes, allowing for example to predict dangerous side effects of new drug designs (Passini et al., 2017) or to understand how cardiac ion-channel mutations influence cardiac rhythm disorders (Gima and Rudy, 2002). Moreover, from their results one can compute the corresponding ECG and predict how lab results on subcellular components would translate to everyday practice (Hoogendijk et al., 2010; Keller et al., 2012; Zemzemi et al., 2013).

Such realistic models are large and, when run on a single processor, would take days to simulate just one heartbeat. Fortunately the problem can be expressed in such a way that the work may be spread over many processors with little communication between them. Therefore, these computations are said to scale very well, meaning that they run almost twice as fast every time the number of processors is doubled (Vázquez et al., 2011). This makes them suitable for use on large-scale parallel computers, allowing models to run in nearly real time (Niederer et al., 2011b; Richards et al., 2013).

Simulation of a realistic ECG from the results of such a numerical heart model is much harder, because the electrical current generated by the heart meets a different conductivity at each point in the torso. As a result, each point influences the potential everywhere else, so to find the potential anywhere one must solve it everywhere at the same time.

Numerically this means that a large system of linear equations must be solved, one for each point in the torso model. These problems are harder when they are larger and require frequent communication between the processors in a parallel computer. This means that they cannot be solved much faster by using more processors. Therefore, ECG computation is becoming a bottleneck, limiting both the speed and the spatial resolution of our models.

To avoid this problem many researchers have used simplified torso models, resulting in a less accurate ECG. A solution that can avoid such a sacrifice is to simulate the ECG using an electrocardiographic concept named a lead field. This allows the problem to be split into a hard (poorly scaling) part and an easy (well scaling) part. The hard part is solved only once for each ECG lead, while the easy part is run repeatedly for each time step in a simulation and for multiple simulations on the same geometry. This approach has been used by several authors, but generally with simplified heart models (Pezzuto et al., 2017) or, again, with simplified torso models (Horacek, 1973; Miller and Geselowitz, 1978; Mailloux and Gulrajani, 1982; Aoki et al., 1987).

The purpose of this paper is to show that a lead-field approach can greatly improve scalability in a high-performance computing (HPC) context without sacrificing accuracy. This is not obvious, because the method requires a large set of transfer coefficients (the lead field) to be stored between the two phases of the computation. The efficiency of the method depends on the accuracy with which the lead field must be computed and the degree to which it can be downsampled without affecting the accuracy of the ECG too much. Finally, to provide answers to these questions an accurate reference solution is needed.

Using a reference solution computed on a full torso model at 0.2 mm resolution this study shows that the lead field can indeed be downsampled enough to achieve an efficient and scalable computation, providing roughly two orders of magnitude speedup with negligible loss in accuracy.

The results of this study make it possible to build more realistic heart models with higher spatial resolution, without spending much more time to compute the ECG.

### 2. METHODS

#### 2.1. Model Equations

The methods in this study are based on the bidomain model of cardiac electrophysiology (Miller and Geselowitz, 1978; Tung, 1978), on which most of the current modeling work in this area is based (Niederer et al., 2011a; Henriquez, 2014). The bidomain model is a continuum approximation of the heart muscle, which in reality consists of a network of interconnected muscle cells embedded in an extracellular matrix and other structures such as fibroblasts and capillaries. The bidomain model approximates this as two co-located spaces: the intracellular domain, consisting of the interior of the cells and the gap junctions that connect them, and the extracellular domain, consisting of everything else.

The two domains are characterized by conductivity tensors Gi and Ge, respectively. Their values at each point in the model depend on the fiber direction and account for the partial volume occupation of the two domains. In addition the parameters Cm and β determine the capacitance of the cell membrane and the amount of membrane per unit volume, respectively. The state variables of the model are the potential fields φi in the intracellular and φe in the extracellular domain, and a set of variables Ey describing the state of the membrane model at each location. Using the auxiliary variable Vm = φi − φe and agreeing that all variables are functions of time and position we can express the bidomain model compactly as

$$\boldsymbol{\beta}^{-1}\nabla \cdot (\mathbf{G}\_{\mathrm{i}}\nabla\phi\_{\mathrm{i}}) = \mathbf{C}\_{\mathrm{m}}\partial\_{\mathrm{i}}V\_{\mathrm{m}} + I\_{\mathrm{ion}}(V\_{\mathrm{m}}, \vec{\jmath})\tag{1}$$

$$
\beta^{-1} \nabla \cdot (\mathbf{G\_e} \nabla \phi\_e) = -C\_{\rm m} \partial\_t V\_{\rm m} - I\_{\rm ion} (V\_{\rm m}, \vec{y}) \tag{2}
$$

$$
\partial\_t \vec{\jmath} = F(V\_{\mathfrak{m}}, \vec{\jmath}) \tag{3}
$$

where the term Cm∂tVm represents the capacitive transmembrane current, the function Iion the density of ionic current flowing between the two domains, and F is a nonlinear vector-valued function describing how the membrane state evolves. The pair of functions Iion and F constitutes the membrane model. Suitable boundary conditions are

$$G\_{\mathbf{i}} \nabla \phi\_{\mathbf{i}} \cdot \partial \Omega\_{\mathbf{A}} = 0 \tag{4}$$

on the boundary A of the cardiac muscle and

$$G\_{\mathbf{e}} \nabla \phi\_{\mathbf{e}} \cdot \partial \Omega\_{\Gamma} = 0 \tag{5}$$

on the torso boundary T (Tung, 1978; Krassowska and Neu, 1994).

The electrical activity of the heart can then be simulated by integrating Equations (1), (2), and (3) under the boundary conditions (4) and (5) (Vigmond et al., 2002). This is known as a bidomain reaction-diffusion model. In this study a simplified version, a "monodomain" reaction-diffusion model, was used. This model can be derived by assuming that Gi and Ge are proportional (Leon and Horácek, 1991). Although this is a gross simplification the effect of this assumption is negligible for most purposes if the model parameters are well chosen (Potse et al., 2006; Nielsen et al., 2007; Bishop and Plank, 2011; Coudière et al., 2014). The monodomain model reads

$$\begin{cases} \mathcal{C}\_{\rm m} \partial\_{t} V\_{\rm m} = \beta^{-1} \nabla \cdot \left( G\_{\rm m} \nabla V\_{\rm m} \right) - I\_{\rm ion} (V\_{\rm m}, \vec{\mathcal{y}}) \\ \partial\_{t} \vec{\mathcal{y}} = F(V\_{\rm m}, \vec{\mathcal{y}}) \end{cases} \tag{6}$$

The "monodomain conductivity tensor" Gm was computed as the series conductivity of the two domains, Gm = GiGe/(Gi + Ge). With this choice the resistance encountered by a current loop through the cell membrane is the same as in a bidomain model, so that also the conduction velocity of a propagating activation wavefront is almost the same.

An ECG potential V(t) at time t is the difference in φe between two locations on the body surface or, more generally, a linear combination

$$V(t) = \sum\_{i} c\_{i} \phi\_{\mathbf{e}}^{i} \tag{7}$$

where ci are the relative contributions of the two or more electrodes and φ i e are the potentials at the corresponding positions. The coefficients P ci must fulfill charge conservation, ci = 0.

To compute φe we must return to the bidomain model. Equations (1) and (2) can be combined and reorganized to yield

$$\nabla \cdot \left( \left( G\_{\text{i}} + G\_{\text{e}} \right) \nabla \phi\_{\text{e}} \right) = -\nabla \cdot \left( G\_{\text{i}} \nabla V\_{\text{m}} \right). \tag{8}$$

This equation can be solved for φe in the whole torso at once from a given distribution of Vm. However, for the ECG we need to know φe at a few locations only. Therefore, it can be more efficient to use a Green's function of the operator ∇ · ((Gi + Ge)∇.) for each of these locations. Since an ECG lead is a linear combination of φe at two or more points it can also be represented directly by a linear combination of Green's functions. In electrocardiology such linear combinations of Green's functions are named lead fields (McFee and Johnston, 1953; Geselowitz, 1989; Colli-Franzone et al., 2000). A lead field is computed once for each ECG lead. It is then used to evaluate the ECG at each time step of the reaction-diffusion model and, as long as the conductivity parameters are not changed, can be re-used for multiple simulations. In terms of a lead field Z(xE) the ECG potential V(t) at time t is

$$V(t) = \int \nabla Z(\vec{\chi}) \cdot G\_{\rm i} \nabla V\_{\rm m} \, d\vec{\chi} \tag{9}$$

where the integration is over the myocardium. In contrast to the solution of the full system (8) this calculation is simple and a priori highly scalable. The lead field can be computed as the potential field resulting from a unit current applied at the electrode locations xEi (Geselowitz, 1989):

$$\nabla \cdot \left( \left( G\_{\vec{\mathbf{x}}} + G\_{\mathbf{e}} \right) \nabla Z(\vec{\mathbf{x}}) \right) = \sum\_{i} c\_{i} \delta(\vec{\mathbf{x}} - \vec{\mathbf{x}}\_{i}) \tag{10}$$

where the coefficients ci are as in Equation (7) and δ is Dirac's delta function. To avoid a scaling factor in (9) the total injected current must be unitary, P|ci | = 2.

#### 2.2. Model Geometry

In order to run tests on a relevant geometry a model of the heart and torso was used that had been created for a previous study (Kania et al., 2017). The methods to build this geometry, only tersely described before, were as follows. High-resolution cardiac and thoracic computed tomography (CT) images were obtained from a female patient in her thirties. Images were segmented automatically using the MUSIC software (IHU Liryc, Université de Bordeaux and Inria Sophia Antipolis, France), under supervision of an expert operator. The boundaries of the segmented volumes were expressed as triangulated surfaces and meshing errors were manually corrected using Blender (The Blender Foundation, Amsterdam, The Netherlands). The resulting surface mesh defined the volumes of the ventricular myocardium, left and right cavities with parts of the great vessels, lungs, and the whole body. To define hexahedral meshes for the computations the surfaces were overlaid with a 3D cartesian mesh whose elements were assigned types according to the surfaces in which they were contained. The bones were also segmented and meshed but not included in the simulations. The atrial myocardium was not segmented.

The heart mesh was processed to define subendocardial and subepicardial layers and fiber directions using the rule proposed by Beyar and Sideman (1984), as previously described (Potse et al., 2006). The torso mesh was similarly processed to define a layer of 1 cm thickness directly under the skin as skeletal muscle and to define a sheet direction in this layer. Since the true fiber directions of the skeletal muscle layer are too complex to account for the model muscle simply had a low conductivity in the radial direction and a high conductivity in all circumferential directions (**Table 1**).

During the thoracic scan the patient was wearing a vest with 252 embedded electrodes (Tilt et al., 2013; Cochet et al., 2014). The locations of these electrodes were extracted from the CT data using software provided by the manufacturer of the vest. In addition the locations of the 9 standard ECG electrodes were determined by referring to the bone mesh, and two electrode locations on the hips were chosen. The surface mesh with electrode positions is illustrated in **Figure 1**.

#### 2.3. Spatial Discretization

Spatial discretization was done using a finite-difference method. Differential operators of the form ∇ · (G∇.), where G is any of the conductivity tensor fields employed, were computed using an

TABLE 1 | Tissues used in the simulations together with the volumes they occupy in the torso model, the conductivity parameters σ (in mS/cm), and β (cm−1 ); the subscript "i" stands for intracellular, "e" for extracellular, "L" for longitudinal, "T" for transverse (within a tissue sheet), and "C" for across-sheet.

expression proposed by Saleheen and Ng (1997). This expression assumes that G is constant on elements and that potentials are defined on the nodes of the mesh. It produces a 19-point stencil that takes anisotropy and inhomogeneities into account. The simulation code read its geometry in terms of elements, and created a node mesh, assigning node types such that all corners of a myocardial element would have myocardial nodes. In order to treat myocardial boundaries correctly, the β value of each node was the average of those associated with the 8 elements around it, which was zero for non-myocardium (Potse et al., 2006).

#### 2.4. Simulation of Cardiac Activity

To prepare input data for ECG simulation propagating activation was simulated using the monodomain reaction-diffusion model (6) using the membrane model of Ten Tusscher and Panfilov (2006) for the functions F and Iion. A uniform time step of 10 µs was used. At each time step the code

After each 100 time steps results were written to file. Simulations were run on a heart mesh at 0.2 mm resolution. Tissue parameters determining Gm and β are listed in **Table 1**. Gating variables were integrated with the method of Rush and Larsen (1978) and all other variables with a forward Euler method.

Activation was started with a single stimulus at one location, at the beginning of the simulation. Seven simulations were run, each time with the stimulus at a different location. Simulations covered 500 ms to include the full depolarization and repolarization of the ventricles.

### 2.5. ECG Simulation

The ECG was computed with several methods:

$$\nabla \cdot \left( \left( G\_{\text{i}} + G\_{\text{e}} \right) \nabla \phi\_{\text{e}} \right) = -I\_{\text{w}} \tag{11}$$

where Iw is a projection of the term ∇ · (Gi∇Vm) from a 0.2 mm resolution heart mesh onto a 1 mm resolution torso mesh. Each coarse-mesh node received contributions from a cube-shaped area including all fine-mesh nodes within the up to 8 coarse-mesh elements around it, with higher weights attributed to nearby nodes, as in a trilinear interpolation: Let 1x, 1y, 1z be the number of fine-mesh edges between a coarse-mesh node and a fine-mesh node along the x, y, and z axis, respectively. Then the contribution of the fine-mesh node to the coarse-mesh node was

$$\mathcal{W} = \begin{cases} 0, & \text{if } \Delta x \ge 5 \lor \Delta y \ge 5 \lor \Delta z \ge 5\\ (5 - \Delta x)(5 - \Delta y)(5 - \Delta z)/5^6, & \text{otherwise} \end{cases}$$

The coarse mesh was constructed such that a myocardial fine-mesh node was always surrounded by 8 coarse-mesh nodes. Therefore, w added up to unity for each fine-mesh node and charge conservation was ensured.

For the FSC method the monodomain reaction-diffusion model (6) was integrated in a separate run which saved Iw to file. This method has been used routinely in several studies (Nguyên et al., 2015; Meijborg et al., 2016; Duchateau et al., 2017; Kania et al., 2017). The torso mesh in this case consisted of 2.7 · 107 nodes.

The notations LF(C, S) and LFS(C, S) will be used for the LF and LFS methods, respectively, with lead fields computed at a resolution of C millimeters and downsampled to a resolution of S millimeters.

#### 2.6. Computation of Lead Fields