Mechanical Characterization of the Vessel Wall by Data Assimilation of Intravascular Ultrasound Studies

Atherosclerotic plaque rupture and erosion are the most important mechanisms underlying the sudden plaque growth, responsible for acute coronary syndromes and even fatal cardiac events. Advances in the understanding of the culprit plaque structure and composition are already reported in the literature, however, there is still much work to be done toward in-vivo plaque visualization and mechanical characterization to assess plaque stability, patient risk, diagnosis and treatment prognosis. In this work, a methodology for the mechanical characterization of the vessel wall plaque and tissues is proposed based on the combination of intravascular ultrasound (IVUS) imaging processing, data assimilation and continuum mechanics models within a high performance computing (HPC) environment. Initially, the IVUS study is gated to obtain volumes of image sequences corresponding to the vessel of interest at different cardiac phases. These sequences are registered against the sequence of the end-diastolic phase to remove transversal and longitudinal rigid motions prescribed by the moving environment due to the heartbeat. Then, optical flow between the image sequences is computed to obtain the displacement fields of the vessel (each associated to a certain pressure level). The obtained displacement fields are regarded as observations within a data assimilation paradigm, which aims to estimate the material parameters of the tissues within the vessel wall. Specifically, a reduced order unscented Kalman filter is employed, endowed with a forward operator which amounts to address the solution of a hyperelastic solid mechanics model in the finite strain regime taking into account the axially stretched state of the vessel, as well as the effect of internal and external forces acting on the arterial wall. Due to the computational burden, a HPC approach is mandatory. Hence, the data assimilation and computational solid mechanics computations are parallelized at three levels: (i) a Kalman filter level; (ii) a cardiac phase level; and (iii) a mesh partitioning level. To illustrate the capabilities of this novel methodology toward the in-vivo analysis of patient-specific vessel constituents, mechanical material parameters are estimated using in-silico and in-vivo data retrieved from IVUS studies. Limitations and potentials of this approach are exposed and discussed.

1. Introduction

Cardiovascular diseases are the principal cause of death and morbidity worldwide (Mathers et al., 2016). The two principal causes of death, cardiac ischemia and stroke, are intrinsically related with the onset and progress and destabilization processes of atherosclerotic plaque, which are still largely unknown (Crea and Liuzzo, 2013; Bentzon et al., 2014). At the final stage of the destabilization process, the plaque ruptures releasing thrombotic components into the blood stream which in turn generate thrombi that block the vessel lumen causing ischemia. Thus, the prediction of rupture events and the identification of the so-called culprit plaques is of the utmost importance for diagnostics and therapeutics. Through computational simulations is it possible to study the arterial wall stress state, which may compromise plaque integrity and induce rupture. Moreover, computational models also allow to recreate different physiological and pathophysiological conditions (hypertension, hyperemia, exercise, stenosis) (Taylor et al., 1999, 2013; Torii et al., 2007; Blanco et al., 2015), as well as interventions (angioplasty balloon inflation, stent deployment, stent-plaque interaction, among others) (Conway et al., 2012, 2014) that are valuable resources for diagnosis, treatment and surgical risk assessment.

In order to accurately simulate patient specific conditions, three kinds of input data are required: (i) patient-specific anatomical models of the vasculature, (ii) the loads to which the anatomical structures are subjected to, and (iii) the patient-specific distribution of the arterial-wall constituents and their corresponding material parameters. As far as anatomical data of the arteries is concerned, it can be straighforwardly extracted from different medical imaging modalities (Wahle et al., 1995; Milner et al., 1998; Bulant et al., 2017). Regarding the force exerted by the blood pressure, it can be accurately estimated from cuff-pressure measurements (O'brien et al., 2001; Miyashita, 2012). Thus, we are left to the problem of setting patient-specific material parameters for the models of the arterial wall. This has long been the Achilles tendon in numerical simulations, most of them relying in material parameters acquired from ex-vivo material experimentation in cadaveric specimens (Walsh et al., 2014; Karimi et al., 2015). In this sense, the in-vivo identification of material parameters for the arterial-wall is still an open research topic.

Toward covering the aforementioned gap, specifically in the coronary artery disease domain, intravascular ultrasound (IVUS) emerges as an suitable imaging modality to make the attempt to retrieve the material parameters and distribution of the vessel materials under in-vivo conditions due to its high temporal and spatial resolution. The acquired images, when coherently ordered, are capable of delivering the motion of the vascular structures. Some works (Kawasaki et al., 2002; Nair et al., 2002; Sathyanarayana et al., 2009) have successfully classified the materials in few discrete categories (e.g., necrotic core, fibrotic, fibro-fatty or lipid-pool, calcified) based on the acoustic impedance response of the tissues in a determined frame of the IVUS study. It has then been demonstrated that there is a notorious variability of the stress-strain response of tissues within the same category (Loree et al., 1994; Holzapfel et al., 2005; Walsh et al., 2014) of such classification. Therefore, this information is not specific enough for simulation purposes. As anticipated above, the temporal resolution of the IVUS study can be exploited to retrieve the motion (displacement field) of the vessel wall along the cardiac cycle (for example by using optical flow techniques or large deformation diffeomorphic metric mapping). Using the displacement field as input, data assimilation techniques can be supplied to estimate the material parameters.

Data assimilation techniques make use of measurable quantities to adjust a physical model whose goal is to represent the reality posed by the in-vivo scenario. In that manner, these techniques permit not only to estimate specific quantities of interest, but also to explore the underlying physical phenomena. Also, measurement errors can be filtered by the physical model being a quid pro quo benefit: the measurements instantiate the model and the model filters the measurements. Such techniques can be classified in two categories: (i) variational approaches and (ii) sequential filtering approaches.

In the variational approach, a cost functional that measures the difference between the observed measures and the model prediction is constructed. The cost functional depends on the parameters of interest (among other parameters required by the model) to render a model prediction of the measured variable. Then, the estimated parameters are those such that minimize the cost functional. The more popular approach is to solve the Karush-Kuhn-Tucker (KKT) necessary conditions which is employed in several works for mechanical parameter estimation (Lagrée, 2000; Martin et al., 2005; Sermesant et al., 2006; Perego et al., 2011; D'Elia et al., 2012; Bertagna and Veneziani, 2014; Ares, 2016). In Lagrée (2000), the viscoelastic parameters of large arteries were estimated using displacement fields of the vessel wall generated by computational models. Similarly, Martin et al. (2005) explored the estimation of the vessel compliance in a 1D model using a 3D fluid-structure interaction (FSI) model to generate the measured displacement of the vessel wall. Using medical data of blood pressure and inner radius of the arteries, Stålhand (2009) also used 1D models to estimate the material parameters according to the model proposed in Holzapfel et al. (2000). The works of Perego et al. (2011) and D'Elia et al. (2012) formulate the inverse problem from 3D FSI models and analyze the sensitivity in the identification of Young modulus to noise in the measurements of arterial wall displacements. In the latter, data assimilation is performed from flow velocity data as well. The main drawback of these variational approaches is the large number of evaluations of the cost functional (or its derivative) which are required in the minimization problem (Lassila et al., 2013). Furthermore, the use of more realistic models such as 3D FSI models or complex heterogeneous anisotropic solid models are many times mandatory to render accurate results, increasing the computational effort. In some cases, reduced order strategies combined with statistical approaches can be applied to reduce the burden behind cost functional evaluations, as shown in Lassila et al. (2013). Other approach is proposed in Bertagna and Veneziani (2014), based on the application of model reduction techniques coupled with a proper orthogonal decomposition to accomplish the solution of 3D FSI in a computationally efficient way. Efficient implementations for solid mechanics problems have also been proposed in Avril et al. (2010) and Pérez Zerpa and Canelas (2016) using a virtual fields method and a constitutive equation gap functional, respectively.

In turn, and for problems involving a small-to-moderate number of unknown parameters, the sequential filtering approach (also known as filtering methods) is less computationally demanding and, at the same time, embarrassingly parallel. These features make the filtering approach an appealing strategy for the kind of problems addressed in the present work. Conceptually, given a set of observations, the method realizes a prediction for each observation and, then, introduces corrections in the model parameters based on the discrepancies between the model estimation and the observed data. For each prediction-correction step, several variations of the parameters are tested in the model and, through statistical analysis of the model predictions, a suitable correction is performed over the parameters. Several methods based on the Kalman filter have been developed to deal with linear and non-linear dynamic problems. As examples, a non-linear extended Kalman filter (EKF) with collocation feedback is applied to identify the Young modulus of different regions of a heart model in Moireau et al. (2008), Moireau et al. (2009), and Chapelle et al. (2009). The observations used varied between the myocardium velocity (Moireau et al., 2008), displacement (Moireau et al., 2009) and velocity of the heart boundaries (Chapelle et al., 2009). The stability of such methods was studied (Moireau et al., 2008) and in terms of accuracy it is reported that Kalman filtering is optimal for linear systems only, while extended algorithms based on linearized operators may lead to efficient, albeit non-optimal, filtering procedures. In Lipponen et al. (2010), the EKF is also applied to estimate parameters of a reduced order Navier-Stokes model (through an orthogonal decomposition of the velocity field) through observations acquired from electrical impedance tomography. In more recent works, Moireau and Chapelle (2011) presented a reduced order Kalman filter based on the unscented transform (abbreviated as ROUKF) that offers an interesting alternative to the EKF method. Such an approach does not require neither linearization nor calculation of the tangent operator of the non-linear model, which substantially eases its implementation. Noteworthy, the ROUKF features a higher order approximation of the system states statistics, delivering more accurate outcomes than EKF. In Bertoglio et al. (2012) and Bertoglio et al. (2014), ROUKF was successfully applied for estimation of Young modulus in arteries with tests in-vivo and in-vitro, showing a simpler and more efficient implementation than EKF. Recently in Caiazzo et al. (2017), terminal resistances and vessel wall properties of a 1D vascular network were estimated via ROUKF using blood flow and/or pressure measurements as observations.

In this work, we present a novel approach to construct patient-specific mechanical models of the arterial wall using in-vivo data from IVUS studies. In a nutshell, this approach integrates the realms of image processing, optical flow, continuum mechanics, and filtering data assimilation to effectively merge patient-specific data with mechanical models, toward the in-vivo estimation of material properties. From the IVUS study, a frame of interest is selected and the corresponding arterial wall is demarcated. For the mechanical model a finite strain framework is considered, and the constituent tissues are assumed to behave as isotropic Neo-Hookean materials. Importantly, it is considered that the arterial vessel corresponding to the diastolic phase is at equilibrium with a certain diastolic pressure level, and it is further subjected to a given axial stretch at that phase. By using gating, registration and optical flow methods developed in previous works (Maso Talou et al., 2015, 2017; Maso Talou, 2017), the displacement field of the vessel wall is estimated along the cardiac cycle. Then, the ROUKF is exploited as a data assimilation procedure in which the previously obtained displacement field is supplied as observational data, while the material parameter of the Neo-Hookean models are the target parameters to be estimated.

The manuscript is structured as follows. In section 2, the proposed methodology is detailed, presenting image processing techniques (section 2.1), the mechanical model for the arterial wall (section 2.2), and, at last, the data assimilation process for the estimation of the material parameters (section 2.3). In section 3, the sensitivity of the data assimilation parameters (section 3.1) and boundary conditions (section 3.2) and baseline stress state (section 3.3) for the mechanical problems are studied to assess their impact on the data assimilation outcomes. Hence, the mechanical characterization is performed for four in-vivo atherosclerotic lesions to analyze the performance of the method in real case scenarios (section 3.4). Insights, strengths and weaknesses of the methodology are then discussed in section 4 and final remarks are outlined in section 5.

2. Methods

This section is divided in four parts. First, the IVUS imaging processing methods are described, where we present the procedures to obtain the displacement field of a specific vessel cross-section along the cardiac cycle (see Figure 1). Second, the mathematical model for the arterial mechanics is formulated, defining the mechanical equilibrium and the material constitutive behavior. Third, the data assimilation algorithm is presented as a tool to estimate unknown material properties in the mechanical models using the displacement field retrieved from the IVUS images. Finally, an efficient three-level parallelization scheme is described for high performance computing environments.

FIGURE 1

Figure 1. Proposed pipeline to estimate the patient-specific mechanical properties from IVUS medical images.

2.1. Image Processing

The goal of the image processing stage is to deliver the displacement field of the vessel wall along the cardiac cycle at a particular site of interest within the artery. As this new methodology is a proof of concept, the data from the in-vivo cases will be extracted from a standard IVUS pullback as a retrospective study. As the transducer is axially displaced from frame to frame, only images corresponding to a single cardiac cycle can be extracted for each cross-section to obtain small topological variations between the images (spatial consistency). Hence, the extraction of the frames at a particular location is hindered due to the motion of the IVUS transducer exerted by the myocardium contraction. To overcome this issue, gating and registration procedures are performed using methods previously presented in Maso Talou et al. (2015, 2017). To retrieve the displacement field, a modified optical flow method is applied to the extracted frames at the site of interest. As follows, the treatment given to the IVUS images is briefly described.

2.1.1. Gating

The gating method aims to recover the cardiac phase at each cross-sectional image of the study. To achieve this, a signal that measures the total motion of each frame is generated as

$\begin{array}{l}s(n)={\alpha }_g\left[1-\frac{{\sum }_{i=1}^H{\sum }_{j=1}^W\left(I_n(i,j)-{\mu }_n\right)\left(I_{n+1}(i,j)-{\mu }_{n+1}\right)}{{\sigma }_n{\sigma }_{n+1}}\right]\\ +(1-{\alpha }_g){\displaystyle \sum _{i=1}^H{\displaystyle \sum _{j=1}^W-\left|\nabla I_n(i,j)\right|}},& (1)\end{array}$

where I_n is the n-th image of the study with a resolution of H × W pixels, μ_n and σ_n are the mean and standard deviation of the intensity at I_n and α_g a mixture parameter. The principal frequency mode of the signal s(n) at the physiological heart-frequency range (i.e., between 0.75 and 1.66 Hz) is extracted to obtain the mean cardiac frequency of the study, f_m. Then, a low frequency signal s_l(n) is generated by low-pass filtering s(n) with cut-frequency f_c = 1.4f_m. If there is not severe arrhythmia during the IVUS acquisition, s_l presents one minimum per cardiac cycle related to the end-diastolic phase, thus, all frames for this phase are easily and directly extracted. Due to heartbeat period variability along the study, some of these minima can be displaced between s and s_l, because of the lack of high frequencies contributions. To avoid such inconsistencies, we iteratively modify $f_c^k=(k+0.4)f_m$ (k is the current iteration number), recompute $s_l^k$ with the new cut-frequency $f_c^k$ and adjust each minimum of iteration k−1 to its nearest local minimum in $s_l^k$. Interestingly, the iterative scheme aids in cases with mild arrhythmia, i.e., where only few heartbeats of the study (not contiguous) present delay or omission of the P-wave. In those cases, the adjustment of the minima identified correctly the P-waves or collapsed the two minima to the same time position (this is the case when a P-wave did not occur and the heartbeat elapsed twice its period). In both of the previous cases, the minima are correct. In cases with severe arrhythmia, it is recommended the use of ECG signal and manual segmentation of minima for a proper gating. In the so-obtained phase, the cardiac contraction is at its minimum, and so, it corresponds to the beginning of the cardiac cycle, more precisely the beginning of the cardiac P-wave.

Since the heart frequency changes along the study, the heartbeats are sampled with a variable amount of frames. This variability in the heartbeat frequency affects mainly the relaxation process of the heart and, consequently, the length of the T-P interval. Despite this, the P-T interval remains almost invariant. Taking this fact into consideration, the end-diastolic instant for each cardiac cycle will be regarded as a reference for the definition of S cardiac phases. Each available frame of the P-T interval is then associated to a specific cardiac phase, obtaining phase-coherent volume datasets. Further details of the gating method, setup of the mixture parameter α_g and validation with in-vivo studies are described in Maso Talou et al. (2015).

2.1.2. Registration

All phase-coherent volumes are registered (axially and transversally) against the volume dataset corresponding to the end-diastolic phase. This procedure is performed for each phase-coherent volume. The transversal registration is achieved by finding the in-plane rigid motion for each image in the current phase that best matches the frame image in the end-diastolic phase. To quantify the matching between two images, we use a maximum likelihood estimator presented in Cohen and Dinstein (2002) and Wachinger et al. (2008),

$\begin{array}{ll}c(I_n,I_m)={\displaystyle \sum _{i=1}^H{\displaystyle \sum _{j=1}^W[I_n(i,j)-I_m(i,j)-\log (e^{2(I_n(i,j)-I_m(i,j))}+1)].}}& (2)\end{array}$

The rigid motion Ξ_n for each cross-section is then estimated by solving the following optimization problem

$\begin{array}{ll}{\Xi }_n=\underset{{\Xi }^{*}}{\text{arg}\text{max }}c\left(I_n^{\text{D}},I_n^s(x({\Xi }^{*}),y({\Xi }^{*}))\right),& (3)\end{array}$

where $I_n^s$ is the n-th cross-section of the phase-coherent volume corresponding to the s-th phase, D denotes the end-diastolic phase, and I(x(Ξ*), y(Ξ*)) is the image I after applying the rigid transformation defined by Ξ* which is composed by an in-plane translation plus a rotation with respect to the image center.

By virtue of the myocardium contraction, the same cross-sections site at the different phases may be longitudinally displaced. Therefore, it is necessary to perform an axial registration to find the corresponding frames at different phases for the same transversal site. Thus, after transversal registration of all phase-coherent volumes, an axial registration against the end-diastolic phase is applied. For each frame of each phase-coherent volume (now transversally registered), the best matching frame in the end-diastolic volume is sought out. To diminish the computational burden, the search is limited to the 14 adjacent frames in the end-diastolic volume which is within the range of axial displacements of a transducer during the IVUS study (Arbab-Zadeh et al., 1999). To quantify the matching between two images, we use a neighborhood likelihood estimator defined as

$\begin{array}{ll}{c}_w(I_n^s,I_m^{\text{D}})=\frac{{\displaystyle \sum _{d=-w}^w}{\varphi }_{{\sigma }_{\text{G}}}(d)c(I_{n+d}^s,I_{m+d}^{\text{D}})}{{\displaystyle \sum _{d=-w}^w}{\varphi }_{{\sigma }_{\text{G}}}(d)},& (4)\end{array}$

where ϕ_σG is a Gaussian weight function with σ_G standard deviation and w is the amount of adjacent frames used to establish the matching between the two sites centered at $I_n^s$ and $I_m^{\text{D}}$ respectively. It is important to note that w is not the search range fixed at 14 frames, but is the size of the neighborhood used for each comparison between two frames. Then, the position for axial registration, i. e., frame of the end-diastolic phase that best matches the current frame $I_n^s$ is given by

$\begin{array}{ll}m={\displaystyle \underset{k=n-7,\dots ,n+7}{\text{arg max}}}{c}_w(I_n^s,I_k^{\text{D}}).& (5)\end{array}$

Finally, given the site of interest at the n-th frame of the end-diastolic phase volume, the set of frames that constitutes a sequence along the cardiac cycle at this site is ${I}=\left\{{Ĩ}_n^s,s=1,\dots ,S\right\}$, where ${Ĩ}_n^s$ is the n-th frame of the phase-coherent volume corresponding to phase s after transversal and axial registration. The reader is directed to Maso Talou et al. (2017) and Maso Talou (2017) for further details of the registration methods.

2.1.3. Optical Flow

For a pair (or sequence) of images, optical flow techniques aim at determining the displacement vector field that relates the points of both images (Horn and Schunck, 1981). Because optical flow strategies rely on the gray constancy assumption, a denoising procedure is performed over the sequence. The applied denoising method is a variation of the TV-L1 method (Rudin et al., 1992; Chan et al., 1999) which modifies the data term (absolute difference measurement of the images) by the negative maximum likelihood estimator assuming one image with gamma distributed noise and another noiseless image. Thus, the denoised image I corresponding to the noisy image J is obtained as

$\begin{array}{cc}I=\underset{Ĩ}{\text{arg}\text{min}}{\displaystyle {\int }_{\Omega }[-{\gamma }_d{\nu }_d(J-Ĩ)+{\delta }_d^{-{\gamma }_d}{e}^{{\gamma }_d(J-Ĩ)}+{\alpha }_d|\nabla Ĩ|]d\Omega }.& (6)\end{array}$

where Ω is the image domain, γ_d, ν_d, δ_d are parameters of the generalized gamma distribution that models the noise and α_d the regularization parameter for denoising.

Then, the optical flow is estimated for the denoised sequence of images using the method proposed in Brox et al. (2004). Particularly, the flow (i.e., the displacement field) is computed between the end-diastolic frame of the sequence to the other frames, corresponding to the different cardiac phases. Thus, the displacement field u^OF = (u^OF, υ^OF) between the end-diastolic frame I^D and the s-phase frame I^s is given by

$\begin{array}{ll}{\text{u}}^{\text{OF}}={\displaystyle \sum _{r=1}^R}\delta {\text{u}}^r,& (7)\end{array}$

where δu^r is the flow component corresponding to the image resolution r that is obtained as

$\begin{array}{l}\delta {\text{u}}^r=\underset{\delta \text{u}}{\text{arg min}}\underset{\Omega }{\underset{}{{\displaystyle \int }}}[\psi \left({\Vert \frac{\partial I^r}{\partial t}+\nabla I^r·\delta \text{u}\Vert }_{G_{\rho }}^2\right)\\ + {\alpha }_o\psi {(\Vert \nabla ({\text{u}}^{r-1}+\delta \text{u})\Vert }_F^2]d\Omega ,& (8)\end{array}$

where ${\text{u}}^{r-1}=\sum _{t=1}^{r-1}\delta {\text{u}}^t$, ∥·∥_F is the Frobenius norm, α_o is the regularization parameter for optical flow. The function ψ and the weighted norm ∥·∥_Gρ are defined by

$\begin{array}{l}\psi (x)=2{\kappa }^2\sqrt{1+\frac{x}{{\kappa }^2}},\\ {\Vert \frac{\partial I^r}{\partial t}+\nabla I^r·\delta \mathbf{\text{u}}\Vert }_{G_{\rho }}^2=G_{\rho }*{\left(\frac{\partial I^r}{\partial x}\right)}^2\delta u^2+G_{\rho }*{\left(\frac{\partial I^r}{\partial y}\right)}^2\delta {\upsilon }^2\\ +G_{\rho }*{\left(\frac{\partial I^r}{\partial t}\right)}^2+2G_{\rho }*\left(\frac{\partial I^r}{\partial x}\frac{\partial I^r}{\partial y}\right)\delta u\delta \upsilon \\ +2G_{\rho }*\left(\frac{\partial I^r}{\partial x}\frac{\partial I^r}{\partial t}\right)\delta u+2G_{\rho }*\left(\frac{\partial I^r}{\partial y}\frac{\partial I^r}{\partial t}\right)\delta \upsilon \end{array}$

where G_ρ is the Gaussian kernel with ρ standard deviation and * is the convolution operator. Note that the flow u^OF is the displacement field between I^D and I^s, then the temporal derivative is estimated as the variation of the intensity between such frames.

Such strategy defines all displacement fields along the cardiac sequence at the same reference phase (the end-diastolic phase), which eases the integrability of the data into the assimilation process introduced in section 2.3.

2.1.4. Patient-Specific Geometric Model

Using an IVUS study gated at the end-diastolic phase, a geometrical model for a frame of interest is constructed (see Figure 2). First, the intima-media area is manually segmented by a specialist from the image using cubic splines to obtain a 2D patient-specific geometry. Then, the 2D geometry is extruded 0.05 mm in the axial direction to render a 3D slice of the arterial vessel. The mesh generation from this geometry is described later in section 2.2.5 when the numerical scheme for the mechanical problem is introduced.

FIGURE 2

Figure 2. Vessel segmentation and geometric model generation: (Left) segmentation of the intima-media area over the IVUS image; (Middle) extruded volume of the intima-media area; (Right) final 3D mesh.

2.2. Mechanical Setup for the Arterial Wall

In this section, the main ingredients from continuum mechanics required to describe the mathematical models are briefly summarized. For further details the reader may refer to Ares (2016) and Blanco et al. (2016).

Let us consider the domain of a cross-sectional slice of the vessel wall. Its spatial configuration in the Euclidean space is denoted by Ω_s, with boundary $\partial {\Omega }_s=\partial {\Omega }_s^W\cup \partial {\Omega }_s^E\cup \partial {\Omega }_s^A$, where $\partial {\Omega }_s^W$ represents the interface between the vessel and the blood, $\partial {\Omega }_s^E$ the external surface, and $\partial {\Omega }_s^A={\bigcup }_{i=1}^2\partial {\Omega }_s^{A,i}$ stands for the set of 2 cross-sectional (non-physical) axial boundaries for the vessel slice (see Figure 3). The unit outward normal vector is denoted by n_s. The coordinates at this configuration are denoted by x_s. A material configuration, used as a reference configuration, is denoted by Ω_m, with coordinates x_m. In the present context, Ω_s stands for the configuration at which mechanical equilibrium is achieved for a given load condition (diastolic, systolic or any other loaded state of the arterial wall). Residual stresses are neglected, therefore, the material configuration Ω_m is both load-free and stress-free.

FIGURE 3

Figure 3. Cross-sectional slice of an arterial vessel. Description of the spatial domain Ω_s and boundaries $\partial {\Omega }_s=\partial {\Omega }_s^W\cup \partial {\Omega }_s^E\cup \partial {\Omega }_s^A$ of the arterial wall.

The displacement field mapping points from the material into the spatial configuration is denoted by u. Then, we characterize the deformation mapping from Ω_m onto Ω_s and its inverse by the following expressions,

$\begin{array}{ll}{\text{x}}_s={\chi }_m({\text{x}}_m)={\text{x}}_m+{\text{u}}_m,& (9)\end{array}$

$\begin{array}{ll}{\text{x}}_m={\chi }_s({\text{x}}_s)={\chi }_m^{-1}({\text{x}}_s)={\text{x}}_s-{\text{u}}_s,& (10)\end{array}$

where subscripts m and s denote the descriptions of the fields in the material and spatial configurations, respectively. Thus, the displacement vector field is given by

$\begin{array}{ll}{\text{u}}_s({\text{x}}_s)={({\text{u}}_m({\text{x}}_m))}_s={\text{u}}_m({\chi }_m^{-1}({\text{x}}_s)),& (11)\end{array}$

and its gradients with respect to material and spatial coordinates are, respectively, obtained as

$\begin{array}{ll}{\text{F}}_m={\nabla }_m{\chi }_m=\text{I}+{\nabla }_m{\text{u}}_m,& (12)\end{array}$

$\begin{array}{ll}{\text{f}}_s={\nabla }_s{\chi }_s={\nabla }_s{\chi }_m^{-1}=\text{I}-{\nabla }_s{\text{u}}_s.& (13)\end{array}$

Observe that ${\left[{\text{F}}_m^{-1}\right]}_s={\text{f}}_s$ and ${\left[{\text{f}}_s^{-1}\right]}_m={\text{F}}_m$. Arterial wall tissues are assumed to behave as incompressible materials, which is mathematically represented by the following kinematic condition

$\begin{array}{ll}det{\text{F}}_m=1.& (14)\end{array}$

In a general case the load state of the model of an arterial cross-section is characterized as follows. Neumann boundary conditions are considered to be given by the forces exerted by the blood flow over $\partial {\Omega }_s^W$, i.e., through a traction field ${\text{t}}_s^W$ which is considered to be characterized as ${\text{t}}_s^W=p_s{\text{n}}_s$ (here we only consider the pressure load, and neglect the shear forces imprinted by the blood flow on the vessel wall), and by the tethering tractions ${\text{t}}_s^{A,i}$ acting over $\partial {\Omega }_s^{A,i}$, i = 1, 2. For ease of notation, the tethering tractions are grouped into ${\text{t}}_s^A$, which is defined over the whole $\partial {\Omega }_s^A$. The action of the surrounding tissues is introduced as an elastic traction over the external boundary, which is characterized by the elastic parameter τ and depends on the displacement field at this boundary (see further details in section 2.2.1). These tractions, representing the external tissues influence, only act over the external surface $\partial {\Omega }_s^E$ in the physiological pressure range, i.e., at end-diastolic pressure or higher. That is, during the preload problem (see section 2.2.2), the boundary ∂Ω_E is an homogeneous Neumann boundary (except for a small region of arc length Δ = 0.1mm which is fixed to remove rigid motions).

The mechanical problem in variational form is framed as a saddle-point problem to accommodate the incompressibility constraint through the corresponding Lagrange multiplier, i.e., the pressure field in the solid domain.

Next, two variational formulations are presented which formalize the concept of mechanical equilibrium for the so-called preload and forward problems. In the preload problem, the known configuration is that one at which the body is at equilibrium (the spatial domain), and the unknown configuration is the material configuration used to define the constitutive equations. In the forward problem, the known configuration is the material one, while the unknown configuration is the one where equilibrium actually occurs.

2.2.1. Forward Problem

When the material (load- and stress-free) configuration Ω_m is known, the variational Equation (16) can be cast in the material domain, yielding what we define as the forward problem. The variational formulation then reads: given the material description of the loads, p_m and ${\text{t}}_m^{A,i}$, find (u_m, λ_m) in ${{U}}_m\times {{L}}_m$ such that

$\begin{array}{l}{\displaystyle \underset{{\Omega }_m}{\int }}(1-det{\text{F}}_m){\widehat{\lambda }}_{m}d{\Omega }_m-{\displaystyle \underset{{\Omega }_m}{\int }}{\lambda }_m({\text{F}}_m^{-T}·{\nabla }_m{\widehat{\text{u}}}_m)det{\text{F}}_{m}d{\Omega }_m\\ +{\displaystyle \underset{{\Omega }_m}{\int }}({\text{S}}_m({\text{E}}_m))·\stackrel{·}{\text{E}}({\widehat{\text{u}}}_m)d{\Omega }_m\hfill \\ ={\displaystyle \underset{\partial {\Omega }_m^E}{\int }}\tau ({\text{u}}_m-{({\text{u}}_d+{\text{u}}^{\text{OF}})}_m)·{\widehat{\text{u}}}_m\left|{\text{F}}_m^{-T}{\text{n}}_0^E\right|det{\text{F}}_{m}d\partial {\Omega }_m^E\hfill \\ +{\displaystyle \underset{\partial {\Omega }_m^W}{\int }}\left(p_m{\text{F}}_m^{-T}{\text{n}}_0^W·{\widehat{\text{u}}}_m\right)det{\text{F}}_{m}d\partial {\Omega }_m^W\hfill \\ +{\displaystyle \sum _{i=1}^2}{\displaystyle \underset{\partial {\Omega }_m^{A,i}}{\int }}({\text{t}}_m^{A,i}·{\widehat{\text{u}}}_m)\left|{\text{F}}_m^{-T}{\text{n}}_0^{A,i}\right|det{\text{F}}_{m}d\partial {\Omega }_m^{A,i}\hfill \\ \forall ({\widehat{\text{u}}}_m,{\widehat{\lambda }}_m)\in {{V}}_m\times {{L}}_m,\hfill & (15)\end{array}$

where $\stackrel{·}{\text{E}}({\widehat{\text{u}}}_m)=\frac{1}{2}\left[{\text{F}}_m^T({\nabla }_m{\widehat{\text{u}}}_m)+{({\nabla }_m{\widehat{\text{u}}}_m)}^T{\text{F}}_m\right]$, n₀ is the unit outward normal vector in the material configuration. Recall that τ is the elastic parameter characterizing the response of the surrounding media, u^OF is the displacement field which maps the end-diastolic to the spatial configuration where equilibrium is achieved (see Equation 7), and u_d is the displacement field which maps points from the material to the end-diastolic configuration. Also, ${{U}}_m$, ${{V}}_m$, and ${{L}}_m$ are the counterparts of ${{U}}_s$, ${{V}}_s$, and ${{L}}_s$, respectively, with functions defined in Ω_m.

Acceleration terms have also been neglected, since the stresses associated to such inertial forces is much smaller than those of constitutive origin (Ares, 2016; Blanco et al., 2016).

2.2.2. Preload Problem

Given the equilibrium configuration Ω_s, the variational formulation reads: given the loads ${\text{t}}_s^{W,n}$ and ${\text{t}}_s^A$, find $({\text{u}}_s,{\lambda }_s)\in {{U}}_s\times {{L}}_s$ such that

$\begin{array}{l}{\displaystyle \underset{{\Omega }_s}{\int }}\left[-{\lambda }_{s}div{\widehat{\text{u}}}_s+{\sigma }_s·{\epsilon }_s({\widehat{\text{u}}}_s)\right]d{\Omega }_s-{\displaystyle \underset{{\Omega }_s}{\int }}\left[1-det{\text{F}}_s^{-1}\right]{\widehat{\lambda }}_{s}d{\Omega }_s=\\ {\displaystyle \underset{\partial {\Omega }_s^W}{\int }}{t}_s^{W,n}{\text{n}}_s·{\widehat{\text{u}}}_{s}d\partial {\Omega }_s^W+{\displaystyle \sum _{i=1}^2}{\displaystyle \underset{\partial {\Omega }_s^{A,i}}{\int }}{\text{t}}_s^{A,i}·{\widehat{\text{u}}}_{s}d\partial {\Omega }_s^{A,i}\\ \forall ({\widehat{\text{u}}}_s,{\widehat{\lambda }}_s)\in {{V}}_s\times {{L}}_s,\hfill & (16)\end{array}$

where ${\epsilon }_s(\widehat{\text{u}})=\frac{1}{2}({\nabla }_s\widehat{\text{u}}+{\nabla }_s{\widehat{\text{u}}}^T)$ is the strain rate tensor, ${{L}}_s=L^2({\Omega }_s)$ and ${{U}}_s=\left\{{\text{u}}_s\in {\text{H}}^1({\Omega }_s),{\text{u}}_s\mathrm{\text{satisfies essential b.c.}}\right\}$ are, respectively, the linear space for pressures and the linear manifold for kinematically admissible displacements, and ${{V}}_s={\widehat{\text{u}}}_s\in {\text{H}}^1({\Omega }_s),{\widehat{\text{u}}}_s$ satisfies homogeneous essential b.c.} is the space of kinematically admissible variations. Also, σ_s is related to the second Piola-Kirchhoff stress tensor S_m through

$\begin{array}{ll}{\sigma }_s=\frac{1}{det{\text{F}}_s}{\text{F}}_s{({\text{S}}_m({\text{E}}_m))}_s{\text{F}}_s^T.& (17)\end{array}$

where S_m is a function of the Green-Lagrange deformation tensor ${\text{E}}_m=\frac{1}{2}({\text{F}}_m^T{\text{F}}_m-\text{I})$ via a constitutive equation (see section 2.2.4).

In this work the preload problem is used to obtain the material configuration that enables an appropriate calculation of the stress field, which realizes the equilibrium in the end-diastolic configuration. Note that in this case the action of the surrounding media is omitted. This is due to the fact that our hypothesis considers the end-diastolic configuration as a reference configuration for the elastic response of the external tissues.

2.2.3. Equilibrium Problems for a Given Set of Material Parameters

The preload problem is a mandatory step toward characterizing the mechanical state (the stress state) of the arterial wall in a geometry obtained from medical images (e.g., the end-diastolic geometry) with given baseline hemodynamics loads. In this context, these loads are given by the end-diastolic pressure, also called preload pressure, and by the axial stretch caused by tethering forces. The material configuration is required because it is used to define constitutive equations, without which the forward problem cannot properly be formulated. In our case, such baseline geometry is obtained from IVUS study, while the baseline hemodynamics loads (the blood pressure) are estimated from patient specific data. Just after solving the preload problem, the baseline mechanical state, that is the stress state due to the preload pressure (i.e., pressure at diastole), is adequately determined and the displacement field u_s—that maps the material (load-free) configuration to the diastolic configuration– is recalled as u_d. Then, the forward problem is solved to determine the equilibrium configuration for other hemodynamics loads occurring during the cardiac cycle. In that manner both problems are synergically coupled to solve a forward problem from an adequately preloaded configuration.

In practice, a set of physiological loads for the vessel will be given. Individualizing the diastolic pressure level as p_s, a set of pressure loads between diastole and systole can be listed as {p_s1, …, p_sS}. Through the forward problem, each load p_si will be in correspondence with an unknown spatial configuration Ω_si, i = 1, …, S. Notice then that, for a given set of material parameters, the preload problem is solved only once and so the forward problem is solved for each load p_si in the set of physiological loads.

2.2.4. Constitutive Models

The main components of the atherosclerotic plaque, i.e., fibrotic, lipidic and calcified tissues, are modeled as isotropic Neo-Hookean materials. In Walsh et al. (2014), it is shown that fibrotic tissue in illiac plaque presents a quasi isotropic behavior. Different from the fibrotic tissue, the lipidic and calcified tissues do not display any contribution of smooth muscle cells or oriented fibers that may endow their structures with anisotropic behavior, what suggests that an isotropic hypothesis for these materials is reasonable.

The isotropic Neo-Hookean model is suitable for materials under large deformations where the stress-strain relationship behaves as non-linear, elastic, isotropic and independent of strain rate. Also, the model assumes an ideal elastic material at every strain level which, for physiological ranges, is satisfied by many biological tissues. The stress-strain relationship for a Neo-Hookean material derives from the strain energy function

$\begin{array}{ll}\psi =\frac{c}{2}({\overline{I}}_1-3),& (18)\end{array}$

where c is the material parameter that characterizes the stiffness of the material and ${\overline{I}}_1$ is the first isochoric invariant of the Cauchy-Green tensor

$\begin{array}{ll}{\overline{I}}_1=Tr{(\text{C}}_m{(det{\text{F}}_m)}^{-2/3}),& (19)\end{array}$

with ${\text{C}}_m={\text{F}}_m^T{\text{F}}_m$. Then, the second Piola-Kirchhoff stress tensor (and the σ_m through Equation 17) is obtained as

$\begin{array}{ll}{\text{S}}_m({\text{E}}_m)=\frac{\partial \psi }{\partial {\text{E}}_m}.& (20)\end{array}$

2.2.5. Numerical Methods

The preload and forward problems are linearized using the Newton-Raphson method. Linear tetrahedral finite elements for both displacement and pressure fields are used for the spatial discretization of the corresponding linearized problems. To stabilize the problem in the sense of the inf-sup condition, the linearized (forward and preload) problems are modified adding a diffusive term in the pressure equation. For the analysis of the proposed approach, four patient-specific 3D geometries were obtained using the technique described in section 2.1.4. These geometries were discretized using Netgen 3D using a characteristic element size ranging from 10μm to 40μm, resulting in meshes with 6,521, 7,516, 4,835, and 3,808 nodes for the cases 1–4, respectively. All these steps are performed using an in-house solver. The resulting systems of linear equations are solved using a direct solver based on LU factorization from the SuperLU library (Li and Demmel, 2003). Further details regarding the linearization and numerical schemes can be found in Ares (2016) and Blanco et al. (2016).

The Newton iterative scheme in both equilibrium problems finishes when $\parallel {\text{u}}_s^{m+1}-{\text{u}}_s^m{\parallel }_{L^{\infty }}<10^{-4}$ mm and $\parallel {\lambda }_s^{n+1}-{\lambda }_s^n{\parallel }_{L^{\infty }}<1$ Pa. Such convergence criterion was chosen to yield a higher precision than the optical flow processing applied to IVUS images (16·10⁻³ mm assuming pixel precision).

2.3. Data Assimilation

In the data assimilation process, the displacement field u^OF obtained using the optical flow technique as explained in section 2.1 and the mechanical models presented in the previous section (section 2.2) are integrated by an unscented Kalman filter. Let us define a partition for the domain of analysis ${\Omega }_s={\bigcup }_{j=1}^M{\Omega }_s^j$ composed by M disjoint regions. Each region ${\Omega }_s^i$ is characterized by its own material parameter, say c_i, see Equation (18). The axial loads ${\text{t}}_{s_i}^{A,n}$, the pressure level p_si and the displacement fields ${\text{u}}_{s_i}^{\text{OF}}$ (obtained by optical flow techniques) are known at S cardiac phases (i = 1, …, S). Since our mechanical problem is time-independent, the time instants in the context of the Kalman filter simply correspond to filter iterations, while at each iteration all forward problems must be solved. By using the mechanical constitutive models, the material parameters grouped as θ = (c₁, …, c_M), are estimated such that

$\begin{array}{ll}\theta ={\displaystyle \underset{\widehat{\theta }}{\text{arg min}}}{\displaystyle \sum _{i=1}^S}\parallel {\text{u}}_{s_i}^{\text{MO}}(\widehat{\theta })-{\text{u}}_{s_i}^{\text{OF}}{\parallel }_{L_2}^2,& (21)\end{array}$

where ${\text{u}}_{s_i}^{\text{MO}}(\widehat{\theta })$ is the displacement field at the configuration s_i obtained by solving the preload and forward problems (described in section 2.2) with pressure level p_si and material parameters $\widehat{\theta }$.

The solution of the parameter identification problem eqution (21), satisfies the discrete dynamic nonlinear system presented as follows

$\begin{array}{l}\begin{array}{c}{X}_k^a=f(X_{k-1}^a,t_{k-1})+W_k,\\ Z_k=h(X_k^a,t_k)+V_k,\hfill & (22)\end{array}\end{array}$

where $X_k^a$ is the augmented state vector

$\begin{array}{ll}{X}_k^a={\left[{\text{u}}_{s_1}^k(\text{x}),\dots ,{\text{u}}_{s_S}^k(\text{x}),{\lambda }_{s_1}^k(\text{x}),\dots ,{\lambda }_{s_S}^k(\text{x}),c_1,\dots ,c_M\right]}^T,& (23)\end{array}$

which contains the displacement u_si and pressure λ_si fields for all forward problems i = 1, …, S, and the material parameters of all regions of the domain θ = (c₁, …, c_M); $f(X_k^a,t_k)$ is the operator that sequentially solves the preload and all forward problems for parameters and initial state conditions in $X_k^a$ at filter iteration t_k (recall that these problems are time-independent, and so the dependence on time is ruled out in practice); W_k are the model errors at the k-th step; $h(X_k^a,t_k)=\text{H}{X}_k^a$ is a linear observation operator represented by the block matrix

$\begin{array}{ll}\text{H}=\left[\begin{array}{ccc}{I}_{uu}& 0_{u\lambda }& 0_{u\theta }\end{array}\right],& (24)\end{array}$

where block matrix indexes indicate the corresponding dimensions; Z is the set of optical flow observations at each cardiac phase, described by the column vector

$\begin{array}{ll}Z={\left[{\text{u}}_{s_1}^{\text{OF}}(\text{x}),\dots ,{\text{u}}_{s_S}^{\text{OF}}(\text{x})\right]}^T,& (25)\end{array}$

where ${\text{u}}_{s_i}^{\text{OF}}(\text{x})$ is the displacement field obtained by the optical flow technique for the cardiac phase i, i = 1, …, S (observe that for the present case of static problems, the observations are fixed concerning the dynamics of the data assimilation process); V is the vector of optical flow and interpolation errors for the observation vector Z.

To obtain an estimate of the parameters θ, a reduced ordered unscented Kalman filter (ROUKF) (Julier and Uhlmann, 2002, 2004) is applied to the system described in Equation (22). The filter comprises the following steps

1. Spherical sigma-points generation ${\sigma }_i^{(n)},i=1,\dots ,M+1$ with their corresponding weights w⁽ⁱ⁾ (see Julier, 2003) and initialization of the variables

$\begin{array}{l}{R}_0={\sigma }_{\text{OF}}{I}_{uu};L_0=\left[\begin{array}{c}{L}_0^X\\ L_0^{\theta }\end{array}\right]=\left[\begin{array}{c}{L}_0^u\\ L_0^{\lambda }\\ L_0^{\theta }\end{array}\right]=\left[\begin{array}{c}{0}_{u\theta }\\ 0_{\lambda \theta }\\ I_{\theta \theta }\end{array}\right];\\ U_0^{-1}=\left[\begin{array}{ccc}{\sigma }_{{\widehat{c}}_1}& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & {\sigma }_{{\widehat{c}}_M}\end{array}\right],& (26)\end{array}$

$\begin{array}{ll}{X}_0^a={\left[{\widehat{X}}_0^{+},{\widehat{\theta }}_0^{+}\right]}^T={\left[0_{\text{u}},0_{\lambda },{\widehat{\theta }}_0\right]}^T,& (27)\end{array}$

$\begin{array}{ll}{\text{P}}_0^{+}={\text{L}}_0{\text{U}}_0^{-1}{\text{L}}_0^T,& (28)\end{array}$

where σ_OF is the uncertainty of the computed optical flow and σ_ĉi is the uncertainty of the parameter c_i, i = 1, …, M. The sensitivity analysis of the uncertainty value is studied in section 3.1.

2. The prediction step

$\begin{array}{l}{\widehat{X}}_{k-1}^{(i)}={\widehat{X}}_{k-1}^{+}+L_{k-1}^X\sqrt{U_{k-1}^{-1}}{\sigma }_i^{(n)},i=1,\dots ,M+1,\\ {\widehat{\theta }}_{k-1}^{(i)}={\widehat{\theta }}_{k-1}^{+}+L_{k-1}^{\theta }\sqrt{U_{k-1}^{-1}}{\sigma }_i^{(n)},i=1,\dots ,M+1,\\ \left[\begin{array}{c}({\widehat{X}}_k^{(i)})\\ ({\widehat{\theta }}_k^{(i)})\end{array}\right]=f\left(\left[\begin{array}{c}({\widehat{X}}_{k-1}^{(i)})\\ ({\widehat{\theta }}_{k-1}^{(i)})\end{array}\right],t_{k-1}\right),\\ {\widehat{X}}_k^{-}={\displaystyle \sum _{i=1}^{M+1}{w}^{(i)}}{\widehat{X}}_k^{(i)},{\widehat{\theta }}_k^{-}={\displaystyle \sum _{i=1}^{M+1}{w}^{(i)}}{\widehat{\theta }}_k^{(i)},{\widehat{Z}}_k={\displaystyle \sum _{i=1}^{M+1}{w}^{(i)}}{\widehat{Z}}_k^{(i)}.& (29)\end{array}$

3. The correction step

$\begin{array}{l}{L}_k^X={\widehat{X}}_k^{(*)}{D}_w{\left({\sigma }^{(*)}\right)}^T,L_k^{\theta }={\widehat{\theta }}_k^{(*)}{D}_w{\left({\sigma }^{(*)}\right)}^T,\\ {\left\{HL\right\}}_k={\widehat{Z}}_k^{(*)}{D}_w{\left({\sigma }^{(*)}\right)}^T,\\ P_w={\sigma }^{(*)}{D}_w{\left({\sigma }^{(*)}\right)}^T,\\ U_k=P_w+{\left\{HL\right\}}_k^T{R}_k^{-1}{\left\{HL\right\}}_k,\\ {\widehat{X}}_k^{+}={\widehat{X}}_k^{-}+L_k^X{U}_k^{-1}{\left\{HL\right\}}_k^T{R}_k^{-1}\left(Z-{\widehat{Z}}_k\right),\\ {\widehat{\theta }}_k^{+}={\widehat{\theta }}_k^{-}+L_k^{\theta }{U}_k^{-1}{\left\{HL\right\}}_k^T{R}_k^{-1}\left(Z-{\widehat{Z}}_k\right).& (30)\end{array}$

The matrices ${\sigma }^{(*)},{\widehat{\text{X}}}_k^{(*)},{\widehat{\text{Z}}}_k^{(*)},{\widehat{\theta }}_k^{(*)}$ are the M × (M+1) matrices whose columns are the vectors ${\sigma }^{(i)},{\widehat{X}}_k^{(i)},{Ẑ}_k^{(i)},{\widehat{\theta }}_k^{(i)}$ with i = 1, …, M+1, respectively. D_w is the diagonal (M+1) × (M+1) matrix with values $D_{ii}=w^{(i)},i=1,\dots ,M+1$, i.e., the sigma-point weights.

4. If stop criteria is not achieved, go to step 2 and k = k+1.

In this iterative scheme, the model errors W_k (inaccuracies in the solution of the preload and forward problems) have been neglected. The stop criteria used in this work is a fixed number of iterations that is reported for each study case in section 3.

In this work, c was reparametrized as $c=2^{\widehat{\theta }}$ (this approach was introduced in Bertoglio et al., 2012) allowing $\widehat{\theta }$ to vary in the whole ℝ (as occurs in the presented formulation 29, 30) without delivering invalid values for c.¹

2.4. Parallelization Scheme

The data assimilation scheme is a computationally demanding task. However, it presents many independent or low dependent tasks. Firstly, notice that all sigma point predictions can be computed in parallel. As the forward problem is static, all forward problems (one per load, for S different loads) are computed in parallel and an extended observation vector Ẑ⁽ⁱ⁾ = [Ẑ^{(i), 1}, …, Ẑ^{(i), S}]T for the i-th sigma-point is created by appending the predicted displacements Ẑ^{(i), j} of the j-th forward problem corresponding to the pressure load p_sj. In that manner, at each Kalman iteration, the observations of all frames, are processed at once. In turn, the forward problem itself is parallelized by partitioning the mesh and communicating among subdomains the results of the local operations in both assembling and solving stages. Partitioning is accomplished using ParMETIS (Karypis and Kumar, 1998), and the solution is achieved using the SuperLU library (Li and Demmel, 2003). Following such parallelization scheme, and assuming there are enough computational resources, the cost per iteration of the data assimilation process equals the cost of the computation of one preload problem plus one forward problem, regardless the number of cardiac phases or sigma-points employed (i.e., regardless the number of parameters to be estimated). Note that the cost of the Kalman filter increases as more parameters are estimated, although when compared to the computations required for solving the mechanical equilibrium problems this increment is insignificant (only a few dozens of parameters will be required in the worst case). In Figure 4, the activity diagram for the proposed parallel scheme is presented. Thus, the data assimilation process is HPC ready and, even, capable to handle large scale FEM problems.

FIGURE 4

Figure 4. Activity diagram depicting the parallel workflow for data assimilation. The three levels of parallelism are highlighted on the right side: (i) Parallelization of the sigma points that are solved at the same time without communication among the threads; (ii) Parallelization of each load state of the artery (i.e., one load per cardiac phase) that is fully parallelized without communication among the threads; and (iii) Parallelization of the FEM problem by mesh partitioning.

3. Results

In what follows, sensitivity analyses are carried out to study the variation of the parameter estimation with respect to: the parameter uncertainties, boundary conditions and baseline stress state of the mechanical model (sections 3.1, 3.2, and 3.3, respectively). From these analyses, a reasonable setup of the data assimilation parameters and mechanical conditions is obtained for the present context of material identification in patient-specific models. Finally, in section 3.4, 4 patient-specific mechanical models are derived from in-vivo IVUS studies and the obtained displacement errors between the model predictions—with its parameters adjusted by data assimilation—and the optical flow observations are assessed.

3.1. Uncertainty Parameters Sensitivity

Let us define a homogeneous ring-shaped domain Ω_s with Neo-Hookean constitutive behavior (see Equation 18). The inner and outer radius of the ring are 2 and 2.71mm, respectively. The size and proportions are chosen to approximate an idealized coronary artery. Loads of t^{W, n} = 80 mmHg and t^{W, n} = 120 mmHg are applied over the inner surface for the preload and forward problems, respectively, and tethering tractions ${\text{t}}_s^{A,i}$ are considered such that an axial stretch of 10% is prescribed. At the outer surface, homogeneous Neumann boundary conditions are assumed (τ = 0 in Equation 16). To avoid rigid movements in this idealized geometry, only radial displacement is allowed for 4 equidistant nodes at the luminar perimeter. The forward operator f, which comprises the preload and forward problems (see Equations 15, 16), is solved at each filter iteration with an iterative scheme where a Newton-Raphson linearization procedure is applied as described in section 2.2.5 (further details in Blanco et al., 2016).

Using this setting, we create an in-silico experiment to analyze: (i) the sensitivity of the parameter estimates $\widehat{\theta }$ with respect to the σ_Z (the observations uncertainty, previously referred to as σ_OF); and (ii) the sensitivity of the parameter estimates $\widehat{\theta }$ with respect to the σ_θ (the estimate uncertainty). Thus, the observations are generated by computing Z = h(f(X^t)) where $X^t=\left[0_{\text{u}},0_{\lambda },{\theta }^t\right]$ is the true augmented state vector with the solution parameters c^t = 2^θt for the experiment. In this particular case, the domain is homogeneous and the constitutive model has only one parameter (c), then, only one parameter is estimated.

To analyze the sensitivity of $\widehat{\theta }$ with respect to the observation uncertainty σ_Z, the estimation of the parameter is performed assuming different values σ_Z, ranging from 10⁻¹ to 10⁻⁵ mm. Also, three different materials are used for the ring, mimicking: cellular fibrotic tissue (c^t = 5·10⁵Pa), lipidic tissue (c^t = 1·10⁵Pa) and calcified tissue (c^t = 4·10⁶Pa). The estimation of the Kalman filter for all the 15 cases is presented in Figure 5. The results showed that in all cases the parameter uncertainty interval $[2^{\widehat{\theta }-\sqrt{\text{diag}\left(U^{-1}\right)}};2^{\widehat{\theta }+\sqrt{\text{diag}\left((U^{-1})\right)}}]$ encloses the true parameter value c^t. Even though, a closer estimate across the three materials is obtained for ${\sigma }_Z=10^{-3}$ mm which seems reasonable as it is the precision of the displacements delivered by the convergence process in solving the nonlinear operator f.

FIGURE 5

Figure 5. Sensitivity analysis of the parameter estimation with respect to σ_Z for the experiment of the 1-material ring fixing the parameter uncertainty σ_θ = 2: (Top) Estimate using observation uncertainties of ${\sigma }_Z=10^i$ mm, i = −5, …, −1 where each dot corresponds to a data assimilation process (after 200 iterations). The color indicates the material of the ring at each experiment, the dashed line is c^t value and the whiskers denote the parameter uncertainty interval. (Bottom) Convergence of the Kalman filter using observation uncertainties of ${\sigma }_Z=10^i$ mm, i = −3, −2, −1 for c^t = 4·10⁻⁶ Pa. The dashed line is c^t value, the solid line the Kalman filter estimate $2^{\widehat{\theta }}$ and the colored ribbon denotes the parameter uncertainty interval. In both cases, the uncertainty interval is estimated as $[2^{\widehat{\theta }-\sqrt{\text{diag}(U^{-1})}};2^{\widehat{\theta }+\sqrt{\text{diag}(U^{-1})}}]$.

Regarding the filter convergence, it is observed that as the uncertainty in the observations decreases, the method converges faster. In Figure 5, it is shown that as the σ_Z increases its value, the convergence is slower. Note that the estimator gain matrix is computed as ${\text{K}}_k={\text{L}}_k^{\theta }{\text{U}}_k^{-1}{\left\{\text{H}\text{L}\right\}}_k^T{\text{R}}_k^{-1}$ and the only operator that varies in the first iteration of the presented cases is ${\text{R}}_0^{-1}$. As the spectral radius of ${\text{R}}_0^{-1}$ diminishes as σ_Z increases then K₀ spectral radius diminishes as well, yielding a smaller correction of ${\widehat{\theta }}_k^{+}$ as presented in the plot. At the same time, since P_w is constant, the update of ${\text{U}}_k={\text{P}}_w+{\left\{\text{H}\text{L}\right\}}_k^T{\text{R}}_k^{-1}{\left\{\text{H}\text{L}\right\}}_k$ is damped by R_k. This damping effect is evidenced in the evolution of the parameter uncertainty intervals plotted in Figure 5. In statistical terms, the lack of confidence in the new observations leads to reducing its weight at the correction step.

An analogous analysis was performed to study the sensitivity of $\widehat{\theta }$ with respect to the parameter uncertainty σ_θ. The uncertainty levels for σ_θ ranged from 0.25 to 4 and the experiment was repeated for the three different ring materials (fibrotic, lipidic and calcified tissues). The results showed that the bigger σ_θ, the wider the search space for the parameter, and the faster the method converges when the initial value is far from the true parameter value (see Figure 6). On the other hand, high values of σ_θ may cause an overshooting in the estimation and a slower convergence. In this scenario, the reparametrization deteriorates the convergence even more. The reparametrization imposes an estimation bias to stiffer values due to the fact that displacements are less sensitive with respect to small variations in stiffer than softer materials. Then, the mean observation error (used as correction term in Equation 30) is biased to the sigma points associated with stiffer materials. This is clearly evidenced in Figure 7, where the initial overshooting delays the estimation of the parameter. Hence, note that the initial uncertainty interval does not necessarily has to contain c^t to estimate its correct value. In fact, the uncertainty parameter values are also iteratively updated and similar values are obtained for all three σ_θ illustrated in Figure 6. The role of the initial value of σ_θ is the dispersion of sigma points around the mean initial guess, and large values may accelerate convergence when the initial guess c_i is far from c^t.

FIGURE 6

Figure 6. Sensitivity analysis of the parameter estimation with respect to σ_θ for the experiment of the 1-material ring fixing the parameter uncertainty ${\sigma }_Z=10^{-2}$: (Top) Estimate using parameter uncertainties of σ_θ = 0.25, 0.5, 1, 2, 4 where each dot corresponds to a data assimilation process (after 200 iterations). The color indicates the material of the ring at each experiment, the dashed line is c^t value and the whiskers denote the parameter uncertainty interval. (Bottom) Convergence of the Kalman filter using parameter uncertainties of σ_θ = 0.25, 1, 4 and c^t = 4·10⁶ Pa. The dashed line is c^t value, the solid line the Kalman filter estimate $2^{\widehat{\theta }}$ and the colored ribbon denotes the parameter uncertainty interval. In both cases, the uncertainty interval is estimated as $[2^{\widehat{\theta }-\sqrt{\text{diag}(U^{-1})}};2^{\widehat{\theta }+\sqrt{\text{diag}(U^{-1})}}]$.

FIGURE 7

Figure 7. Convergence of the Kalman filter for the experiment of the 1-material ring using σ_θ = 4, 16, ${\sigma }_Z=10^{-2}$ mm and 2^θt = 1·10⁵ Pa. The dashed line is 2^θt value, the solid line the Kalman filter estimate $\widehat{\theta }$ and the colored ribbon denotes the parameter uncertainty interval $[2^{\widehat{\theta }-\sqrt{\text{diag}(U^{-1})}};2^{\widehat{\theta }+\sqrt{\text{diag}(U^{-1})}}]$.

Overall, a good agreement is found in term of accuracy and convergence for parameters ${\sigma }_Z=10^{-3}$ mm and σ_θ = 4. These parameters identify clearly the three different kinds of tissues in this idealized problem. Also, the observations generated in-silico present an accuracy of similar order than the obtained (assuming no error carried by the optical flow) through the IVUS image processing. For this reason, ${\sigma }_Z=10^{-3}$ mm is used in cases analyzed in forthcoming sections. The value of σ_θ cannot be straightforwardly assigned because parameter overshooting using in-vivo complex geometries may lead to excessively soft materials which could cause contact at the inner surface in when solving the preload equilibrium, yielding non-free material configurations. Since stress-free configurations have been assumed, a more conservative value of σ_θ = 0.5 is used to avoid such problem.

3.2. Boundary Conditions Sensitivity

As described in section 2.2.1, the observational datum u^OF is considered as an additional information over ∂Ω^E through a penalization factor τ (i.e., a Robin boundary condition). This strategy is an attempt to incorporate the contribution of surrounding tissues through a surrogate surface model. Moreover, since u^OF can be exposed to errors caused by brightness variations, image artifacts or non-physical optical flow regularization issues, the use of a Robin boundary condition allows the model to naturally filter out the field u^OF similarly as a surface spring model. Then, a characterization of the surrounding tissues provided by τ in the parameter estimation is addressed in this section.

The in-silico study case used for this sensitivity analysis was generated from the cross-section IVUS image depicted in Figure 8 by considering the configurations corresponding to two cardiac phases: end-diastole and systole. The geometrical model was constructed for the end-diastolic configuration following the pipeline described in sections 2.1 and 2.2.5. The configurations at each one of the cardiac phases are related to an end-diastolic load (i.e., the preload) of t^{W, n} = 80mmHg and to a systolic load t^{W, n} = 120mmHg, accordingly. The loads are applied over the inner surface of the vessel in the preload and forward problems, respectively. Finally, tethering tractions ${\text{t}}_s^{A,i}$ are considered such that an axial stretch of 10% is prescribed in the end-diastolic configuration. The remaining setup of boundary conditions is defined for each of the following analyses: (i) parameter estimation sensitivity as τ decreases from a large value (almost Dirichlet condition) to a small value (almost Neumann condition); (ii) parameter estimation robustness when observation u^OF features errors at the boundaries.

FIGURE 8

Figure 8. Optical flow and proposed ground truth for in-silico test: (Top) IVUS image at end-diastolic phase; (Bottom-left) optical flow u^OF used in the data assimilation process extracted from the in-vivo IVUS pair of images between end-diastole and systole phases; (Bottom-right) proposed ground truth ${\text{u}}_2^{\text{GT}}$ generated using homogeneous Neumann boundary condition over the solid red line. The green dashed line indicates the position of the guidewire artifact in the image, and therefore the area in which the optical flow displacement field can be largely affected.

3.2.1. Test 1: Sensitivity of τ for Error-Free Observations

For this analysis, the observations for the ROUKF were generated by solving the mechanical equilibrium with our model, avoiding observational and modeling errors. Thus, a Robin boundary condition was imposed at the outer surface in the forward problem with τ = 10⁶ (practically yielding a Dirichlet boundary condition). This setting rendered a ground truth displacement field ${\text{u}}_1^{\text{GT}}$ for this test. Using observations ${\text{u}}_1^{\text{GT}}$, the data assimlation algorithm was executed for τ ∈ {10⁶, 10⁴, 10²} (higher values of τ were not analyzed since τ = 10⁶ is already almost a Dirichlet boundary condition). The geometric model was partitioned in sextants with two concentric layer yielding 12 regions each with its own material parameter c_i.

The results are presented in Figure 9, depicting the parameter estimation and predicted observations variations as the Robin boundary condition moves toward a Neumann boundary condition. The decrease of forces at the boundaries caused by the decreasing value of τ is compensated by the estimation of softer materials (which experiment higher strains) to match the ${\text{u}}_1^{\text{GT}}$ observations. Particularly, the method recovers the correct material parameters when the penalization value is the true value used to generate the observations. i.e., τ = 10⁶. For the parameter estimation with τ = 10⁴, a qualitatively similar distribution of materials is observed with an uniform reduction in the magnitude of the material parameter. The lowest penalization value, τ = 10², delivers a totally different arrangement of materials. This result emphasizes the important contribution of the surrounding tissues for a correct estimation of material parameters, which is clearly retrieved when sufficient large values of the penalization parameter τ are employed.

FIGURE 9

Figure 9. Sensitivity of the parameter estimation with respect to the penalization factor τ. (Top) In-silico model with the ground truth c^t parameters and boundary condition model with τ = 10⁶ for the generation of $Z=h(f(\left[0_{\text{u}},0_{\lambda },{\theta }^t\right]))$. (Middle) Estimated parameters $\widehat{\theta }$ from the observation Z by using τ = 10⁶, 10⁴, 10² (from left to right) in the forward operator f during the data assimilation process. (Bottom) Displacement error (ε_Z = Z−Ẑ_k mm) for the parameter estimation process using τ = 10⁶, 10⁴, 10² (from left to right) in the forward operator f.

The observation error |ε_Z|, which is defined as the Euclidean distance between the observations Z and mean filter observation Ẑ_k at the last iteration, increases as τ is decreased. Specifically, the mean values of ε_Z are 7.16·10⁻⁵, 1.03·10⁻³ and 9.79·10⁻³ for τ = 10⁶, 10⁴ and 10², respectively. Clearly, should the observations ${\text{u}}_1^{\text{GT}}$ be error-free at the boundaries, a Dirichlet boundary condition (a higher value for τ) would be the correct choice. Notwithstanding this, the observations from in-vivo scenarios are degraded by diverse sources of errors and, as it will be shown next, an excessively stringent boundary condition (of Dirichlet type) may not be the best option.

3.2.2. Test 2: Sensitivity of τ for Realistic Observations

The current analysis aims at assessing the robustness of the parameter estimation process when the u^OF at the outer boundary differs from the real in-vivo displacements. For this purpose, an allegedly ground truth ${\text{u}}_2^{\text{GT}}$ is generated by altering the observation u^OF in a certain region (untrusted region) using the mechanical model with Neumann conditions. Finally, the assimilation process is performed with the observation u^OF and different values of τ, to assess if it is capable to approximate ${\text{u}}_2^{\text{GT}}$ despite the observation errors.

Thus, an IVUS sequence with a swinging artifact (induced by the guidewire) was chosen to perform our analysis. The IVUS cross-section depicted in Figure 8 presents an image artifact from the IVUS guidewire at the bottom-right quadrant of the frame. The guidewire projects a shadow that hides the arterial wall and, as consequence, the optical flow is polluted with a swinging movement not related with the true arterial-wall motion. Thus, a displacement field, denoted by ${\text{u}}_2^{\text{GT}}$, is generated from the in-vivo data removing the guidewire influence, with the purpose of comparing this ground truth against the Kalman predictions Ẑ_k when the polluted optical flow u^OF is used as observations. In that manner, the difference ${\epsilon }^{\text{GT}}={Ẑ}_k-{\text{u}}_2^{\text{GT}}$ can be regarded as an estimate of the error in the Kalman prediction due to the artifact in the image processing data. At last, ε^GT is computed for different values of τ to assess the discrepancies in the predictions as the external Robin boundary condition is characterized differently.

The displacement ${\text{u}}_2^{\text{GT}}$ is generated by solving the equilibrium problems with a model constituted by a single material. To define a reasonable value for this constitutive property, a data assimilation process was performed using u^OF as observation and τ = 10⁴, yielding to c = 33.52 kPa. Note that the c is biased by the image artifact among other errors in the displacement field and it cannot be regarded as an estimate of the real material, thus, it is analyzed the ranges among which the estimated ĉ varies. At the boundary $\partial {\Omega }_m^E$, a Neumann homogeneous condition (τ = 0) was applied in the area affected by the guidewire (see red line in Figure 8) and a Robin boundary condition with τ = 10⁴ was applied to the remaining part of the boundary. The obtained displacement field ${\text{u}}_2^{\text{GT}}$ is displayed in Figure 8.

The sensitivity of ε^GT with respect to τ is then studied. For each value of τ ∈ {10ⁱ, i = 5, 4, …, 0}, the data assimilation process is executed using u^OF as observation. The relative difference between the generated ground truth ${\text{u}}_2^{\text{GT}}$ and the Kalman predicted observation Z_k, for each τ, is reported in Figure 10. For τ greater than 10³, the Robin condition guarantees that the artifact-related displacements are preserved regardless the impact on the induced internal stresses. When τ varies from 10³ to 10², the relative error difference significantly drops at the guidewire locus, from 1.33 to 0.58. As τ decreases even more, the resulting force induced by the Robin boundary condition diminishes its magnitude, yielding lower internal stresses, and spreading the error outside the region of the guidewire shadow. For values lower than 10², the error in the displacement field is concentrated at the bottom area of the artery. Particularly, this concentration of the error is explained by the fact the continuum model is enforced to behave as incompressible, while the optical flow is not divergence-free. In terms of the parameter estimation, the value of c was of 123.96, 33.52, 27.93, 68.51, 144.88, and 125.04 kPa for τ = 10⁵, 10⁴, 10³, 10², 10¹, and 10⁰ respectively, presenting mean and standard deviation value of 87.31±50.70 kPa, all close to a cellular fibrotic tissue. Moreover, there is a large sensitivity in the estimated parameter with respect to the chosen value of τ. In comparison with the ground truth, the closest matching prediction in terms of the displacement field (i.e., the prediction for τ = 10², see Figure 10) presents an estimation of c two times higher. This is a clear demonstration of the large sensitivity in the estimated parameter with respect to the setting of models for the external tissues. Even more, it indicates that the minimization of the displacement field is not directly related to the best parameter estimation.

FIGURE 10

Figure 10. Relative difference ${\epsilon }_r^{\text{GT}}=\frac{{\epsilon }^{\text{GT}}}{\langle \parallel {\text{u}}_2^{\text{GT}}\parallel \rangle }$ with respect to the boundary parameter τ. 〈·〉 Denotes the mean value in Ω_s. From top-to-bottom and left-to-right, ${\epsilon }_r^{\text{GT}}$ is presented for τ = 10⁵, 10⁴, …, 10⁰.The estimation of c was of 123.96, 33.52, 27.93, 68.51, 144.88 and 125.04 kPa for τ = 10⁵, 10⁴, 10³, 10², 10¹, and 10⁰ respectively, while the value of c in the generated ground truth is c = 33.52 kPa.

3.3. Effects of Preload and Axial Stretch

An appropriate baseline stress state of the vessel is key toward an accurate characterization of the stress state in arterial tissues. In fact, as reported in Ares (2016), a preloaded and axially stretched artery features notoriously different stress patterns compared to the case when such loads are neglected. Therefore, it is important to quantify the change in the parameter estimation when the initial stress state is either considered or not in the analysis. To quantify such disagreement, the parameters of an in-vivo study were estimated assuming three different conditions, namely: (i) the diastolic configuration is neither preloaded nor axially stretched; (ii) the diastolic configuration is preloaded but not axially stretched; (iii) the diastolic configuration is preloaded and 5% axially stretched; and (iv) the diastolic configuration is preloaded and 10% axially stretched. The choice for the last two cases is based on the experimental observations of Holzapfel et al. (2005) where it is reported a physiological range for axial stretch in coronary arteries ranging between 5 and 10%.

The geometrical model and the optical flow u^OF used for this study are the ones previously presented in Figure 8. The geometric model was partitioned in sextants with a unique concentric layer leading to the estimation of 6 material parameters (the same partition used in Figure 12. The remaining parameters for the mechanical problems and data assimilation process are described in Table 1 along with the estimated values c_i. The results showed different trends for soft (c < 200 kPa, i.e., c₁, c₂, c₃, and c₅) and stiff materials (c ≥ 200 kPa, i.e., c₄ and c₆). The obtained parameter c increases in the soft tissues and decreases in stiff tissues as the baseline stress increases from a preload-free to a preloaded state. Further increments in the baseline stress due to the axial stretch result in material stiffening for these two categories of tissues. Interestingly, the increment of the parameter uncertainty σ_θ or the decrement of the observation uncertainty σ_Z in the stiff tissues increases the estimate of parameter c more in the preload-free state than in the preloaded cases. In fact, in the cases 2 and 3, the preloaded and 5% axially stretched model (cases 2.C and 3.C) featured lower c values in the stiff tissues than the preload-free model (cases 2.A and 3.A), contrarily to case 1. Some of these findings may appear counter-intuitive at first glance because as the baseline stress state increases it would be expected that all tissues soften to maintain the same deformation for the given load. Thus, the following paragraphs address the role of assimilation uncertainties, image artifacts, and the very mechanical model in the assimilation.

TABLE 1

Table 1. Sensitivity of the parameter estimation with respect to the baseline stress conditions and uncertainty parameters.

Firstly, as the baseline stress at the diastolic configuration rises, the parameter estimation is less sensitive with respect to variations between the predicted and the observed displacements i.e., Z−Ẑ_k. For the different baseline stress states, it was assumed the same observation uncertainty which is analog to establish an uncertainty interval for the observed strains. As the Neo-Hookean model consists of a quadratic stress-strain relation, the increment of the baseline stress yields an increase in the uncertainty interval of the stresses. And because the stress is linear to the material parameter c, the estimated parameters undergo the same increase of their uncertainties diminishing the accuracy of their estimation. Moreover, the estimated value of c increased as the baseline state is subjected to a more significant preload condition, turning the data assimilation process even less sensitive. In short, this implies that dealing with the real problem –for which preload is definitely a condition of the vessels– is even more challenging than the case where initial stress conditions are neglected.

Secondly, the gap between the observations and the predicted displacements, hereafter simply discrepancy, in the data assimilation process is in part given by some observed displacement components generated by errors in the image processing stage and by physical phenomena which is not recoverable by the proposed mechanical and material models (e.g., external tissues, off-plane displacements, compressible materials or, even, misrepresentation of the constitutive law). These discrepancies could be referred to as out-of-model components, introducing a bias in the predicted displacement field and in the parameter estimation as well. Comparing the estimations with different baseline assumptions, it is observed that the discrepancies of the identified parameter value remain below 37% and 10% for soft and stiff tissues, respectively. Particularly, we choose to use the more complex model (preloaded and axially stretched) in the following in-vivo studies because it endows the mechanical setting with more relevant physical features when compared to the other models.

3.4. In-Vivo Cases

The proposed methodology is now applied to 4 in-vivo cases featuring atherosclerotic lesions to derive their specific mechanical models. The goal is to analyze the accuracy of the mechanical models to predict the optical flow observations, as well as, to assess the usage of multiple (more than two) cardiac phases (and then more than one optical flow displacement field as observational data) in the parameter estimation. For each lesion, the IVUS frames that are involved in the data assimilation correspond to end-diastole, 50%-systole and full-systole, as dictated by the ECG signal of the IVUS study. Optical flow was estimated between end-diastole and 50%-systole frames and end-diastole and full-systole frames, denoted by ${\text{u}}_1^{\text{OF}}$ and ${\text{u}}_2^{\text{OF}}$ respectively (see Figure 12). Then, we compare the resulting estimated parameter for two cases: when the assimilation is performed using a single optical flow displacement field as observation ($Z={\left[{\text{u}}_2^{\text{OF}}\right]}^T$); and when two optical flow displacement fields are utilized as observations ($Z={\left[{\text{u}}_1^{\text{OF}},{\text{u}}_2^{\text{OF}}\right]}^T$). Note that the observed displacement field for maximum load, i.e., ${\text{u}}_2^{\text{OF}}$, is employed in both cases because the displacement between end-diastole and systole is expected to yield higher strains.

The geometric model was partitioned in sextants with a unique concentric layer (see Figure 11). Each partition contains only a single type of material leading to a data assimilation process with 6 material parameters. The diastolic configuration is preloaded and 10% axially stretched for all cases and the blood pressure at each phase was assumed to be 80, 100, and 120 mmHg for the end-diastole, 50%-systole and full-systole, respectively. The parameter τ was set to 100 for lesions 1, 3, and 4 and 50 in case 2, the latter avoided contact at the luminar surfaces during in the preload problem. The ROUKF uncertainties were fixed to σ_θ = 1 and ${\sigma }_Z=10^{-2}$ mm.

FIGURE 11

Figure 11. Optical flow estimated between end-diastole and 50%-systole frames, ${\text{u}}_1^{\text{OF}}$, and end-diastole and full-systole frames, ${\text{u}}_2^{\text{OF}}$ for 4 in-vivo atherosclerotic lesions: (row 1) IVUS end-diastolic frame depicting the atherosclerotic lesions; (row 2) displacement field ${\text{u}}_1^{\text{OF}}$; (row 3) displacement field ${\text{u}}_2^{\text{OF}}$.

The proposed data assimilation process rendered the results depicted in Figure 12. The material parameters estimated in all cases remained within the physiological range (between 1 kPa to 10 MPa, see Walsh et al., 2014). Also, the addition of an extra displacement field as observation showed no considerable effect for cases 2 and 3. The reliability of the results can be assessed in terms of the model prediction error presented in Table 2. Due to intrinsic sources of errors in the observations (motion artifacts, spatial incoherence between cross-sections in the cardiac cycle and optical flow model artifacts), it is expected an observation error of few pixel spacing units (recall that the image discretization spacing is 16μm). Thus, model prediction errors for cases 3 and 4, and even case 1 for a single optical flow field per cardiac phase, seems highly reliable in terms of our observation precision since the error results 26±14μm (1.625±0.875 pixel spacing units), while case 2 seems to be the less reliable estimation with an average error of 43±24μm. Overall, the average model prediction error was below 43μm and 61μm for the observation with one or two observational data, respectively.

FIGURE 12

Figure 12. Parameter estimation and discrepancies between the model prediction and observations for 4 in-vivo IVUS frames featuring atherosclerotic plaques (one per column). (Row 1) Estimation of parameters c_i using $Z={\left[{\text{u}}_2^{\text{OF}}\right]}^T$; (row 2) ε_r for the observations related to ${\text{u}}_2^{\text{OF}}$; (row 3) Estimation of parameters c_i using $Z={\left[{\text{u}}_1^{\text{OF}},{\text{u}}_2^{\text{OF}}\right]}^T$; (row 4) ε_r for the observations related to ${\text{u}}_1^{\text{OF}}$; (row 5) ε_r for the observations related to ${\text{u}}_2^{\text{OF}}$. The relative displacement discrepancy between the model predictions and the observations was defined as ε_r = ∥Ẑ_k−Z∥/〈∥Z∥〉 where Ẑ_k are the model predictions at the last Kalman iteration, Z are the optical flow observations and 〈·〉 denotes the mean value in Ω_s.

TABLE 2

Table 2. Model prediction error after data assimilation process for the 4 in-vivo cases using 1 or 2 loading conditions.

In case 1, it is observed that the material parameters estimated with 1 and 2 optical flow displacement data are significantly different. The flow ${\text{u}}_1^{\text{OF}}$ presents larger displacements than ${\text{u}}_2^{\text{OF}}$, which seems counter-intuitive since the blood pressure variation is smaller for the former condition. However, the motion exerted by the cardiac contraction is higher, in fact the larger component of displacement is rigid (a rotation of the structures). Thus, as ${\text{u}}_1^{\text{OF}}$ presents the observation components with higher norm, it features a larger contribution than ${\text{u}}_2^{\text{OF}}$ during the data assimilation process (see Equation 30). In that manner, the parameters estimated with 1 flow datum minimize discrepancies against ${\text{u}}_2^{\text{OF}}$ while the ones estimated with 2 flow data minimize mainly discrepancies against ${\text{u}}_1^{\text{OF}}$.

Conversely in cases 2 and 3, the observation ${\text{u}}_1^{\text{OF}}$ is the one with smaller displacements (≈4 and 2 times smaller for cases 2 and 3, respectively), yielding a small contribution to the data assimilation. This implies that the minimization of the discrepancies between the model predictions and the observations (i.e., Ẑ_k−Z) related to ${\text{u}}_2^{\text{OF}}$ dominates over the discrepancies associated to ${\text{u}}_1^{\text{OF}}$. In fact, Figure 12 shows that the discrepancies represented by ε_r for ${\text{u}}_2^{\text{OF}}$ remained almost invariant using 1 or 2 flows in the observation.

In case 4, the discrepancies between the model predictions and the observations related ${\text{u}}_2^{\text{OF}}$ also remained invariant using 1 and 2 flows data in the observation, although the parameters estimated in the lower part of the geometry varied significantly (see Figure 12). As previously studied in section 3.2.2, the guidewire artifact in the lower part of these images features a swinging movement not related with the arterial-wall motion. To approximate the artifact's rigid motion, the local tissue is stiffened during the assimilation process when 1 single flow was employed as observation. Conversely when ${\text{u}}_1^{\text{OF}}$ is added to the observations, the spurious motion of the guidewire is negligible, and the data assimilation is not affected by this artifact.

To determine the applicability of the current approach in clinical practice, we execute the in-vivo cases in a single server with 2 Intel Xeon CPU E5-2620 at 2.00 GHz processor (each with 12 threads) and Kingston 99U5471-031.A00LF at 1333 MHz (latency of 27 ns) RAM memory. For data assimilation of these in-vivo cases, mesh (3 threads per mechanical problem) and sigma parallelism (1 thread per sigma point) were applied because it delivered the best speed up for our 24 threads (actually only 21 were employed). The wall clock time reported in Table 2 for each execution showed that the current methodology is appropriate for offline medical applications because the processing times elapsed from 0.5 to 3 days. The use of clusters would allow further processing speed up exploiting the load parallelism as well as a more massive parallelization at the mesh level.

4. Discussions

The presented methodology offers a workflow to estimate material parameters for mechanical model of coronary arteries. The strategy is composed by three key components: the image processing, the mechanical model and the data assimilation algorithm. The most appealing aspect of this proposal is that the three components are loosely coupled as black boxes which allowed us to modify, as required, each component without the need for altering the remaining ones. In fact, the image processing renders observations for the data assimilation, regardless the imaging technique employed and the nature of the displacement field. In turn, the mechanical model can also be modified without influencing in the other components, it simply must receive a set of parameters and return back the internal state variables to the data assimilation strategy. Due to this architectural design, this initial biomechanical characterization approach can be further refined by improving aspects of these individual components. Some identified hot-spots for improvement are discussed in what follows.

The data assimilation showed high sensitivity with respect to variations in the model boundary conditions which aimed at mimicking the external tissues. As the displacement over the boundary was increasingly constrained (large τ) the model was less sensitive to variations in the material parameters, hindering the parameter estimation and, even, causing divergence of the Kalman iterative process in some situations. Also, the disagreement in the spatial arrangement of model forces and the in-vivo (unknown) forces at the boundary notoriously affects the outcome of the estimation. This was exposed in section 3.2.2 when an image artifact (the IVUS guidewire) induced a spurious tangential displacement in the observation and the boundary condition. It was also showed that if a homogeneous Neumann condition is assumed at the site of such artifact, the parameter estimation varies significantly (from 4 to 15 fold reduction of parameter c). Improving the capabilities of the model in this sense requires to incorporate the estimation of these forces exerted by external tissues in the data assimilation process. In short, parameter τ could be a further variable to be estimated.

It is also important to highlight that this approach can be directly extended to account for more geometrically and physically complex models. The set of here reported results constitute a solid proof of concept toward the extension of this methodology. Here, we derived a patient-specific mechanical model for an orthogonal slice of the vessel assuming plane strain state with an homogeneous axial traction force. However, there are some assumptions that imply neglecting certain physical components that may be necessary to increase the accuracy of the estimated stress/strain state of the vessel. To list some of them: (i) shear forces exerted by the blood flow which are expected to be key in the study of plaque development (Stone et al., 2003; Chatzizisis et al., 2008); (ii) out-of-plane forces produced by the blood pressure due to the heterogeneous constitution of the vessel wall and the tilting of the transducer tip with respect to the cross-section; and (iii) variable axial tractions along the cross-section due to the heterogeneous composition of the vessel wall. These issues can be tackled at once by making use of 3D models. In fact, the image processing strategy allows the gating and registration of the whole arterial 3D volume of the study. Also, the extension of the optical flow techniques to 3D domains is straightforward by a proper adaptation of the differential operators and Gaussian kernel within the formulation. A further issue to address is the spatial reconstruction to obtain the proper 3D geometrical description of the vessel instead of its rectified representation in intrinsic coordinates delivered by the IVUS study. The integration of IVUS with angiographic images enabled us to perform such 3D reconstruction, as reported in Maso Talou (2013). These extensions imply in heavier computational cost and complementary implementation aspects, yet, they present no further conceptual differences regarding the methodology presented in this work.

Extension to 3D problems discussed above, as said, becomes computationally more demanding. Associated to the image processing, the cost scales with the number of cross-sections extracted from the IVUS dataset. However, the registration stage, which is the most computationally intensive task, is fully parallel (see Maso Talou et al., 2017) and gating cost is negligible. Thus, the performance of the optical flow and the spatial reconstruction process through the integration with angiographic images, turn out to be key for the efficiency of the methodology in 3D cases. Regarding the data assimilation procedure, the computational cost continues to be the approximate solution of the mechanical problem. As a significant increase in the number of degrees of freedom is expected, the computational cost would raise as well.

A first limitation in the present scheme is that the displacement field retrieved from medical images is naively used as observation from our model without further processing. This implies that the performance of the method can be improved by extracting the observation components that are spurious (such as artifacts or unreliable regions of the optical flow displacement field) or even incompatible with our model (e.g., use only the divergence-free component of the field because the mechanical model is incompressible).

Regarding the baseline stress state in our model, the residual stresses produced during the arterial tissue genesis and growth have clearly been neglected. In Wang et al. (2017), an experimental test showed that the omission of these residual stresses may produce a significant overestimation of internal stresses (from 2- to 4-fold the actual stress). Furthermore, it has been observed (Guo et al., 2017) that accounting for residual stresses is also relevant for the proper material parameter estimation. This seems to be natural, as residual stresses can be considered as a subproduct of existing residual deformations (that, in general, may not be kinematically compatible, i.e., they cannot be derived from continuous displacement fields) of the elastin matrix. Consequently, not only the stresses are not properly assessed, but the actual deformations observed at the equilibrium states are misguided. These facts highlight the need for further research to tackle simultaneously the estimation of both material parameters and residual deformations in arterial walls. In recent works (Ares, 2016; Ares et al., 2017), models and methods for the estimation of such residual stresses were proposed, with a similar spirit to the one developed in this work.

At last, it is worthwhile to remark that no validation techniques are currently available for the assessment of stress-strain state in in-vivo conditions. Even though, approaches for an indirect in-vivo or ex-vivo validation can be discussed. Techniques such as elastography and palpography (Ophir et al., 1991; Shapo et al., 1996; Céspedes et al., 1997; de Korte et al., 1998; Céspedes et al., 2000) deliver with some degree of reliability the stresses in the innermost part of the vessel. In these cases a Bland-Altman analysis can be applied to assess the similarity between the prediction of our approach and elastographic solutions. A more controlled experimental setup can be planned for ex-vivo condition using coronary specimens. For each specimens, an IVUS study can be acquired and a specimen-specific model can be constructed employing the proposed methodology. Finally, several mechanical tests can be carried out with the specimens comparing their mechanical response with predictions given by our specimen-specific models. Another in-vivo alternative is to associate ranges of the estimated material parameters to the underlying tissue composition, delivering a histological description of the vessel (usually referred to as virtual histology). As there are already methods that estimate the vessel histology from IVUS images (e.g., Kawasaki et al., 2002; Nair et al., 2002; Sathyanarayana et al., 2009), a comparative analysis can be performed to evaluate the degree of agreement between the proposed method and these virtual histologies. An appealing aspect of this last validation is that the techniques presented in those works are already validated with cadaveric specimens of coronary arteries. The experimental settings suggested above should serve to bridge the world of computational models and methods with the experimental realm, toward gaining insight into the complex mechanisms underlying the development of cardiovascular diseases.

5. Final Remarks

A data assimilation environment for analysis of arterial models and material characterization was described. The proposed methodology delivers the necessary tools to construct patient-specific mechanical models of an arterial site using data from standard IVUS studies. A complete sensitivity analysis of the biomechanical characterization with respect to numerical and physical parameters was reported to aid the methodology setup, as well as the interpretation of data assimilation outcomes. Validation in controlled scenarios was provided to demonstrate the capabilities of the present approach.

The potential and limitations of this approach were exposed and discussed in the previous section, delineating future research to enhance the image processing stage and the mechanical model of the arterial wall for this problem.

The applicability of this methodology on in-vivo scenarios was proven in the characterization of the arterial tissue for 4 in-vivo atherosclerotic lesions. After data assimilation, the obtained mechanical models predicted the displacement field between diastole and systole with errors below 43μm using frames of only two cardiac phases. Although no validation was performed with the in-vivo cases, the estimated material parameters remained within the expected range for this kind of tissue.

The development of this tool for the biomechanical analysis allows the indirect estimation of the internal stress state of the arterial wall. Such information combined with the vessel histology (that can be inferred from the material parameters here estimated) enables the assessment of the structural integrity of the atherosclerotic plaque to aid medical decisions and research. In summary the proposed strategy provides an imaging-assimilation-mechanics integrated environment to characterize, within a truly in-vivo and patient-specific setting, the behavior of the materials that compose the arterial vessels, specifically coronary vessels, which is of the utmost importance in assessing risk of plaque progress and rupture.

Author Contributions

GM, GA, PB, and RF designed the model and the computational framework. GM and GA carried out the computational implementation. GM, PB, and RF planned the experiments. GM and PB performed the calculations and wrote the manuscript with input from all authors. CG and PL performed measurements and contributed to sample preparation. All authors contributed to the discussion and interpretation of the results, and helped shape the final version of the manuscript.

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.

Acknowledgments

This work was partially supported by the Brazilian agencies CNPq, FAPERJ, and CAPES. The support of these agencies is gratefully acknowledged.

Footnotes

1. ^The reparametrization $2^{\widehat{\theta }}$ modifies the assumed distribution of the parameter from normal distribution to log-normal distribution. Mean and covariance of θ are propagated in a different fashion than mean and covariance of c. Nevertheless, there is a similar statistical interpretation for c using these descriptors in an exponential space of coordinates. For example, a covariance of σ = 1 must be understood as giving the same probability of c being half or twice its initial value.

References

Arbab-Zadeh, A., DeMaria, A. N., Penny, W. F., Russo, R. J., Kimura, B. J., and Bhargava, V. (1999). Axial movement of the intravascular ultrasound probe during the cardiac cycle: implications for three-dimensional reconstruction and measurements of coronary dimensions. Amer. Heart J. 138(5 Pt 1), 865–872. doi: 10.1016/S0002-8703(99)70011-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ares, G. D. (2016). Integrative Computational Modeling & In-vivo Characterization of Residual Deformations in Hemodynamics. Ph.D. thesis, National Laboratory for Scientific Computing.

Ares, G. D., Blanco, P. J., Urquiza, S. A., and Feijóo, R. A. (2017). Identification of residual stresses in multi-layered arterial wall tissues using a variational framework. Comput. Methods Appl. Mech. Eng. 319, 287–313. doi: 10.1016/j.cma.2017.02.028

CrossRef Full Text | Google Scholar

Avril, S., Badel, P., and Duprey, A. (2010). Anisotropic and hyperelastic identification of in vitro human arteries from full-field optical measurements. J. Biomech. 43, 2978–2985. doi: 10.1016/j.jbiomech.2010.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bentzon, J. F., Otsuka, F., Virmani, R., and Falk, E. (2014). Mechanisms of plaque formation and rupture. Circ. Res. 114, 1852–1866. doi: 10.1161/CIRCRESAHA.114.302721

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertagna, L., and Veneziani, A. (2014). A model reduction approach for the variational estimation of vascular compliance by solving an inverse fluid–structure interaction problem. Inverse Prob. 30:055006. doi: 10.1088/0266-5611/30/5/055006

CrossRef Full Text | Google Scholar

Bertoglio, C., Barber, D., Gaddum, N., Valverde, I., Rutten, M., Beerbaum, P., et al. (2014). Identification of artery wall stiffness: In vitro validation and in vivo results of a data assimilation procedure applied to a 3d fluid–structure interaction model. J. Biomech. 47, 1027–1034. doi: 10.1016/j.jbiomech.2013.12.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertoglio, C., Moireau, P., and Gerbeau, J.-F. (2012). Sequential parameter estimation for fluid–structure problems: Application to hemodynamics. Int. J. Numer. Methods Biomed. Eng. 28, 434–455. doi: 10.1002/cnm.1476

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanco, P. J., Ares, G. D., Urquiza, S. A., and Feijóo, R. A. (2016). On the effect of preload and pre-stretch on hemodynamic simulations: an integrative approach. Biomech. Model. Mechanobiol. 15, 593–627. doi: 10.1007/s10237-015-0712-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanco, P. J., Watanabe, S. M., Passos, M. A. R., Lemos, P. A., and Feijóo, R. A. (2015). An anatomically detailed arterial network model for one-dimensional computational hemodynamics. IEEE Trans. Biomed. Eng. 62, 736–753. doi: 10.1109/TBME.2014.2364522

PubMed Abstract | CrossRef Full Text | Google Scholar

Brox, T., Bruhn, A., Papenberg, N., and Weickert, J. (2004). “High accuracy optical flow estimation based on a theory for warping,” in Computer Vision-ECCV 2004 (Prague: Springer), 25–36. doi: 10.1007/978-3-540-24673-2_3

CrossRef Full Text | Google Scholar

Bulant, C. A., Blanco, P. J., Maso Talou, G. D., Guedes Bezerra, C., Lemos, P. A., and Feijóo, R. A. (2017). A head-to-head comparison between ct-and ivus-derived coronary blood flow models. J. Biomech. 51, 65–76. doi: 10.1016/j.jbiomech.2016.11.070

PubMed Abstract | CrossRef Full Text | Google Scholar

Caiazzo, A., Caforio, F., Montecinos, G., Muller, L. O., Blanco, P. J., and Toro, E. F. (2017). Assessment of reduced-order unscented kalman filter for parameter identification in 1-dimensional blood flow models using experimental data. Int. J. Numer. Methods Biomed. Eng. 33:e2843. doi: 10.1002/cnm.2843

CrossRef Full Text | Google Scholar

Céspedes, E., de Korte, C., van der Steen, A., von Birgelen, C., and Lance, C. (1997). Intravascular elastography: principles and potentials. Semin. Int. Cardiol. 2, 55–62.

PubMed Abstract | Google Scholar

Céspedes, E. I., de Korte, C. L., and van der Steen, A. F. (2000). Intraluminal ultrasonic palpation: assessment of local and cross-sectional tissue stiffness. Ultrasound Med. Biol. 26, 385–396. doi: 10.1016/S0301-5629(99)00169-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, T. F., Golub, G. H., and Mulet, P. (1999). A nonlinear primal-dual method for total variation-based image restoration. SIAM J. Sci. Comput. 20, 1964–1977. doi: 10.1137/S1064827596299767

CrossRef Full Text | Google Scholar

Chapelle, D., Moireau, P., and Le Tallec, P. (2009). Robust filtering for joint state parameter estimation for distributed mechanical systems. Discrete Contin. Dyn. Syst. Ser. A 23, 65–84. doi: 10.3934/dcds.2009.23.65

CrossRef Full Text | Google Scholar

Chatzizisis, Y. S., Jonas, M., Coskun, A. U., Beigel, R., Stone, B. V., Maynard, C., et al. (2008). Prediction of the localization of high-risk coronary atherosclerotic plaques on the basis of low endothelial shear stress: an intravascular ultrasound and histopathology natural history study. Circulation 117, 993–1002. doi: 10.1161/CIRCULATIONAHA.107.695254

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, B., and Dinstein, I. (2002). New maximum likelihood motion estimation schemes for noisy ultrasound images. Patt. Recogn. 35, 455–463. doi: 10.1016/S0031-3203(01)00053-X

CrossRef Full Text | Google Scholar

Conway, C., McGarry, P., and McHugh, P. E. (2014). Modelling of atherosclerotic plaque for use in a computational test-bed for stent angioplasty. Ann. Biomed. Eng. 42, 2425–2439. doi: 10.1007/s10439-014-1107-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Conway, C., Sharif, F., McGarry, P., and McHugh, P. E. (2012). A computational test-bed to assess coronary stent implantation mechanics using a population-specific approach. Cardiovasc. Eng. Technol. 3, 374–387. doi: 10.1007/s13239-012-0104-8

CrossRef Full Text | Google Scholar

Crea, F., and Liuzzo, G. (2013). Pathogenesis of acute coronary syndromes. J. Amer. Coll. Cardiol. 61, 1–11. doi: 10.1016/j.jacc.2012.07.064

PubMed Abstract | CrossRef Full Text | Google Scholar

de Korte, C. L., Cespedes, E., van der Steen, A. F., Pasterkamp, G., and Bom, N. (1998). Intravascular ultrasound elastography: assessment and imaging of elastic properties of diseased arteries and vulnerable plaque. Eur. J. Ultras. 7, 219–224. doi: 10.1016/S0929-8266(98)00043-3

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Elia, M., Mirabella, L., Passerini, T., Perego, M., Piccinelli, M., Vergara, C., et al. (2012). “Applications of variational data assimilation in computational hemodynamics,” in Modeling of Physiological Flows, eds D. Ambrosi, A. Quarteroni, and G. Rozza (Milano: Springer), 363–394.

Google Scholar

Guo, X., Zhu, J., Maehara, A., Monoly, D., Samady, H., Wang, L., et al. (2017). Quantify patient-specific coronary material property and its impact on stress/strain calculations using in vivo ivus data and 3d fsi models: a pilot study. Biomech. Model. Mechanobiol. 16, 333–344. doi: 10.1007/s10237-016-0820-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Holzapfel, G. A., Gasser, T. C., and Ogden, R. W. (2000). A new constitutive framework for arterial wall mechanics and a comparative study of material models. J. Elast. Phys. Sci. Solids 61, 1–48. doi: 10.1023/A:1010835316564

CrossRef Full Text | Google Scholar

Holzapfel, G. A., Sommer, G., Gasser, C. T., and Regitnig, P. (2005). Determination of layer-specific mechanical properties of human coronary arteries with nonatherosclerotic intimal thickening and related constitutive modeling. Amer. J. Physiol. Heart Circul. Physiol. 289, H2048–H2058. doi: 10.1152/ajpheart.00934.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Horn, B. K. P., and Schunck, B. G. (1981). Determining optical flow. Artif. Intell. 17, 185–203. doi: 10.1016/0004-3702(81)90024-2

CrossRef Full Text | Google Scholar

Julier, S. J. (2003). “The spherical simplex unscented transformation,” in Proceedings of the American Control Conference, Vol. 3 (Denver, Co: IEEE), 2430–2434. doi: 10.1109/ACC.2003.1243439

CrossRef Full Text | Google Scholar

Julier, S. J., and Uhlmann, J. K. (2002). “Reduced sigma point filters for the propagation of means and covariances through nonlinear transformations,” in Proceedings of the 2002 American Control Conference (IEEE Cat. No. CH37301), vol. 2 (Anchorage, AK: IEEE), 887–892. doi: 10.1109/ACC.2002.1023128

CrossRef Full Text | Google Scholar

Julier, S. J., and Uhlmann, J. K. (2004). Unscented filtering and nonlinear estimation. Proc. IEEE 92, 401–422. doi: 10.1109/JPROC.2003.823141

CrossRef Full Text | Google Scholar

Karimi, A., Navidbakhsh, M., and Shojaei, A. (2015). A combination of histological analyses and uniaxial tensile tests to determine the material coefficients of the healthy and atherosclerotic human coronary arteries. Tissue Cell 47, 152–158. doi: 10.1016/j.tice.2015.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Karypis, G., and Kumar, V. (1998). Multilevelk-way partitioning scheme for irregular graphs. J. Parall. Distrib. Comput. 48, 96–129. doi: 10.1006/jpdc.1997.1404

CrossRef Full Text | Google Scholar

Kawasaki, M., Takatsu, H., Noda, T., Sano, K., Ito, Y., Hayakawa, K., et al. (2002). In vivo quantitative tissue characterization of human coronary arterial plaques by use of integrated backscatter intravascular ultrasound and comparison with angioscopic findings. Circulation 105, 2487–2492. doi: 10.1161/01.CIR.0000017200.47342.10

PubMed Abstract | CrossRef Full Text | Google Scholar

Lagrée, P.-Y. (2000). An inverse technique to deduce the elasticity of a large artery. Eur. Phys. J. Appl. Phys. 9, 153–163. doi: 10.1051/epjap:2000213

CrossRef Full Text | Google Scholar

Lassila, T., Manzoni, A., Quarteroni, A., and Rozza, G. (2013). A reduced computational and geometrical framework for inverse problems in hemodynamics. Int. J. Numer. Methods Biomed. Eng. 29, 741–776. doi: 10.1002/cnm.2559

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X. S., and Demmel, J. W. (2003). SuperLU_DIST: a scalable distributed-memory sparse direct solver for unsymmetric linear systems. ACM Trans. Math. Softw. 29, 110–140. doi: 10.1145/779359.779361

CrossRef Full Text | Google Scholar

Lipponen, A., Seppänen, A., and Kaipio, J. P. (2010). Reduced-order estimation of nonstationary flows with electrical impedance tomography. Inverse Probl. 26:074010. doi: 10.1088/0266-5611/26/7/074010

CrossRef Full Text | Google Scholar

Loree, H. M., Grodzinsky, A. J., Park, S. Y., Gibson, L. J., and Lee, R. T. (1994). Static circumferential tangential modulus of human atherosclerotic tissue. J. Biomech. 27, 195–204. doi: 10.1016/0021-9290(94)90209-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, V., Clément, F., Decoene, A., and Gerbeau, J.-F. (2005). “Parameter identification for a one-dimensional blood flow model,” in ESAIM: Proceedings, vol. 14 (Marseille: EDP Sciences), 174–200.

Google Scholar

Maso Talou, G. D. (2013). Ivus Images Segmentation Driven by Active Contours and Spacio-temporal Reconstruction of the Coronary Vessels Aided by Angiographies. Master's thesis, National Laboratory for Scientific Computing.

Maso Talou, G. D. (2017). From Medical Image Processing to In-vivo Mechanical Characterization: A Framework Based on IVUS Studies. Ph.D. thesis, Laboratório Nacional de Computação Científica.

Maso Talou, G. D., Blanco, P. J., Larrabide, I., Bezerra, C. G., Lemos, P. A., and Feijóo, R. A. (2017). Registration methods for IVUS: transversal and longitudinal transducer motion compensation. IEEE Trans. Biomed. Eng. 64, 890–903. doi: 10.1109/TBME.2016.2581583

CrossRef Full Text | Google Scholar

Maso Talou, G. D., Larrabide, I., Blanco, P. J., Bezerra, C. G., Lemos, P. A., and Feijóo, R. A. (2015). Improving cardiac phase extraction in ivus studies by integration of gating methods. IEEE Trans. Biomed. Eng. 62, 2867–2877. doi: 10.1109/TBME.2015.2449232

PubMed Abstract | CrossRef Full Text | Google Scholar

Mathers, C., Stevens, G., Mahanani, W. R., Ho, J., Fat, D. M., and Hogan, D. (2016). Who Methods and Data Sources for Country-level Causes of Death 2000-2015. Technical Report, World Health Organization.

Milner, J. S., Moore, J. A., Rutt, B. K., and Steinman, D. A. (1998). Hemodynamics of human carotid artery bifurcations: computational studies with models reconstructed from magnetic resonance imaging of normal subjects. J. Vasc. Surg. 28, 143–156. doi: 10.1016/S0741-5214(98)70210-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyashita, H. (2012). Clinical assessment of central blood pressure. Curr. Hypertens. Rev. 8, 80–90. doi: 10.2174/157340212800840708

PubMed Abstract | CrossRef Full Text | Google Scholar

Moireau, P., and Chapelle, D. (2011). Reduced-order unscented kalman filtering with application to parameter identification in large-dimensional systems. ESAIM Control Optim. Calc. Variat. 17, 380–405. doi: 10.1051/cocv/2010006

CrossRef Full Text | Google Scholar

Moireau, P., Chapelle, D., and Le Tallec, P. (2008). Joint state and parameter estimation for distributed mechanical systems. Comput. Methods Appl. Mech. Eng. 197, 659–677. doi: 10.1016/j.cma.2007.08.021

CrossRef Full Text | Google Scholar

Moireau, P., Chapelle, D., and Le Tallec, P. (2009). Filtering for distributed mechanical systems using position measurements: perspectives in medical imaging. Inverse Prob. 25:035010. doi: 10.1088/0266-5611/25/3/035010

CrossRef Full Text | Google Scholar

Nair, A., Kuban, B. D., Tuzcu, E. M., Schoenhagen, P., Nissen, S. E., and Vince, D. G. (2002). Coronary plaque classification with intravascular ultrasound radiofrequency data analysis. Circulation 106, 2200–2206. doi: 10.1161/01.CIR.0000035654.18341.5E

PubMed Abstract | CrossRef Full Text | Google Scholar

O'brien, E., Waeber, B., Parati, G., Staessen, J., and Myers, M. G. (2001). Blood pressure measuring devices: recommendations of the european society of hypertension. Brit. Med. J. 322:531. doi: 10.1136/bmj.322.7285.531

PubMed Abstract | CrossRef Full Text | Google Scholar

Ophir, J., Céspedes, I., Ponnekanti, H., Yazdi, Y., and Li, X. (1991). Elastography: a method for imaging the elasticity in biological tissues. Ultrasonic Imaging 13, 111–134. doi: 10.1177/016173469101300201

PubMed Abstract | CrossRef Full Text | Google Scholar

Perego, M., Veneziani, A., and Vergara, C. (2011). A variational approach for estimating the compliance of the cardiovascular tissue: An inverse fluid-structure interaction problem. SIAM J. Sci. Comput. 33, 1181–1211. doi: 10.1137/100808277

CrossRef Full Text | Google Scholar

Pérez Zerpa, J. M., and Canelas, A. (2016). Efficient formulations of the material identification problem using full-field measurements. Comput. Mech. 58, 1–21. doi: 10.1007/s00466-016-1291-1

CrossRef Full Text | Google Scholar

Rudin, L. I., Osher, S., and Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms. Phys. D 60, 259–268. doi: 10.1016/0167-2789(92)90242-F

CrossRef Full Text | Google Scholar

Sathyanarayana, S., Carlier, S., Li, W., and Thomas, L. (2009). Characterisation of atherosclerotic plaque by spectral similarity of radiofrequency intravascular ultrasound signals. EuroIntervention 5, 133–139. doi: 10.4244/EIJV5I1A21

PubMed Abstract | CrossRef Full Text | Google Scholar

Sermesant, M., Moireau, P., Camara, O., Sainte-Marie, J., Andriantsimiavona, R., Cimrman, R., et al. (2006). Cardiac function estimation from mri using a heart model and data assimilation: advances and difficulties. Med. Image Anal. 10, 642–656. doi: 10.1016/j.media.2006.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Shapo, B. M., Crowe, J. R., Erkamp, R., Emelianov, S. Y., Eberle, M. J., and O'Donnell, M. (1996). Strain imaging of coronary arteries with intraluminal ultrasound: experiments on an inhomogeneous phantom. Ultrasonic Imaging 18, 173–191. doi: 10.1177/016173469601800302

CrossRef Full Text | Google Scholar

Stålhand, J. (2009). Determination of human arterial wall parameters from clinical data. Biomech. Model. Mechanobiol. 8, 141–148. doi: 10.1007/s10237-008-0124-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Stone, P. H., Coskun, A. U., Kinlay, S., Clark, M. E., Sonka, M., Wahle, A., et al. (2003). Effect of endothelial shear stress on the progression of coronary artery disease, vascular remodeling, and in-stent restenosis in humans: in vivo 6-month follow-up study. Circulation 108, 438–444. doi: 10.1161/01.CIR.0000080882.35274.AD

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, C. A., Fonte, T. A., and Min, J. K. (2013). Computational fluid dynamics applied to cardiac computed tomography for noninvasive quantification of fractional flow reserve. J. Amer. Coll. Cardiol. 61, 2233–2241. doi: 10.1016/j.jacc.2012.11.083

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, C. A., Hughes, T. J., and Zarins, C. K. (1999). Effect of exercise on hemodynamic conditions in the abdominal aorta. J. Vasc. Surg. 29, 1077–1089. doi: 10.1016/S0741-5214(99)70249-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Torii, R., Oshima, M., Kobayashi, T., Takagi, K., and Tezduyar, T. E. (2007). Numerical investigation of the effect of hypertensive blood pressure on cerebral aneurysmdependence of the effect on the aneurysm shape. Int. J. Numer. Methods Fluids 54, 995–1009. doi: 10.1002/fld.1497

CrossRef Full Text | Google Scholar

Wachinger, C., Wein, W., and Navab, N. (2008). Registration strategies and similarity measures for three-dimensional ultrasound mosaicing. Acad. Radiol. 15, 1404–1415. doi: 10.1016/j.acra.2008.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Wahle, A., Wellnhofer, E., Mugaragu, I., Saner, H., Oswald, H., and Fleck, E. (1995). Assessment of diffuse coronary artery disease by quantitative analysis of coronary morphology based upon 3-d reconstruction from biplane angiograms. IEEE Trans. Med. Imaging 14, 230–241. doi: 10.1109/42.387704

PubMed Abstract | CrossRef Full Text | Google Scholar

Walsh, M. T., Cunnane, E. M., Mulvihill, J. J., Akyildiz, A. C., Gijsen, F. J. H., and Holzapfel, G. A. (2014). Uniaxial tensile testing approaches for characterisation of atherosclerotic plaques. J. Biomech. 47, 793–804. doi: 10.1016/j.jbiomech.2014.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Zhu, J., Samady, H., Monoly, D., Zheng, J., Guo, X., et al. (2017). Effects of residual stress, axial stretch, and circumferential shrinkage on coronary plaque stress and strain calculations: a modeling study using ivus-based near-idealized geometries. J. Biomech. Eng. 139:014501. doi: 10.1115/1.4034867

PubMed Abstract | CrossRef Full Text | Google Scholar