An image-based computational hemodynamics study of the Systolic Anterior Motion of the mitral valve

Systolic Anterior Motion (SAM) of the mitral valve – often associated with Hypertrophic Obstructive Cardiomyopathy (HOCM) – is a cardiac pathology in which a functional subaortic stenosis is induced during systole by the mitral leaflets partially obstructing the outflow tract of the left ventricle. Its assessment by diagnostic tests is often difficult, possibly underestimating its severity and thus increasing the risk of heart failure. In this paper, we propose a new computational pipeline, based on cardiac cine Magnetic Resonance Imaging (cine-MRI) data, for the assessment of SAM. The pipeline encompasses image processing of the left ventricle and the mitral valve, and numerical investigation of cardiac hemodynamics by means of Computational Fluid Dynamics (CFD) in a moving domain with image-based prescribed displacement. Patient-specific geometry and motion of the left ventricle are considered in view of an Arbitrary Lagrangian–Eulerian approach for CFD, while the reconstructed mitral valve is immersed in the computational domain by means of a resistive method. We assess clinically relevant flow and pressure indicators in a parametric study for different degrees of SAM severity, in order to provide a better quantitative evaluation of the pathological condition. Moreover, we provide specific indications for its possible surgical treatment, i.e. septal myectomy.


Introduction
Systolic Anterior Motion (SAM) is a pathological condition of the heart where the anterior leaflet of the mitral valve is displaced towards the septum narrowing the Left Ventricle Outflow Tract (LVOT). This may lead to an unphysiologically high blood velocity through the LVOT and the aortic valve orifice, other than to an elevated intraventricular pressure gradient [1,2]. SAM is often related to Hypertrophic Cardiomyopathy (HCM), which is a pathological condition of the myocardium, characterized by marked myocardial hypertrophy (more than 15 mm of wall thickness), with a prevalence of the 0.2÷0.6% and an overall annual mortality rate of 1% in the western world. HCM is typically a genetic disorder [3], but the severity of its effects on the patient can be worsened by hypertension or other pathologies that provoke an increase of the ventricles afterload. When it affects the medio-basal portion of the septum, this condition takes the name of • the blood flow in the ventricle and ascending aorta in the presence of SAM will be computationally analyzed, in order to assess the hemodynamic effects of this condition; • preliminary clinical indications will be provided, that can help the surgeon in the preoperative design of septal myectomy; • in order to fulfill the points above, a novel reconstruction pipeline will be proposed, to include diagnostic imaging data in the computational model.
Computational studies regarding the mitral valve dynamics and its interaction with blood flow in realistic scenarios have been proposed since the early 2000's. Fluid-Structure Interaction (FSI) problems were considered for example in [9][10][11][12][13][14][15][16][17]. On the other side, computational hemodynamics where the domain displacement is prescribed -thus avoiding the complexity of FSI modeling for ventricle and valve dynamics -has shown its suitability in studying blood flow patterns, as done in [18][19][20][21]. In this direction, the last decade has seen computational hemodynamics based on kinetic medical images become an effective tool to provide quantitative insights about cardiovascular diseases and useful indications to design clinical practices. Among the others, we mention: the work [22], in which kinetic images of the aorta are used to reconstruct the motion of the vessel interface and to provide suitable boundary conditions at the vessel wall; the MRI-based study in [23], where a discussion on turbulence in the left ventricle is carried on; the works [24][25][26][27], where CFD in the left ventricle is driven by the motion reconstructed from three-dimensional (3D) echocardiography data, see also [28] for a review on ventricular hemodynamics. However, CFD studies prescribing cardiac valve kinematics directly reconstructed by medical images have been scarcely used so far to describe ventricular blood flow, see [25,29]. To the best of our knowledge, exploitation of cine-MRI in this context has not yet been performed, despite this technique is becoming increasingly important to study the cardiac function in clinics.
Concerning the analysis of blood flow in presence of SAM, patientspecific CFD studies are rare. At the best of our knowledge, this topic has been addressed so far only in [8] , which has considered an FSI model to inspect the effectiveness of septal myectomy, only in the early stage of systole. We point out that an FSI modeling of the mitral valve leaflets for the whole systolic phase would require to account for its complex structure, including the subvalvular apparatus. Furthermore, calibrating the valve rheological properties would represent a very challenging and subject-dependent task, especially in pathological conditions. On the other hand, diagnostic kinetic imaging can provide the leaflets motion without the need of a mechanical model: for this reason, we decide here to avoid an FSI description and, rather, to use such imaging data to prescribe the valve dynamics.
The aim of this work is to provide a CFD investigation and a numerical assessment of SAM based on kinetic imaging of both the Mitral Valve (MV) and the Left Ventricle (LV). In particular, the motion of the ventricle endocardium and of the mitral leaflets are reconstructed from cine-MRI data and employed to deform the computational domain and to prescribe boundary conditions to the blood flow model. In this way, the proposed computational framework: (i) avoids modeling the valve dynamics and myocardium mechanics, thus considerably reducing the computational effort; (ii) embeds patient-specific geometric and functional data without any prior assumption on the model parameters; (iii) avoids the modeling of the leaflets mechanics and the complex numerical treatment of an FSI problem; (iv) hinges upon standard image data routinely acquired in current clinical practice.
Towards this goal, we propose a complete pipeline for image processing, consisting in the segmentation of left ventricle and mitral valve, the reconstruction of the corresponding motion through registration algorithms, the fulfillment of missing data by means of suitable interpolation procedures, and the use of a template geometry of the human heart to integrate the anatomical data. Such a template allows completing the geometric data by including, e.g., the aorta and full 3D leaflets, that are usually not acquired during standard clinical imaging procedures.
We consider a mathematical model based on the Navier-Stokes equations defined in a moving domain, wherein we immerse the mitral valve leaflets according to a resistive method [30,31]. In particular, we adopt the resistive immersed implicit surface (RIIS) model [32] that enables dealing with a moving immersed surface in an Eulerian framework without explicitly building a surface-conforming mesh. In this framework, we perform a parametric study to assess the sensitivity of clinically meaningful indicators with respect to a parameter representative of the SAM-induced obstruction severity.
This numerical study has the potential to significantly impact on the clinical investigation and treatment of SAM. In particular, we provide preliminary indications about the best region to be removed during myectomy, supporting the surgeon's decision during this procedure.
A chartflow of the proposed computational pipeline, encompassing image-processing, simulation and clinical indications, is displayed in Fig. 1.
The outline of the paper is as follows. In Section 2 we detail the pipeline of image processing that, starting from MRI images, allows reconstructing the geometry and motion of the ventricle (Section 2.1) and of the mitral valve (Section 2.2), and to include the aorta in the model (Section 2.3). In Section 3 we recall the resistive method for the mathematical description. Finally, in Section 4, we show several results obtained in a patient-specific setting, analyzing some useful hemodynamic quantities for different parametrizations of the SAM degree. Conclusions follow.

From image processing to geometric and functional data
In this section we illustrate the methods proposed to reconstruct the geometry and the movement of both the left ventricle endocardium (Section 2.1) and the mitral valve (Section 2.2), starting from medical images. Then, we delineate the strategy adopted to introduce the ascending aorta in the computational domain, by resorting to a template geometry (Section 2.3).
As input to our reconstruction procedure, we consider cardiac Magnetic Resonance Imaging (MRI). In the last years, cardiac MRI, including cine-MRI, has become the gold standard imaging technique to assess myocardial anatomy and global heart function [33]. This non-invasive technique has been widely used both in the diagnosis of hypertrophic cardiomyopathy [34][35][36] and in the investigation of heart valve diseases [37,38]. Therefore MRI is the reference technique to diagnose SAM.
In Section 4.1, the procedure described in the present section will be applied to cine-MRI data of a patient provided by L. Sacco Hospital in Milan. These data come from standard clinical acquisitions, without the need of ad hoc imaging procedures for the purposes of the present work.

Reconstruction of the geometry and motion of the left ventricle endocardium
Cine-MRI has become a standard cardiac image data acquisition procedure, and it features a very accurate time resolution (at least 20 images per heartbeat). Concerning its space resolution, the so-called short axis view is typically the sole 3D representation available, whereas long-axis images are acquired only for few two-dimensional (2D) slices. More precisely, a long-axis image consists in a single 2D slice where the left ventricle is visible in all its extension from the apex to the base (see Fig. 2-b); on the contrary, short-axis views are made of a set of 2D images acquired orthogonally to the long-axis images at various equally spaced positions (see Fig. 2-a).
Short-axis cine-MRI are suitable for the left ventricle anatomy reconstruction for each available time of the heartbeat. However, due to the low resolution along the long axis, they are not sufficient to I. Fumagalli et al. (b) the endocardium geometry is extracted from the segmentations, and its motion is obtained by registration algorithms; (c) from long-axis MRI views, the position of the SAM-affected mitral valve is reconstructed; (d) a 3D leaflet model is obtained combining the data of point (c) with a template geometry; (e) in the geometrical settings reconstructed in the previous steps, CFD simulations provide information on the blood flow and the pressure distribution in the ventricle and ascending aorta; (f) hinging upon the simulations outputs, we can provide indications on the best location where to perform septal myectomy.
accurately reconstruct the movement of the left ventricle. Indeed, the shortening/stretching of the left ventricle along the long axis represents a relevant component of the ventricular contraction/relaxation. Moving from these considerations, we propose a pipeline to reconstruct both the geometry and the motion of the left ventricle endocardium, by combining information between short-axis and long-axis cine-MRI views.
We point out that such a pipeline is designed to be applicable to standard cine-MRI data that are routinely acquired in the clinical procedure. More advanced MRI acquisitions, that may include additional 3D series, are not considered because they are usually available only at the research level, but not in during the standard clinical practice.
The proposed pipeline, which is sketched in Fig. 2, is comprised of the following steps: 1. for each acquisition instant, segmenting the 3D short-axis image (see Fig. 2-a) to produce a 3D surface of the left ventricle endocardium; 2. measuring the apex-to-base distance as a function of time in the long-axis image (see Fig. 2-b); 3. cutting at the base the 3D reconstructions of the endocardium obtained at step 1, coherently with the long-axis measurements obtained at step 2 (see Fig. 2-c); 4. modeling the cut surfaces as artificial high resolution level-set images ( Fig. 2-d); 5. introducing a reference configuration by fixing a time instant (e.g. the end of systole) and apply a registration algorithm among the artificial level-set images obtained at step 4 (see Fig. 2-e); 6. applying the transformation obtained by the registration to the 3D endocardium at the reference configuration (see Fig. 2-e).
Concerning the segmentation (step 1), different automatic or semiautomatic methods for the left ventricle reconstruction have been proposed in the literature [39,40], from the more classical methods based on thresholding and active contours [41] to the more recent techniques based on convolutional neural networks [42,43]. Our proposed pipeline is independent of the segmentation method adopted, the unique requirement being that the reconstructed endocardial surface accurately describes the pathological left ventricle of the patient. For this purpose, we adopt the semi-manual segmentation algorithm proposed in [44] and implemented in the Medical Image Toolkit (MITK) (www.mitk.org) open source software [45,46]. In particular, this algorithm is based on a manual identification of the endocardial 2D contours for a subset of the short-axis slices. Then, it automatically creates a 3D smooth surface exploiting a 3D radial basis function interpolation based on the intensity of the image and on the 2D contours. In Fig. 2-a we show the results of this method applied to our end-diastolic short-axis data.
In order to better reconstruct the motion of the ventricle, we aim at enriching our reconstruction with information about the ventricular shortening and stretching. These components of the motion can be measured from the long-axis images, in particular from the so-called three chamber view, where the left-ventricle is visible together with the left atrium and the aortic root. More in details, we measure the apex-tobase distance at each instants of this view (step 2, see Fig. 2-b) and we use these measures to select the portion of the segmented endocardium accordingly (step 3, see Fig. 2-c).
Then, we aim at registering all the reconstructed images at different instants. Roughly speaking, registration means recovering the geometric transformations or local displacements able to transform the configurations at each instant into the reference one. Registration algorithms are usually applied to images and their accuracy depends also on the resolution of the input image. However, we are not interested in the registration of the whole short-axis image with all the visible human organs (e.g. lungs, bones, . . . ), but only in the left ventricle endocardium. For this purpose, we model the endocardial surfaces of Fig. 2-c as artificial level-set images (step 4) as shown in Fig. 2-d. This choice allows creating images with the same resolution along each axis direction by extending the standard short-axis resolution to the long-axis direction.
Thus, the registration algorithm is applied to the artificial level-set images (step 5). As for the segmentation, our pipeline is independent of the particular registration algorithm chosen. Moreover, since our input are the artificial level-set images which are of binary type, these I. Fumagalli et al.

Fig. 2.
Pipeline to reconstruct the geometry and displacement of the left ventricle endocardium from cine-MRI: (a) short-axis 3D segmentation and corresponding endocardial surface reconstruction; (b) long-axis apex-to-base distance; (c) 3D surfaces obtained by merging the short-and long-axis information; (d) artificial level-set images created from the 3D surfaces; (e) final displacement field obtained by registering the artificial images with respect to the end-systolic instant.
can be easily registered using standard non-rigid registration algorithms [47][48][49]. In particular, here we adopt the non-affine B-splines registration algorithm implemented in the Elastix open source library (http://elastix.isi.uu.nl, [50]), which consists in the minimization of a cost functional depending on the difference between two images and on some regularization terms, under the assumption that the images share the same bounding box. This assumption prevents the registration problem from being ill-posed.
For each time instant, the output of the registration is a 3D vector field defined in the whole 3D domain of the artificial image. This field can be evaluated at all the points of the reference endocardial surface (step 6). In this way, we are able to reconstruct, in the reference configuration of our left ventricle domain, a vector field describing the displacement of the endocardium in the whole heartbeat. In the following, by denoting with MRI the time resolution of cine-MRI, we refer to such displacement field as MRI ( ), being = MRI , = 0, … , = ∕ MRI , the + 1 acquisition times of the cine-MRI and the heartbeat duration. An example of this field, representing the displacement between the end-systolic configuration and the end-diastolic one, is shown in Fig. 2 We remark that, despite the relatively fine temporal resolution of cine-MRI data if compared to other imaging techniques (e.g. CTscans), the image sampling does not match the resolution needed for an accurate time discretization of the PDEs since the time step length of the latter is at least 10 times smaller than the cine-MRI resolution MRI . For this reason, we need to interpolate the displacement field MRI ( ) on a finer time grid, see Section 3.

Reconstruction of the geometry and motion of the mitral valve
As mentioned in the previous section, short-axis views of cine-MRI are well suited for the reconstruction of a chamber like the left ventricle. However, their resolution in the ventricular apex-to-base direction is much coarser (typically 8 mm) than along the other directions. On the other hand, typical long-axis views are time-dependent 2D images made of a single slice, orthogonal to the short axis of the ventricle, and they have a 1-mm resolution in the apex-to-base direction. Since the physiological thickness of valve leaflets is about 2 mm, we base our reconstruction of the mitral valve on long-axis views. In particular, we use the three-chambers view (see Fig. 3-a), that allows appreciating the distance between the anterior leaflet and I. Fumagalli et al. the septum and thus the degree of stenosis induced by SAM. Since imaging data of a three-chambers view are made of a single 2D slice, the sole information concerning the mitral valve that can be extracted is a single line for each time of the MRI acquisition. Thus, in order to have a complete 3D leaflet reconstruction, we complement this data by exploiting a template geometry (see also Section 2.3). For this purpose, we use the Zygote solid 3D heart model [51], a complete heart geometry reconstructed from CT-scans representing an average healthy heart, including the mitral valve leaflets, see Fig. 3-b, top. Moreover, since SAM involves the displacement of the anterior mitral leaflet, we focus our reconstruction on such leaflet. Indeed, during systole the posterior leaflet affects the blood flow only locally, in the posterior basal portion of the ventricle; therefore it should only have a secondary effect on the hemodynamics in the LVOT, where the SAM obstruction takes place.
The pipeline to merge the information of the template geometry into the MRI data aims at obtaining the end-diastolic configuration of the anterior mitral leaflet, including the available patient-specific data, and goes through the following steps: 1. a surface templ describing the anterior mitral leaflet is extracted from the Zygote template mitral valve (see Fig. 3-b, top); 2. the leaflet longitudinal centerline segm is segmented from the end-diastolic three-chambers view shown in Fig. 3-a, left; 3. the centerline obtained at step 2 is extruded in the normal direction to the three-chambers view plane, to obtain a surface extr ; 4. the template leaflet templ is projected onto the extruded surface extr to obtain the surface representing the patient's MV anterior leaflet. The final result is shown in Fig. 3-b, bottom, together with the segmented leaflet centerline segm .
It is worth to notice that the three-chambers view is the main plane on which the severity of SAM is assessed in the clinical practice. To this aim, during its acquisition this plane is set to be orthogonal to the short-axis view and to the septum, thus cutting in half the ventricle and passing through the mitral leaflets centerline. The segmentation step 2 has been performed by means of MITK, while for the other steps, we developed ad hoc semi-automatic tools, based on the Visualization Toolkit (VTK, vtk.org) and the Vascular Modeling Toolkit (VMTK, vmtk. org, [52]). At the end of this procedure, the reconstructed surface = (0) represents the configuration of the anterior mitral leaflet at the end of the diastolic phase, that is at = 0.
In order to include the effects of the mitral valve in a timedependent CFD simulation, the motion of the leaflet needs to be reconstructed. As our target is SAM, we focus on the systolic phase, that spans from the end-diastolic time = 0 to the end-systolic time = S , thus accounting for of the frames available during the complete heartbeat.
Observing the three-chambers view time series, it is possible to notice that during systole, the leaflet progressively approaches the septum, without fluttering nor deforming significantly. Accordingly, the time evolution of the leaflet is then defined as a rigid rotation of its end-diastolic configuration (0) around the normal direction of the three-chambers view. Thus, we introduce the angles ( ), = 0, … , , between the leaflet longitudinal centerline at the current configuration and the one at 0 , as shown in Fig. 3-c, and the associated displacement fields ,MRI ( , ) that map (0) to ( ). Aiming at representing the actual patient condition, the end-systolic angle ( = ) is reconstructed by segmenting the leaflet centerline also at = , and then measuring the angle it makes with the end-diastolic configuration segm described above.
In view of our sensitivity study, starting from the actual value ( ) characterizing the patient at hand, we also consider other two virtual values of ( ), representing a milder and a stronger level of SAM severity, respectively. This will allow us to compare the results in terms of clinical outputs among different SAM scenarios.

Extension of the computational domain using a template geometry
In order to fully capture the hemodynamics effects of SAM, and to prevent boundary conditions from affecting our numerical results, we need to consider an extended computational domain that comprises the left ventricle and the ascending aorta. Unfortunately, standard cardiac cine-MRI are not sufficient for a 3D reconstruction of the ascending aorta, that is only partially visible in long-axis images (e.g. in the three-chambers view, see Fig. 3-a). For this reason, we propose here a strategy to properly merge, at the level of the LVOT, the patientspecific ventricular geometry and motion with a template ascending aorta obtained by the Zygote geometry [51]. Such a strategy consists of the following steps: 1. we start from two input surfaces: the reference end-diastolic patient-specific left ventricle endocardium -denoted hereafter as LV ps -and the upper part of the Zygote left heart -denoted as LH up zyg . The former is the output of the segmentation cut at the mitral annulus and at the outflow tract. The latter is created by joining the internal surface of the Zygote left atrium with the ascending aorta; 2. we perform a rigid registration between the LV ps and LH up zyg ventricular rings using the Iterative Closest Point (ICP) algorithm [53] (see Fig. 4-a); 3. we merge the two surfaces by deforming LH up zyg so that it conforms to LV ps . The applied deformation is obtained as the solution of a Laplace-Beltrami problem [54,55] on a portion of LH up zyg near the ventricular ring. In particular, the vectorial distance between the two rings is used as a Dirichlet boundary condition (see Fig. 4-b); 4. since we are interested in the simulation of the systolic flow, we exclude the left atrium from the domain capping the mitral annulus; 5. eventually, we harmonically extend the displacement field MRI ( ) at the LVOT reconstructed from cine-MRI (see Section 2.1) in the ascending aorta and in the mitral annulus cap, using the same technique of step 3 (see Fig. 4-c).
All the described steps are performed by exploiting the VMTK library [52] and the additional software tools for the cardiac meshing recently proposed in [56], to which we refer the interested reader for further details.

Mathematical and numerical modeling
We describe here the mathematical and numerical framework for the blood flow problem in a moving domain in presence of the mitral valve, which is reconstructed as described in Section 2. Let ( ) be the moving computational domain made of the left ventricle and the ascending aorta; its boundary is partitioned into the following subsets, displayed in Fig. 5: the ventricular endocardium wall , ( ), the aortic wall , ( ), the mitral valve orifice ( ), and the outflow section ( ) located at the end of the aortic root. For the sake of convenience, we set ( ) = , ( ) ∪ , ( ).
The anterior leaflet of the mitral valve is described by a surface ( ) immersed in ( ). Notice that changes in time only due to the effect of the (prescribed) ventricular wall motion, which induces the movement of the aortic root (see Section 2). The contribution of the immersed mitral valve to blood dynamics is accounted for by means of a resistive method [31,32], as explained below. The motion of the domain and of the mitral leaflet are described by the displacement fieldŝM RI and̂, MRI referred to the end-diastolic configuration, reconstructed in Section 2. Since these fields are available only at the MRI acquisition times , = 0, … , S , time interpolation is employed to define such displacements for all ∈ [0, S ) of the systolic phase. Indeed, being SAM a pathology that pertains to systole, we focus our study only on this phase of the heartbeat.
We consider blood as an incompressible Newtonian fluid with uniform density and viscosity . The ventricle contraction can be taken into account by solving the fluid problem at each time > 0 by means of an Arbitrary Lagrangian-Eulerian (ALE) approach (cf., e.g., [57,58]). The presence of the mitral valve is accounted for by the Resistive Immersed Implicit Surface (RIIS) method, introduced in [32]; see also [30,31] for other resistive methods. This method consists in the introduction of an additional term in the momentum equation, which penalizes the kinematic condition, i.e. the adherence of the blood to the moving leaflet. We point out that in the present work the RIIS method is applied for the first time in the ALE context.
In the framework, the continuous problem reads: where , are the blood velocity and pressure, = (∇ + ∇ ) − the corresponding Cauchy stress tensor, ALE is the time derivative in the ALE framework, and ALE the fluid domain velocity. The latter is defined in ( ) as ALE = (̂) •̂− 1 , namely as the time derivative of the fluid domain displacement̂, which is the solution to the following harmonic extension problem, for each ∈ (0, S ]: The velocity fields MRI , m are defined in a similar way, resorting directly to the reconstructed field̂M RI . The resistive term penalizing the leaflet velocity, with a resistance coefficient , has support in a narrow layer around represented by a smoothed Dirac delta type function where is a signed distance function that implicitly describes the immersed surface as = { ∶ ( ) = 0}, and > 0 is a suitable parameter representing half of the thickness of the leaflet. The prescribed leaflet velocity may be provided in different ways, with computational costs progressively diminishing along the following order: (i) as the solution to an additional structural problem for the leaflet; (ii) by a reconstruction procedure based on clinical data; (iii) adopting a quasistatic approach, that is setting ≡ . In the present work, we deem the third option as a suitable approximation: indeed, an estimate of the SAM-induced leaflet tip velocity based on clinical observations is on the order of 0.03 m∕s, which is much lower than 1 m∕s, the characteristic blood velocity in the LVOT.
Regarding boundary conditions for (1), we point out that the scalar o ( ) is a prescribed stress at the outflow section, standing for the aortic pressure, whereas the mitral valve orifice m is kept close and follows the wall motion, since only the systolic phase is considered. Notice also that the aortic orifice is always open, thus justifying our choice to neglect the aortic valve leaflets. Indeed, our attention is on the intraventricular pressure drop and on the stenosis induced by the SAM of the mitral valve. (1) penalizes the kinematic condition = on . As a matter of fact, the RIIS method is one of the possible strategies to account for an immersed surface; some of its features are shared with other methods available in the literature. It has many points in common with the immersed boundary method [59][60][61][62][63][64], but its reduced computational cost is comparable with that of the Resistive Immersed Surface (RIS) method [30,31]. The main advantage w.r.t. the latter is that the computational mesh and the surface do not need to be conforming.

Remark 2.
In the context of cardiac hemodynamics with valves, several approaches are proposed in the literature to account for moving surfaces inside a fluid domain. The variety of the methods can be roughly grouped into two categories: those which entail updating the fluid computational mesh through time, near the moving surface, and those in which the fluid mesh is independent of the immersed surface. The first category comprises the ALE approach [65][66][67][68], cutFEM/XFEM methods [69][70][71][72] amd chimera overset grid methods [73,74]. The second category is mainly made of immersed-boundary/fictitious-domain methods [17,[75][76][77][78][79][80][81][82][83][84][85][86][87]. In general, this second approach -shared also by the RIIS method -entails a lower computational cost, and it can be more easily applied to challenging patient-specific geometries. For a more thorough review about the modeling of heart valves, we refer the reader to [28,[88][89][90][91][92].
In order to numerically solve problem (1), a stabilized piecewise linear finite element approximation is adopted, with a semi-implicit is triangulated into the computational mesh  0 ℎ , ℎ being the maximum diameter of the mesh elements, and the piecewise linear finite element spaces 0 ℎ , 0 ℎ are introduced for the velocity and pressure variables ℎ and ℎ . Then, at each time step, a displacement field ℎ , = 1, 2, … is introduced in terms of the reconstructed field̂M RI and employed to update the computational mesh and the space, before solving the fluid problem. We remind that the geometrical datâM RI and then compute the numerical ALE velocitŷ 3. Move the domain and the mesh according to the ALE map  +1 ∶ 0 → R 3 ,  +1 ( ) = +̂+ 1 ℎ : and update the spaces ℎ , ℎ to +1 ℎ , +1 ℎ . 4. Update the position +1 of the immersed mitral leaflet, according to a time interpolation of the displacement̂, MRI ; 5. Solve a time step of the SUPG-stabilized, linear discrete fluid problem, obtained by a semi-implicit first order approach: where (⋅, ⋅) denotes the scalar product in . The scheme is stabilized by the SUPG/PSPG method [93]. The stabilization parameters M and C are chosen according to the variational multiscale approach [94,95]. In particular, they depend on the local mesh size ℎ and on the local velocity magnitude. This numerical method was adapted for the RIIS model in a rigid domain in [32] and it is extended here to the ALE case.

Numerical results
The reconstruction procedures and the numerical method described in the previous sections are applied to the cine-MRI of a female patient provided by L. Sacco Hospital in Milan. The study has been approved by the Ethics Committee according to institutional ethics guidelines. The patient gave her signed consent for the publication of data. This patient suffered from SAM of the mitral valve and Hypertrophic Obstructive Cardiomyopathy (HOCM), particularly significant in the interventricular septum. No MV regurgitation has been reported. The dimensions of the LV blood chamber (end-diastolic volume of 100.8 ml) and its global systolic function (ejection fraction of 63%) lay within the limits of non-pathological conditions. Moreover, an increased mass has been measured (indexed end-diastolic wall mass of 59.9 g) and the thickness and the systolic thickening of the myocardium are typical of hypertrophic cardiomyopathy. This condition, combined with the pathological motion of the MV anterior leaflet during systole, determines an intraventricular stenosis, whose effects on hemodynamic quantities are investigated by using our numerical software.
The results presented in the following sections are threefold: 1. in Section 4.1, we report the outcome of the reconstruction pipeline described in Section 2, that provides the endocardium displacement MRI and the configurations of the mitral valve for different levels of SAM severity; 2. in Section 4.2, the effect of SAM severity on hemodynamic quantities is assessed by means of CFD results; 3. in Section 4.3, we focus on some specific quantities in order to provide useful indications for the design of septal myectomy, a possible surgical treatment for the patient's condition.
In Fig. 6-a we highlight some significant regions of the computational domain, exploited in Sections 4.2 and 4.3 to compute relevant fluid dynamics quantities. The computational results of these sections are obtained by means of the finite element library LifeV [96], www. lifev.org on the tetrahedral mesh shown in Fig. 6-b,c for the initial time = 0 and the end of systole = S . An average mesh size of ℎ = 1.2 mm is chosen for most of the fluid domain, while a refinement to ℎ = 0.6 mm   Table 1. Concerning the resistive method with resistance , the characteristic length is chosen to ensure that the corresponding smeared Dirac delta , spans over two mesh elements [32]. Regarding boundary conditions, since we are interested in the systolic phase, the mitral orifice is closed by a no-flow condition ( m = ALE on m in (1)), whereas a time-dependent pressure o is applied as Neumann condition to the outflow section of the ascending aorta (see (1)). Since no pathological alterations of the systemic pressure is reported for the patient at hand, we consider a physiological behavior of the aortic pressure o ( ), reconstructed from Wiggers diagrams [97,98] and represented in Fig. 6-d. We point out that o , which is the baseline value upon which the pressure gradients develop in the LVOT, attains its maximum at = 0.3 s: in the following, we will refer to this time instant as the systolic peak P = 0.3 s.

Application of the reconstruction pipeline
Employing the segmentation and registration pipelines described in Sections 2.1 and 2.2, the kinematics of the left ventricle and of the mitral valve have been reconstructed from cine-MRI. In this section, we assess the results of this reconstruction, which represents the input to the CFD simulations presented in Section 4.2.
The reconstructed displacement field MRI of the endocardium surface is shown in Fig. 7. Observe that this motion well captures the macroscopic contraction and shortening of the ventricular chamber, as well as the displacement variations characterizing the hypertrophy of the cardiac muscle.
In order to assess the reliability of the LV reconstruction in terms of the available clinical measurements, the time evolution of the LV cavity volume is reported in Fig. 8. The shape of the volume curve well represents the behavior of the LV as known from Wiggers diagram [97,98]. Moreover, comparing the values of the reconstructed End Diastolic Volume (EDV) and Ejection Fraction (EF) with those measured during the clinical acquisition, we notice a substantially good agreement, see Table 2. In this respect, it is worth to point out that the measure of the volume provided in clinical tests is obtained from the length of the two main LV axes and having regarded the LV as an ellipsoid. This represents a coarse approximation of the actual LV geometry, and the discrepancies w.r.t. our reconstructed measures can be ascribed to it.
In order to study the effect of SAM on blood flow quantities in the LVOT (see   Fig. 7. Displacement field MRI of the LV endocardium with respect to its endsystolic configuration: on the left, initial systolic configuration ( 0 = 0 s); on the right, mid-diastolic configuration ( 12 = 0.6 s). In Fig. 9, the reconstructed MV at end diastole and end systole are shown for the three scenarios P-, PS, P+. We observe that our method is able to provide significant MV configurations which could then be successfully used as geometric immersed data in our resistive CFD-ALE model.
In summary, our reconstructions allow us to obtain at each discrete time the computational domain and the moving ventricular wall , together with its velocity MRI to be prescribed as Dirichlet condition, see (1). Regarding the moving aortic wall , and the corresponding velocity, these are obtained by extrapolation of , and MRI . Completing the blood pool geometry, our reconstructions allowed us to identify the mitral valve orifice . Notice that all these data are independent of the parametric end-diastolic angle ( ) used in three scenarios P-, PS, P+. In addition, for each of such scenarios, we reconstruct the immersed MV leaflet at each time .

Hemodynamics of systolic anterior motion
In this section we analyze the hemodynamics in the left ventricle and in the ascending aorta in the presence of SAM. In particular, we focus on the systolic phase and investigate the sensitivity of fluid dynamics quantities on the SAM severity parameter ( S ) introduced in Section 4.1 (cf. Fig. 9). The aim is to obtain possible indications about the development of the pathology and the appearance of secondary effects.
Starting from the geometries and motions reconstructed and reported in Section 4.1, the ALE displacement and velocitŷℎ,̂A LE,ℎ , = 0, 1, …, are introduced. We recall that during the systolic phase ∈ [0, S ] the mitral orifice is closed. Moreover, a physiological pressure o is applied at the outflow section of the ascending aorta, with a peak value at time P = 0.3 s (see Fig. 6-d).
In the following, the three scenarios P-, PS, P+ will be compared, in order to study how the level of severity of SAM affects the blood flow. We point out that the same LV geometry and displacement are considered in the three scenarios, so that the comparison identifies the effects of SAM alone, separating them from those due to the HOCM that concurs in the patient at hand. We point out that, although the flow rate exiting the ventricle should decrease for increasing LVOT narrowings, with an effect on the LV dynamics, the same LV geometry and displacement are considered in the three scenarios. This allows us to isolate the effects of SAM alone, separating them from those due to the HOCM that concurs in the patient at hand.
In Fig. 10 the blood velocity field is represented on a 2D longitudinal slice to capture the SAM-induced hemodynamics near the mitral valve. We can observe that, in all the cases, the LVOT restriction due to the combination of septum hypertrophy and SAM yields a jet-like flow through the aortic orifice. Due to this jet structure, the characteristics of the flow in the rest of the ventricle can be considered having almost negligible influence on the flow patterns, apart from providing a flow rate that is consistent with systolic contraction.
Comparing the flows displayed in Fig. 10 for the different scenarios, one can notice that the intensity of the aforementioned jet grows with the level of SAM severity. This is shown also in Table 3, where increasing values of the peak-systolic maximum value of the velocity magnitude max ∈ | ( P , )| are attained for increasing end-systolic leaflet angles ( S ). In this respect, we notice that the high values in the velocity observed for the P+ case are in accordance with the strong LVOT restriction that the displacement of the mitral leaflet determines.
In order to assess the reliability of our results, we compared them with echo-Doppler evaluations made on the same patient, which report a peak velocity of 5.5 m/s. The value of 6.18 m/s obtained by our simulations in the actual patient-specific scenario PS (see Table 3)  shows quite a good agreement with the clinically measured value. We point out that echo velocity measurements are taken on a specific direction, which is not necessarily the one along which the maximum of | | is attained; thus, the echo-Doppler could slightly underestimate the peak velocity.
About the pressure distribution in the ventricle and proximal aorta, Fig. 11 shows how SAM induces large pressure jumps in the LVOT: in comparison to these variations, the pressure distribution can be considered as uniform in the other regions of the computational domain. In particular, larger values of the pressure gradient are observed where the leaflet is nearer to the septum, thus yielding a stronger obstruction to the blood flow. Moreover, significant differences are noticeable among the pressure jumps of the three scenarios P-, PS, P+, both in the plots of Fig. 11 and in the peak-systolic values of Table 3, last row. This confirms the role of SAM in determining an intraventricular pressure I. Fumagalli et al.

Table 3
Peak-systolic values of flow quantities for the three SAM-severity scenarios considered (P−, PS, P+). Pressure values are averaged over the LV-inner region (see Fig. 6 jump, that is much higher than the subaortic stenosis due to the septum hypertrophy. Comparing the simulation results with echo-Doppler measurements, we observe a good agreement between the peak pressure jump of 106 mmHg for the PS scenario (see Table 3) and an echo-based clinical estimation of 100 mmHg.
One can also compare our results on the pressure jump with the literature, e.g. with the peak value data measured by catheterization on HCM patients (with and without SAM, indistinctively) reported in [7]. Such data lays in a range from 0 to 150 mmHg, thus our actual patient-specific scenario PS is in line with these measures, in particular with a severe HCM condition: see Fig. 11, where our values for PS lay in the range 25-106 mmHg. This is consistent with the fact that SAM can strongly worsen the hypertrophy-induced LVOT obstruction. Accordingly, the P-scenario, in which the SAM severity is artificially damped, shows pressure gradients that are typical of relatively mild HCM. On the other hand, it is worth to point out that the pressure values obtained in the P+ case are much larger than in a physiological or even pathological condition. Indeed, this case Table 4 Peak-systolic values of the enstrophy density  represents an extreme situation, in which the very strong severity of SAM induces an obstruction that would hardly be sustained by the cardiac muscle. Nevertheless, including in our comparison the results obtained in such a condition allows to quantify the clinical relevance of a correct assessment of this pathology, which can degenerate if not promptly identified or evaluated.
In order to provide a more thorough description of the fluid dynamics with SAM, the vorticity of the flow and the Wall Shear Stress (WSS) are also inspected.
Concerning vorticity, Fig. 12 displays the vortical structures according to the Q-criterion [99]: contour lines of the quantity = 1 2 are represented, where and are the skewsymmetric and the symmetric part of the velocity gradient, respectively. We can see that coherent structures develop in the aortic root, as a consequence of the jet flow entering from the aortic orifice, and that the intensity of these structures is strongly related to the SAM severity degree. In order to quantify such intensity, which is related to kinetic energy dissipation and possibly to transition to turbulence, we compute I. Fumagalli et al.   Table 4.The high intensity of in the LVOT is related to the strong jet flow, whose energy is then progressively dissipated while giving rise to the aforementioned structures in the aortic root, and then flowing through the ascending aorta. Comparing the values reported in the table for the different cases P-, PS, P+, we notice that vorticity is more pronounced if a more severe SAM condition is considered.
Regarding WSS, it is known that elevated values could be associated to a weakening of the vessel wall, thus creating conditions that may favor aneurysm formation [100]. Therefore, it is interesting to consider its distribution, displayed in Fig. 13 at the systolic peak = P . We notice that very high values are reached on the septal wall of the LVOT, as a consequence of the high-velocity flow in this region. As expected, the intensity of WSS increases with the severity of the SAM condition. Moreover, it is possible to identify also a region of the ascending aorta characterized by high WSS values. In this location, the jet flow coming from the LVOT impacts on the aortic wall and this impact is stronger in the case of a more severe SAM condition. This suggests that this region should be periodically monitored by clinical tests, in order to assess the risk of damage of the wall structure and thus prevent possible weakening and aneurysm formation [100].

Towards clinical indications
One of the most employed treatment for patients suffering from HOCM and SAM is septal myectomy, that is an open-heart intervention in which a portion of the interventricular septum is surgically removed. Under precise circumstances [3], this procedure is typically characterized by low perioperative mortality and high long-term survival rate [101,102]. Notice that, in critical conditions mitral valvuloplasty can be performed together with myectomy.
In this section, we want to provide preliminary quantitative indications that can help the surgeon in the preoperative design of the myectomy, in particular in determining the region and the size of the myocardium portion to be resected. In common practice, indeed, this decision is taken on the basis of diagnostic images and the surgeon's experience. Our results could provide more quantitative indications in this respect.
We focus on the difference = − o between the pressure field and its aortic value, because this is the quantity by which clinicians I. Fumagalli et al. typically assess the severity of the obstruction. More specifically, a line running along the septum is identified and parametrized by its curvilinear coordinate from the aorta to the apex (cf. Fig. 14-b), and the distribution of on such line is displayed in Fig. 14-a. One can notice that the narrowing of the LVOT due to the hypertrophy of the septum yields a pressure drop in a region near to the aortic orifice (green region in Fig. 14-a). This is consistent with the Venturi effect that is often claimed to be the cause of SAM: the reduced pressure in the outflow tract of the ventricle drives the anterior mitral leaflet towards the septum. In the same figure, we can also observe that intraventricular pressure increases with , with strong gradients in the subaortic region (especially in the grey region in Fig. 14-a) and reaching its maximum near the leaflet tip.
As a clinical indication, the present discussion identifies two regions of the LVOT (highlighted in Fig. 14-a,b) in which the surgical intervention may have positive outcomes: intervening in the green, lowpressure region (between = 28 mm and = 32 mm in our case) can possibly dampen the Venturi effect causing SAM; on the other hand, removing a portion of the septum in the grey region (between = 32 mm and = 44 mm in our case) can reduce the large pressure gradients and the intraventricular stenosis. These sites are those commonly targeted in surgical practice for septal myectomy [3,103]. We remark that the patient-specific indications on the myectomy region that we provide clearly depend on the anatomy of the patient and on his/her level of SAM severity. Nevertheless, the proposed computational procedure does not make strong assumptions on the geometry of the patient's LV, therefore it can be applied with no modification to standard cine-MRI data of any patient.
For completeness, in Fig. 14-c we also provide the pressure profiles for the other two cases P-, P+ considered in our parametric study. The differences among the three degrees of SAM severity are significantly large, therefore we can state that SAM has a major effect on the pressure levels inside the left ventricle. Moreover, this observation highlights the importance of an accurate reconstruction of the patient specific mitral valve configuration, in order to properly assess the severity of the condition. In this respect, we remark that the discussion presented here compares the pressure gradients that different SAM severity levels induce, assuming that the ventricle undergoes the same displacement and thus that the cardiac muscle is strong enough to always provide the same flow rate through the aorta. Although this assumption may lead to an overestimation of the resulting pressure gradients, it allows to isolate the effects of SAM alone on the hemodynamics, separating them from the effects of HOCM and of additional dysfunctions in the ventricle contraction.
We deem that also the analysis of the WSS could give indications about the location where to perform septal myectomy. For this reason, in Fig. 14-d we report its values on a line along the septum. In this representation, we can appreciate how the maximum value of WSS, attained at the narrowing, strongly depends on SAM severity, thus confirming the importance of a proper evaluation of this parameter in the assessment of SAM.

Conclusions and limitations
In the present work, a CFD study on the effects of SAM in presence of HOCM has been carried out. A novel reconstruction pipeline has been proposed to obtain the ventricle and mitral valve geometry and motion from kinetic imaging data. Hinging on this data, we could avoid an FSI approach that would entail a complex modeling of the mechanics of the myocardium and the composite valve structure in pathological conditions. The blood flow in the prescribed moving domain has been modeled with Navier-Stokes equations in ALE form, enhanced with a resistive method to represent the immersed mitral valve. The main outcomes of the present work can been summarized as follows: • The proposed reconstruction pipeline has been effectively applied to standard cine-MRI data, which are routinely acquired in diagnostic procedures; • Through a parametric study, the dependence of hemodynamic quantities on the level of SAM severity has been assessed. These results showed that the subaortic stenosis due to the combination of HOCM and SAM yields a jet-like flow through the LVOT, with sensible pressure jumps and localized high WSS in the ascending aorta. A comparison with echo-Doppler measurements showed a good agreement between the patient-specific simulation results and the measured data and clinical estimations, in terms of peak blood velocity and pressure gradients; • Driven by a commonly employed surgical treatment of HOCM and SAM, namely septal myectomy, special attention has been paid to the distribution of pressure along the septum, in order to provide practical indications for the preoperative design phase.
The present study has the following limitations, that open the perspective to further investigation: • The reconstruction pipeline of the left ventricle geometry and motion is mainly based on the short-axis series of cine-MRI data, with some information coming from the 3-chambers view. A more comprehensive approach encompassing many different acquisition series is currently being developed by the authors. • The mitral leaflets motion is described as a rigid rotation along a fixed axis. Although this motion shows a quite fair consistency with imaging data in the systolic phase that we consider, it does not capture the 3D features of valve dynamics. Considering nonrigid, fully 3D displacements may allow to provide indications also on the latero-lateral extension of myectomy. Along this direction, including also the posterior leaflet in the reconstruction would be advisable, due to its influence on the flow near the lateral walls of the ventricle. • We considered an anatomically non-pathological mitral valve, even though SAM might occur also in valves that show abnormalities. The computational analysis here proposed could be applied with little modifications also to patients with such abnormalities: this would allow to include also their effects.
The results of the present work give preliminary clinical indications that can supplement clinical data in the assessment of the condition and design of the treatment. Directions of further investigation may consider a wider set of patients, with different with different anatomical and functional conditions, including the inspection of the hemodynamic effects of different valve geometries. Then, it would be possible to devise a global score for SAM severity, that can help in classifying the patients and deciding whether and how to surgically intervene on them.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.