Joel Kronborg

Frida Svelander

Samuel Eriksson-Lidbrink

Ludvig Lindstr Öm

Carme Homs-Pons

Didier Lucor

Johan Hoffman

Computational analysis of flow structures in turbulent ventricular blood flow associated with mitral valve intervention

Keywords: patient-specific heart modelling, left ventricle, mitral valve clip, finite element method, turbulent blood flow, triple decomposition of velocity gradient tensor

Cardiac disease and clinical intervention may both lead to an increased risk for thrombosis events due to a modified blood flow in the heart, and thereby a change in the mechanical stimuli of blood cells passing through the chambers of the heart. Specifically, the degree of platelet activation is influenced by the level and type of mechanical stresses in the blood flow. In this article we analyze the blood flow in the left ventricle of the heart through a computational model constructed from patient-specific data. The blood flow in the ventricle is modeled by the Navier-Stokes equations, and the flow through the mitral valve by a parameterized model which represents the projected opening of the valve. A finite element method is used to solve the equations, from which a simulation of the velocity and pressure of the blood flow is constructed.

The intraventricular blood flow is complex, in particular in diastole when the inflow jet from the atrium breaks down into turbulent flow on a range of scales. A triple decomposition of the velocity gradient tensor is then used to distinguish between rigid body rotational flow, irrotational straining flow, and shear flow. The triple decomposition enables the separation of three fundamentally different flow structures, that each generates a distinct type of mechanical stimulus on the blood cells in the flow. We compare the results to a simulation where a mitral valve clip intervention is modelled, which leads to a significant modification of the ventricular flow. Further, we perform a sensitivity study of the results with respect to the positioning of the clip. It was found that the shear in the simulation cases treated with clips increased more compared to the untreated case than the rotation and strain did. A decrease in valve opening area of 64 % in one of the cases led to a 90 % increase in rotation and strain, but a 150 % increase in shear. The computational analysis suggests a process for patient-specific simulation of clinical interventions in the heart with a detailed analysis of the resulting blood flow, which could support clinical risk assessment with respect to platelet activation and thrombosis events.

INTRODUCTION

Heart attack and stroke follow from a loss of blood supply to the heart and brain, respectively, often caused by a blood clot that obstructs a critical blood vessel, referred to as thrombosis. Blood clots are naturally formed as a response to injury to avoid blood loss, but may also form under other circumstances where the platelets (thrombocytes) of the blood are activated which leads to an increased tendency of the platelets to aggregate, see e.g., Previtali et al. (2011). Platelet activation can be triggered by contact with a foreign material, such as a medical device, see e.g., Rojano et al. (2019), or by the mechanical stresses in the blood flow. Specifically, it is associated with elevated levels of shear in the blood. The degree and time extent of the activation is a function of the type of mechanical stimulus experienced by the platelets, see e.g., Rahman et al. (2020). Therefore, a detailed analysis of the structure of the blood flow may help to assess the risk for platelet activation and thrombosis. The heart is supplied blood through the right and left coronary arteries, with inlets through the right and left coronary sinuses which are part of the aortic root attached to the left ventricle (LV) of the heart. Further downstream the aorta, the brain is supplied blood through the carotid arteries that emanate from the aortic arch. As a consequence, the level of platelet activation in the LV blood flow is of interest when assessing the risk for thrombosis events in both the heart and the brain.

Blood flow in the LV is complex, in particular in diastole when the heart relaxes and blood flows into the ventricle from the left atrium through the open mitral valve. The ventricular blood flow is also highly sensitive to the structure and function of the mitral valve, and pathologies such as mitral valve regurgitation or stenosis can drastically change the flow pattern. In the early filling phase of diastole (the E-wave) a jet is formed which then breaks down in the diastasis phase, after which a second jet is formed in the atrial contraction phase (the A-wave). The two successive jets create shear layers and vortex rings in the ventricle, which then break down into turbulent flow structures on a range of scales. The blood flow generates mechanical stimuli on blood cells that trigger different responses in the cells, for example, activation of platelets. Mechanical stimuli on blood cells are difficult or impossible to assess by medical imaging alone, and the turbulent nature of the ventricular flow makes this challenge even greater. In this article we investigate to what degree a computational model can be used to assess the mechanical stimuli on blood cells in the LV flow, building on our previous work on patient specific simulations of the LV blood flow, see e.g., [START_REF] Larsson | Patient-specific left ventricular flow simulations from transthoracic echocardiography: robustness evaluation and validation against ultrasound doppler and magnetic resonance imaging[END_REF] and [START_REF] Larsson | Estimation of left ventricular blood flow parameters: clinical application of patient-specific cfd simulations from 4d echocardiography[END_REF], and analysis of turbulent flow structures, see [START_REF] Hoffman | Energy stability analysis of turbulent incompressible flow based on the triple decomposition of the velocity gradient tensor[END_REF].

Cardiac disease and clinical interventions may both lead to changes in the LV blood flow and the mechanical stimuli experienced by the blood cells in the ventricle. Mitral valve regurgitation is a condition which affects the LV blood flow, where the mitral valve does not close properly and therefore may leak in systole when the valve is supposed to be closed. A common procedure to mitigate mitral valve regurgitation is to suture the two leaflets of the mitral valve together with a a mitral valve clip (MVC). While a MVC intervention reduces or removes the regurgitation, there is a risk for creating a new problem in the form of stenosis, due to a reduced valve opening. Recently, patient-specific computer simulations have been used to evaluate the risks and benefits of a MVC intervention, see [START_REF] Caballero | A comprehensive engineering analysis of left heart dynamics after mitraclip in a functional mitral regurgitation patient[END_REF]. In this article we use a finite element method to simulate the blood flow in the LV before and after a MVC intervention, to analyze the difference in the ventricular blood flow. We also perform a sensitivity study with respect to the positioning of the MVC. A two dimensional model of the mitral valve is used, which represents an approximation of the projected opening between the left atrium and the left ventricle. The mitral valve model is then modified to take the positioning of the MVC into account.

To analyze the simulated blood flow we use the triple decomposition of the velocity gradient tensor, suggested by [START_REF] Kolář | Vortex identification: New requirements and limitations[END_REF] to improve the identification and visualization of vortices. [START_REF] Hoffman | Energy stability analysis of turbulent incompressible flow based on the triple decomposition of the velocity gradient tensor[END_REF] noted that the triple decomposition implies that the flow locally can be approximated by a sum of shear flow, rigid body rotational flow and irrotational straining flow. Our hypothesis is that this refined analysis of the blood flow will lead to an improved quantification of the type and extent of the mechanical stresses experienced by blood cells while inside the LV, which can support a risk assessment of thrombosis events to a greater extent than a blood flow analysis based solely on the velocity and pressure of the blood.

In Section 2 we describe the computational model and the triple decomposition of the velocity gradient tensor, and in Section 3 we present results from the simulations. These results are then discussed in Section 4.

METHOD

Left ventricle model

The LV model is based on 4D transthoracic echocardiography (TTE) images of the left ventricle of one human subject, acquired using the technique described in detail in [START_REF] Larsson | Patient-specific left ventricular flow simulations from transthoracic echocardiography: robustness evaluation and validation against ultrasound doppler and magnetic resonance imaging[END_REF]. Images were captured in 26 frames over a full cardiac cycle which lasted for 830 ms, and the heart wall was segmented in a semi-automated fashion in each image, see [START_REF] Hansegård | Semi-automated quantification of left ventricular volumes and ejection fraction by real-time three-dimensional echocardiography[END_REF]. From this a set of triangulated surface meshed were generated which described the deformation of the endocardium over the cardiac cycle.

The orifice of the mitral valve (MV) was identified by a sonographer in thirteen 2D short axis B-mode views in early diastole, with each view rotated 15 degrees around the LV long axis. Similarly, the aortic valve opening was identified in a single 3-chamber B-mode view in early systole. Both valve openings were marked on the surface meshes.

From one surface mesh corresponding to end diastole, here referred to as the initial surface mesh, a refined version was produced, from which a tetrahedral volume mesh of the interior of the ventricle was generated.

This volume mesh consists of 406,964 vertices. This resolution was chosen based on the convergence studies in [START_REF] Larsson | Patient-specific left ventricular flow simulations from transthoracic echocardiography: robustness evaluation and validation against ultrasound doppler and magnetic resonance imaging[END_REF], which showed less than 2 % deviation in results when increasing mesh resolution above 400,000 vertices.

Simulation of blood flow

The blood in the LV is modelled as an incompressible Newtonian fluid. Although blood is known to express non-Newtonian properties in smaller blood vessels, see e.g., Pedley andFung (1980), and[START_REF] Adjoua | Blood flow reduced-order modeling across macroscopic through mesoscopic scales[END_REF][START_REF] Adjoua | Blood flow reduced-order modeling across macroscopic through mesoscopic scales[END_REF], the blood inside the heart chambers is typically modelled as a Newtonian fluid, see Westerhof et al. (2010). To simulate the LV blood flow, we determine the velocity vector u(x, t) : Ω t → R 3 and scalar pressure p(x, t) : Ω t → R, such that the Navier-Stokes equations are satisfied:

ρ( u + (u • ∇)u) -µ∆u + ∇p = 0, (x, t) ∈ Ω t × [0, T ]
(1)

∇ • u = 0, (x, t) ∈ Ω t × [0, T ], (2)
defined over the time-dependent domain Ω t ⊂ R 3 at time t ∈ [0, T ]. Equation (1) expresses conservation of momentum, and equation (2) conservation of mass. Standard mathematical notation is used, with ∇ the gradient operator and ∆ the Laplacian operator. We use the dynamic viscosity µ = 0.0027 Pa•s and the blood density ρ = 1060 kg/m 3 , see Di [START_REF] Martino | Fluidstructure interaction within realistic three-dimensional models of the aneurysmatic aorta as a guidance to assess the risk of rupture of the aneurysm[END_REF].

Ω t represents the deforming LV, with the boundary ∂Ω t partitioned into the heart wall and the two valve orifices. We use no slip velocity boundary conditions, corresponding to setting the blood flow velocity equal to the deformation velocity of the heart wall and closed valves. In systole, the open aortic valve (AV) is modelled by a standard outflow boundary condition, where the pressure is set to zero. The model of the mitral valve is described below, with and without a MVC intervention. To simulate the LV blood flow over the cardiac cycle, the domain Ω t is approximated by the tetrahedral volume mesh which is deformed over time based on the 26 surface meshes. For a detailed description of this process, see [START_REF] Larsson | Patient-specific left ventricular flow simulations from transthoracic echocardiography: robustness evaluation and validation against ultrasound doppler and magnetic resonance imaging[END_REF].

The Navier-Stokes equations are solved on this deforming volume mesh using an aribtrary Lagrangian-Eulerian (ALE) formulation of a stabilized finite element method, which is described in Hiromi Spühler and Hoffman (2021).

Mitral valve model

We extend the LV model with a simple planar, data-driven parameterized model of the MV dynamics, which acts as a time dependent inflow boundary condition. The MV is modeled as a time-varying orifice representing the area enclosed by the planar projection of the leaflets onto the mitral annulus, accounting for virtual folding/unfolding mitral leaflets, similar to the work in [START_REF] Imanparast | The impact of valve simplifications on left ventricular hemodynamics in a three dimensional simulation based on in vivo mri data[END_REF].

The basic shape of the annulus is modeled as two half ellipses sharing a common long axis. From the ultrasound identification of the MV, the length of the long axis and the two short axes are identified, as well as the spatial 2D location of the coaptation line within the annulus area. An MV center point on the wall of the ventricle mesh is also given from the TTE images, which is used as a spatial reference point on the mesh to keep the MV in place. In figure 1a the MV and AV areas on the simulation mesh are displayed, as well as the contour of the MV model and the projected leaflets for a particular moment in time in figure 1b.

The opening and closing of the MV are given as temporal variations of the short axes of the two half ellipses. During systole, the short axes are equal, effectively keeping the projection of the leaflet edges overlapped along the coaptation line. Throughout diastole, the short axes are scaled to imitate a projection of the moving leaflets onto the surface mesh. Effectively this corresponds to a maximum opening area at E-wave peak, slightly smaller in diastasis, and somewhere in between at A-wave peak. For these points in time in the diastole timeline, the short axis scalings are determined by an idealized MV, see [START_REF] Fyrenius | Major and minor axes of the normal mitral annulus[END_REF]. Patient-specific scalings are difficult to determine from echocardiography images since the resolution is limited. The predetermined scale factors are given in table 1 for E-wave peak, diastasis start and end, and A-wave peak. Although the patient-specific scaling of the short axes could not be determined from the echocardiography images, the timings for the E-wave and A-wave peaks, as well as diastasis start and end, are determined based on the expansion rate of the LV of the particular subject studied here. This can be done patient specifically for other cases as well.

The volumes of the 26 frames extracted from the echocardiography images are displayed in figure 2, plotted against time from the beginning of systole. A cubic smoothing spline curve, V (t) is fitted to the data points. From the first and second derivative (also included in figure 2) of this smooth curve the timings are identified. Fundamentally, the shift from systole to diastole is identified as the point where the first derivative dV /dt (i.e. the rate of the volume change) is 0, shifting from negative (contraction) to positive (expansion). The shift from diastole to systole is given simply as the beginning/end time of this heartbeat.

Note that the isovolumetric phases between diastole and systole are excluded for simplicity.

The E-wave and A-wave peaks are easily identified as the local maxima in dV /dt, with E-wave peak being the former and higher one. To define diastasis, the local minimum between these two peaks is determined. Diastasis is then defined as the time between the two original data frames from the TTE images on either side of this local minimum. These timings are listed in table 1.

With the timings of these four crucial points determined, along with the timings for shifts from diastole to systole and vice versa, the scaling of the short axes of the MV model are given by cubic Hermite interpolation between the scaled axis lengths and times given in table 1. The resulting axes scaling over the full duration of diastole is displayed in figure 3. 1 for exact values). The bottom plot shows the second derivative.

After mitral valve clip intervention

To simulate the inflow of blood after a MVC intervention, part of the coaptation line of the MV was kept stationary in its closed state throughout diastole. This produces a valve model with two smaller openings, one on either side of the clip as displayed in figure 1c. This is achieved by creating two smaller openings consisting of two half ellipses each, similar to the untreated case, and superpositioning this with the coaptation line. The two smaller openings behave much like smaller versions of the untreated MV model throughout diastole. They follow the same pattern for opening and closing, i.e. the same scalings and timings given in table 1. In its closed state, the clipped MV model has the same coaptation line as the untreated valve model.

Since the length of the coaptation line is the same for the clipped and untreated models, the long axes of the two smaller openings are directly given from the remaining unclipped parts. Additionally, the relationship between the anterior and posterior short axis of each opening is the same as the relationship of 1.

the short axes in the untreated case. The length of the two anterior short axes are assumed to have the same relationship as the length of their corresponding long axes, and the same holds for the posterior short axes.

To compute the maximum length of the short axes in the MVC case, the length of the edge of the openings at E-wave peak is assumed to be the same as the length of the edges in the untreated case, i.e. the circumference of the orifice. Along with the relationships between different axes described above, this is used to uniquely determine the maximum length of the short axes. As in the untreated case, this is assumed to imitate a projection of the 3D MV opening onto the surface mesh of the ventricle, see [START_REF] Imanparast | The impact of valve simplifications on left ventricular hemodynamics in a three dimensional simulation based on in vivo mri data[END_REF]. Given the length of the ellipse axes and the MVC, the total area of the two openings were computed and are presented in table 2.

Simulations

Five different cases were simulated, with MVC centers positioned along the coaptation line at s • l from the coaptation line center, with l being the same long axis length as used in the untreated valve, i.e. half the intercommissural distance, and s ∈ (-1, 1) a scaling parameter. The values of s used in these simulations are given in table 2, along with the resulting MV opening areas at mid diastole. Hence, s = 0 corresponds to the clip being in the center of the valve and s = ±1 to it being on the posteromedial/anterolateral commissure of the valve respectively.

Table 2. Position scale constants for the MVC simulations and resulting MV opening areas at E-wave peak. The MVC is positioned at s • l from the center of the coaptation line, with l being the length of the long axis in the MV model. Here s < 0 corresponds to clip positions anterolaterally of the center, and s > 0 posteromedially. For reference, the area of the MV opening at mid diastole with no MVC is 2.54 cm 2 . s -0.4 -0.2 0.0 0.2 0.4 Area (cm 2 ) 1.08 0.93 0.86 0.93 1.08

In addition to the five cases with the MVC in different positions, one case with an untreated valve was also simulated using the original MV model. Note that the untreated valve was modeled without any regurgitation. This is a model simplification motivated by the fact that our focus is to compare the blood flow into the ventricle during diastole. In this model the inflow is not greatly affected by any regurgitation in systole. The significant difference is that between the flow in the different double orifice cases in MVC treated valves and the single orifice of the untreated valve.

The LV used in these simulations has a volume at end diastole of 117 ml, and an ejection fraction of 64 %.

These values, as well as the deformation of the ventricle volume throughout the heartbeats, are the same in all simulations. The untreated MV has a long axis of 12.8 mm, and short axes of 11.6 mm (posterior) and 8.1 mm (anterior) respectively. This results in an opening area of 4.0 cm 2 at E-wave peak, and 2.54 cm 2 at mid diastole. This is a relatively small opening area, so the resulting velocities and pressures in the MVC simulations are likely higher than a typical in-vivo case.

All simulations were run on the KTH supercomputer of Technology (2021), a Cray XC40 system with 67,456 cores distributed over 4,120 Intel Xeon E5-269X 2.1/2.3 GHz CPUs. The simulations were run for six heartbeats, starting from early systole, all with identical wall motion from the captured TTE images.

The time step length was set to 1 ms. One of the 26 surface meshes generated from the TTE data diverged significantly from the others, and was thus considered an outlier and excluded from the simulations.

Between the rest of the captured frames, the surface meshes were temporally interpolated using cubic Hermite interpolation, to yield one surface mesh for each simulated time step.

Blood flow analysis

To analyze the blood flow in the ventricle we can visualize the flow velocity and pressure computed from the Navier-Stokes equations. It is also useful to analyze the velocity gradient ∇u, a tensor which describes the local spatial change in velocity at each time and position, defined by its components

(∇u) ij = ∂u i ∂x j . (3)
Near each point x 0 ∈ Ω t , the velocity gradient can be used to construct a linear approximation

u(x) ≈ u(x 0 ) + ∇u(x 0 )(x -x 0 ), (4)
and the mechanical stresses in the blood can be derived from the velocity gradient.

A standard double decomposition of the velocity gradient constitutes a separation of the flow into a straining flow and a rotational flow, corresponding to its symmetric part S(u), the strain rate tensor, and the skew-symmetric part Ω(u), the spin tensor:

∇u = 1 2 (∇u + (∇u) T ) + 1 2 (∇u -(∇u) T ) = S(u) + Ω(u), (5)
where the superscript T represents the transpose of the tensor. The weakness of the double decomposition is that shear flow is unaccounted for. To address this shortcoming, the triple decomposition was proposed by [START_REF] Kolář | Vortex identification: New requirements and limitations[END_REF], where the velocity gradient tensor is decomposed into a sum of three parts corresponding to irrotational straining flow (compression and elongation) (∇u) EL , rigid body rotational flow (∇u) RR , and shear flow (∇u) SH :

∇u = Q((∇u) EL + (∇u) RR + (∇u) SH )Q T , ( 6
)
where Q is an orthogonal matrix that represents a new frame of reference in which the shear flow can be subtracted from the other two parts. For the sake of brevity, we will refer to the three flow structures as strain, rotation and shear.

Algebraically, the matrix decomposition may be interpreted as a sum of a normal symmetric part, a normal skew-symmetric part, and a non-normal part, which can be computed through the real Schur decomposition of the velocity gradient tensor, see e.g., Van Loan and Golub (1996). In this article, the triple decomposition was computed by using the Matlab function schur, which returns a frame of reference Q that corresponds to a standardized form of the real Schur decomposition from which the shear flow tensor can be extracted, see [START_REF] Hoffman | Energy stability analysis of turbulent incompressible flow based on the triple decomposition of the velocity gradient tensor[END_REF]. To determine the magnitude of each flow component, Frobenius norms of each respective matrix were computed.

RESULTS

Flow structures

Snapshots of the velocity magnitudes in a bicommissural vertical cut through the ventricle for all simulated cases are shown in figure 4. The snapshots are taken at time t = 0.45 s from the start of the analyzed heartbeat (cf. figure 3), i.e. 30 ms after the start of diastole. At this point in time inflow jets through the mitral valve can be observed in all cases. The untreated valve has a single wide jet, whereas the cases with MVC's in different positions show the double jets characteristic of inflow after MVC treatment.

The average velocity in the jet of the untreated case (figure 4a), which has a larger MV opening area than the others, is smaller than the velocities in the clipped cases. This single wide jet also does not reach as far down into the ventricle as the jets in the other cases. Comparing the treated cases, the one with the clip in the center of the coaptation line (figure 4d) has a slightly higher peak velocity than the others. This case is the one with the smallest total opening area, as given in table 2. In the ventricles with asymmetrically placed clips, the wider of the two jets typically reaches further into the ventricle than the smaller jet.

Snapshots of the pressure at the same instant and the same cut as the velocity described above are shown in figure 5. A recurring pattern of high pressure in the apical and mid-cavity parts of the ventricle is apparent in all simulated cases, with slight variations. In the untreated case (figure 5a), where the jet does not extend as far into the ventricle, the overall pressure variations are lower than in the cases with MVC's.

In the MVC cases, the double jets extend to different depths in the ventricle, which can also be seen in the 2.

velocity in figure 4. This makes the line between high and low pressure less straight in the MVC cases, especially when the MVC is asymmetrically positioned. High pressures just below the jets are apparent, but the lowest pressures are found around the front of the jets. This indicates the presence of vortex rings in these parts, since pressure is typically low in vortex centers.

Triple decomposition of velocity gradient tensor

Placing a clip in the MV significantly reduces the opening area of the MV, as can be seen in table 3. A natural result of this is higher inflow velocities, which the table also displays, with higher velocities the smaller the opening is. Both the velocity and the pressure are mainly affected by the reduction in area, but not so much by which side the clip is placed on. This is however not necessarily the case for the triple decomposition components, as is evident from the same table.

With the clip placed in the center of the MV, the opening area is the smallest out of all these simulation cases. For shear and strain, this leads to the highest computed peak values, but the peak rotation magnitude is actually slightly larger with the clip placed slightly posteromedially of the center. It is true for all three 2.

components that the peak values with the clip displaced to either side are not only dependent on how far the clip is displaced, but the direction as well. However, in most cases it does seem like smaller opening area does contribute to higher magnitudes of rotation, shear and strain, although the effect varies. The MV with the clip in the center has an opening area of 34 % of that of the untreated MV. Its peak average rotation and strain are both roughly 190 % of those for the untreated case, but its peak average shear is about 250 % of that for the untreated case. As a reference, the maximum velocity increased about twofold in the same case.

Throughout all simulation cases, shear has the highest relative increase when treated with MVC, whereas rotation and strain increase roughly equally.

There are significant differences in how the flow in the ventricle behaves when dominated by strong inflow jets, i.e. during E-wave, and when there is no single dominating structure dictating most of the flow.

Figure 6 shows snapshots of the three modalities of the triple decomposition, i.e. rotation, shear and strain, in the simulation case without the MVC. The components are displayed at two different times in diastole: t = 0.45 s, which corresponds to the velocities shown in figure 4, and t = 0.71 s, which is the time of the 6c) appears in small connected areas in the lower part of the ventricle, but notably also on either side of the lower part of the jet structure, showing the presence of a vortex ring forming around the jet.

Within the jet itself, where the flow is roughly unidirectional, the rotation is 0.

Shear (figure 6a) at the same moment is high at the outer edge of the jet, and along the lower wall of the ventricle. There are weaker shear structures throughout much of the rest of the ventricle as well, but none within the jet. Strain (figure 6e) is also large along the wall around the apical part of the ventricle.

Additionally, below the jet, at the front of the blood flowing into the ventricle from the atrium, straining flow is high. There are also smaller spots throughout the ventricle with relatively high or low strain.

At t = 71, when the E-wave jet has diminished and there is no single strong feature dominating the flow field, the overall distribution of the different flow modalities is more mixed. Rotation (figure 6d) is loosely connected in structures throughout the ventricle. Shear (figure 6b) and strain (figure 6f) are still relatively large along the ventricle wall around the apex, but the shear now also extends to other parts of the wall. In some structures near the wall but extending into the interior of the ventricle the strain is relatively high.

When it comes to the cases treated with MVC's, some of the appearing flow structures share similarities with those in the untreated case. Snapshots of rotation, shear and strain for the MVC case with position scaling s = 0.2 (corresponding to case c in figures 4 and 5) are displayed in figure 7. The snapshots are taken at the same moments in time as the ones in figure 6, and show the same cut through the ventricle.

Note however the different scales for shear and strain.

With the flow profile dominated by two jets, as in early diastole in figure 7a, c and e, the rotation (c) again shows structures reminiscent of vortex rings around the lower part of the jets. (For the corresponding velocity snapshot, see figure 4c.) However, the magnitude of the vortex rings appear larger for the parts in the interior of the ventricle than where the jet runs close to the wall.

Along the edges of the jets, high shear is visible in the clipped case as well (figure 7a). The shear of the wider jet seems to extend into the center of the jet as well, but this is because the displayed cut is not perfectly through the center of the jet. Like in the untreated case, there is also significant shear along the wall in the lower part of the ventricle. The strain in the clipped case (figure 7e) is once again large at the bottom of the jet where the inflow meets slower moving blood. It is also high at the wall in the apical part of the ventricle, and at the edges of the MV openings. There is also an area of higher strain right in the center of the ventricle, between the two jets. At the time of the second snapshot, i.e. the end of diastasis when the strongest effects of the initial E-wave the apical part of the ventricle where the rotation is relatively high (figure 7d). Much like in the untreated case, the shear (figure 7b) extends along more parts of the ventricle wall in the second snapshot, with weaker effects in the interior part. The strain (figure 7f) retains its high magnitude along the apical part of the ventricular wall. Differences between the triple decomposition components in the treated and untreated cases include, perhaps most notably, significantly higher magnitudes in strain and shear in the clipped case. The magnitude scales here have been adjusted to better show the location and structures of the different modalities, but figure 8 indicates that the claim holds. In the snapshots of the early diastolic flow, interactions between the two jets in the clipped case give rise to effects not seen in the untreated case, such as the high strain in the center of the ventricle. But vortex rings, high shear along the edge of the jets, and high strain at the bottom of the jets appear in both the treated and untreated cases. Shear and strain along the walls are also prevalent in both cases.

To get an idea of the duration and intensity of the different flow modalities from the triple decomposition over a full heartbeat, figure 8 shows the magnitude of the rotation, shear and strain for all simulated cases over time. There are some similarities apparent between all modalities. In systole, they all vary much less than in diastole, with the exception of the untreated case. Across all modalities, the magnitudes in the untreated case peak at much lower values than the MVC cases, despite being at similar levels in systole. It is difficult to distinguish any clear differences between the simulation cases in systole.

In diastole, all modalities follow similar temporal patterns in all MVC cases. With the increasing influx of blood in early diastole, the magnitudes of the rotation, shearand strain rise. Shortly after the E-wave peak, the magnitudes also peak and then start to fall. They continue decreasing towards systolic levels, with no real noticeable effects from the A-wave. However, despite having the same inflow volume of blood as the MitraClip cases, the magnitudes of rotation, shear and strain do not rise high above the systolic levels in the untreated case. In fact, after a small initial rise right at the start of diastole, they drop during some of the time leading up to the E-wave peak.

DISCUSSION

The triple decomposition used here gives an excellent view of the locality and structure of different flow modalities in the LV blood flow. Structures such as strong inflow jets can be visualized and quantified by studying the velocity, and with the addition of pressure, phenomena like the vortex rings around jets can be visualized. However, studying e.g. the shear around the jets requires more sophisticated tools. The triple decomposition captures this, as well as confirms expectations from the velocity and pressure plots.

The pressure in the untreated case is relatively low in the basal part of the ventricle, but slightly higher further down. This pattern repeats in the MVC cases, and some of these have even lower pressure where the vortex ring around the inflow jet is expected. The rotation parts (a and b) of figures 6 and 7 confirm that the rotation magnitude is indeed large in those parts. The high strain in the apical part of the ventricle also corresponds well to the high pressure in the same area. A related interesting feature of the rotation and strain is that areas with high rotation often tend to correspond to low strain, and vice versa. An exception to this observation is within the jets, where rotation is 0 and strain is also relatively low. See e.g. figure 6c and e.

To get a better idea of some of the 3D structures that appear in the LV flow, and how they develop over time, figure 9 displays snapshots of some 3D structures defined by the computed rotation. Both the untreated and MVC case display clear vortex rings forming around the fronts of the jets as connected structures with rotation magnitude greater than a threshold value of 80 s -1 (figure 9a andb). Above these the shear structures are visible at the border around the jet, with significantly higher magnitude in the MVC case. In the rest of the ventricle at this point in time both simulation cases display relatively little rotation, with only small connected volumes reaching above the threshold.

At t = 0.71 s, when the E-wave jets have decomposed, figure 9c andd gives a very qualitative view of the distribution of high rotation and shear structures throughout the ventricle. In the untreated case (c), the rotation appears as connected features, but with very little in the apical part of the ventricle. The shear is most prominent along the wall. In the clipped case (d), the rotation structures appear more numerous, and throughout the entire ventricle in different sizes. The shear is again concentrated on the wall, with higher values as indicated by the color bar. Interestingly, figure 8 at this point in time indicates that the total amount of rotation/shear magnitudes are similar in both the untreated and clipped (s = 0.2) case. We have not investigated this further, but the reason could be that the more intense patterns visible in the MVC case make up a relatively small part of the entire ventricle volume, and that the integrated values in figure 8 are thus more strongly influenced by large regions of smaller values.

From figure 9 the differences between the distributions of high rotation and shear structures at different points in time is evident. During the E-wave, much is centered around the jets. Later in diastole, there is a wider distribution over the ventricle. Since high shear is known to contribute to thrombosis formation, being able to quantify this in the interventricular blood flow is desirable. Our results suggest that the locations of high shear structures in diastole are similar before and after MVC treatment. They appear mainly around the inflow jets and along the ventricle walls. However, after MVC treatment the magnitude of the shear is higher, possibly as a result of the higher inflow velocities caused by the smaller MV opening areas. As seen in figure 8, the integrated magnitude over the ventricle is higher for all MVC cases than for the untreated case throughout most of the E-wave and diastasis. This implies that the risk of platelets experiencing higher shear which could lead to thrombosis formation is increased after MVC treatment, but further studies would be required to quantify this risk.

While shear has been clearly linked to thrombosis formation, the case is more open for strain. Our results indicate that the locations of strain structures, like shear, do not change much in the LV after MVC treatment. But as for shear, strain is higher after treatment, both in maximum magnitude and integrated over the ventricle during the E-wave and diastasis.

Since shear, strain and rotation affect platelets differently, it is useful to separate these modalities using the triple decomposition, and be able to study them separately. Often vorticity has been used to similar ends.

However, vorticity fails to separate rotation from shear, which can lead to misleading conclusions. Rotation only transports the blood cells, and thus likely does not contribute to platelet activation. Our results clearly show that both shear and rotation are present in the simulated blood flow, but in very different locations and of different magnitudes. Rotation is roughly one order of magnitude lower than shear, and is more prevalent in the vortex ring around the jet, and in connected vortex structures throughout the interior of the ventricle after the strong E-wave jet has diminished.

To make a visual comparison between the vorticity and the triple decomposition, figure 10 the large rotation structure in the center of the ventricle. Thus using vorticity as a measure of mechanical stimuli on blood cells may lead to inaccurate conclusions. In this particular case the magnitude of the shear is large in comparison to the rotation, which might limit the errors if vorticity were to be used in studies as a substitute for shear. In other settings this may not be the case, or one might be more interested in studying the rotation without influence of the shear. Either way, the triple decomposition offers a very promising way to distinguish rotation and shear in a way that vorticity fails to do.

The simulated velocity in the untreated case (figure 4) is within a reasonable range of a healthy heart, but the peak velocities in the clipped cases may be indicative of mitral stenosis [START_REF] Otto | acc/aha guideline for the management of patients with valvular heart disease: A report of the[END_REF]. The opening areas of the clipped cases (table 2) are also quite small, and MVC treatment would likely not be clinically considered to mitigate regurgitation if the resulting areas were so small. Nevertheless, increases in inflow velocity as well as the triple decomposition components in the LV following MVC treatment are still likely to appear in clinically realistic cases due to decreased MV area. While our results show that MVC treatment increases the magnitudes of rotation, shear and strain within the LV in early diastole, more information is needed to make any claims about the risk of platelet activation or thrombosis as a consequence of MVC treatment. Not only the magnitude of the shear and strain are needed to assess such risks, but the duration for which platelets experience these higher levels are of interest as well. However, studying the magnitude and locations of different flow component structures are an important first step, which opens up for new opportunities to learn more about mechanical stimuli in blood flow.

Figure 1 .

 1 Figure 1. (a) AV placement (yellow) and MV placement (red) on the LV surface mesh in an axial view. (b) The two half ellipses making up the MV model at different stages throughout the heartbeat cycle. The solid red line is the projection of the anterior leaflet, and the solid blue line the posterior. The dashed lines are the short axes of the corresponding half ellipse. In the top view the valve is fully closed (systole), and the leaflets overlap on the coaptation line. In the bottom view the valve is fully open (E-wave peak), and the center view shows early diastole. The solid black line through the valve model is the long axes of the two half ellipses. (c) The coaptation line (left), the two double half ellipses making up the MV model after MVC intervention (center) and the final shape of the MV model after MVC intervention (right). The part that does not open is the location of the clip.

Figure 2 .

 2 Figure 2. Volume of the acquired surface meshes of the LV, and the first and second derivative of the interpolated curve. The top plot shows the volumes of the 25 used frames directly from the TEE images as blue dots, and the interpolation as a red line. The middle plot shows the first derivative, and the resulting timings of diastole start, E-wave peak, diastasis start and end, A-wave peak, and diastole end as red circles (see table1for exact values). The bottom plot shows the second derivative.

Figure 3 .

 3 Figure 3. Lengths of the short axes of the anterior (blue) and posterior (red) leaflets in the MV model over time. The times marked by red circles are E-wave start (Es), E-wave peak (Et), E-wave end/diastasis start (Ee), A-wave start/diastasis end (As), A-wave peak (At) and A-wave end (Ae). The times and scalings of the short axes are given in table 1.

Figure 4 .

 4 Figure 4. Snapshots of the velocity magnitude at time t = 0.45 s in a vertical bicommissural cut through the left ventricle. This is in early diastole before the E-wave peak, when the intraventricular blood flow is dominated by strong inflow jets. Subfigure (a) shows the untreated case. The other subfigures have MVC's positioned at s • l from the center, with l being the long axis of the MV model. The scale factor s for the different cases are (b) s = -0.2, (c) s = 0.2, (d) s = 0.0, (e) s = -0.4, (f) s = 0.4. The corresponding MV areas for each case in mid diastole are given intable 2.

Figure 5 .

 5 Figure 5. Snapshots of the pressure at time t = 0.45 s in a vertical bicommissural cut through the left ventricle. This is in early diastole before the E-wave peak, when the intraventricular blood flow is dominated by strong inflow jets. The regions just below the jets experience high pressure, while there is a pressure drop around the base of the jet. (a) shows the untreated case. The other subfigures have MVC's positioned at s • l from the center, with l being the long axis of the MV model. The scale factor s for the different cases are (b) s = -0.2, (c) s = 0.2, (d) s = 0.0, (e) s = -0.4, (f) s = 0.4. The corresponding MV areas for each case in mid diastole are given intable 2.

Figure 6 .

 6 Figure 6. Snapshots of shear flow (a and b), rotational flow (c and d) and straining flow (e and f) at times t = 0.45 s (a, c and e) and t = 0.71 s (b, d and f) in a vertical bicommissural cut through the left ventricle. This simulation was done with an untreated valve. At time t = 0.45 s (early diastole), the flow is dominated by the incoming E-wave jet. The velocity and pressure at this time are shown in figures 4a and 5a respectively. At time t = 0.71 s (end diastasis) the inflow jet is weaker, and the initial dominant flow structures of the E-wave have started to break down. Note that the different sizes of the ventricle are due to the diastolic expansion.

Figure 7 .

 7 Figure 7. Snapshots of shear (a and b), rotation (c and d) and strain (e and f) at times t = 0.45 s (a, c and e) and t = 0.71 s (b, d and f) in a vertical bicommissural cut through the left ventricle. This simulation was done with MVC placed at s • l from the center of the coaptation line, with l being the long axis in the MV model, and the position scale factor s = 0.2. At time t = 0.45 s (early diastole), the flow is dominated by the two incoming E-wave jets. The velocity and pressure at this time are shown in figure 4c and 5c respectively. At time t = 0.71 s (end diastasis) the inflow jets are weaker, and the initial dominant flow structures of the E-wave have started to break down. Note that the different sizes of the ventricle are due to the diastolic expansion.

Figure 8 .

 8 Figure8. Magnitudes of rotation, shear and strain averaged over the entire ventricle volume and plotted against time over one heartbeat, starting with systole, for all simulated cases. The legend indicates the location of the MVC as relative shift from the center of the coaptation line, except for the healthy (untreated) case. The gray dashed lines mark the time of the start of diastole, E-wave peak, diastasis start, and A-wave peak, also given in table 1. The black dashed lines mark t = 0.45 s and t = 0.71 s, i.e. the times at which the blood flow is visualized in figures 4, 5, 6, 7, 9 and 10.

Figure 9 .

 9 Figure 9. Snapshots of 3D flow structures in the LV blood flow, as determined by the rotation magnitude. The opaque blue structures are all cells that have a rotation magnitude greater than 80 s -1 . The transparent structures are cells with shear magnitude greater than 500 s -1 , colored in accordance with the color bars. (a) Untreated case, t = 0.45 s. (b) MVC case, s = 0.2, t = 0.45 s. (c) Untreated case, t = 0.71 s. (d) MVC case, s = 0.2, t = 0.71 s. Corresponding 2D snapshots are displayed in figures 6 and 7.

  shows a snapshot of the vorticity for one of the simulated MVC cases, along with the three components from the triple decomposition at the same instant (same as in figure7a, c and e, but with a logarithmic scale). From this figure it is clear that the vorticity shares similarities with the shear. They reach similar magnitudes, and form structures mainly around the inflow jets and along the walls of the ventricle near the apex. But the logarithmic scale allows for observing the rotation contribution to the flow in the vorticity as well. E.g. the vortex rings show quite clearly as connected areas of high rotation magnitude around the bottom of the jets in figure10d. (Note also the low strain in the vortex rings in figure10c.) Since the vortex rings are rotation dominated, they do not show up as clear structures in the shear part (figure10b). They do however clearly influence the vorticity magnitude around the lower parts of the jets in figure10a, especially

Figure 10 .

 10 Figure 10. Snapshots of (a) vorticity, (b) shear, (c) strain and (d) rotation at time t = 0.45 s in a bicommissural vertical cut through the left ventricle. This simulation was done with a MVC placed at s • l from the center of the coaptation line, with l being the long axis in the MV model, and the position scale factor s = 0.2. The velocity and pressure at this time are shown in figure 4c and 5c respectively.

Table 1 .

 1 Timings of different parts of diastole, and scalings of the short axis anterior (SAA) and short axis posterior (SAP) of the MV model at those times. The timings are determined by analyzing changes in the rate of expansion of the LV. Note that the coaptation line does not coincide with the intercommissural line (long axis), but the scalings at diastole start and end correspond to both leaflet edges coinciding with the coaptation line, i.e. the MV is closed at these times.

		Diastole start E-wave peak Diastasis start Diastasis end A-wave peak Diastole end
	t (s)	0.420	0.525	0.644	0.711	0.779	0.830
	SAA scale	0.40	1.00	0.50	0.50	0.85	0.40
	SAP scale	-0.57	1.00	0.70	0.70	0.85	-0.57

Table 3 .

 3 Peak values in the untreated simulation case, and the corresponding values for each of the MVC cases. Area is given at mid diastole. Velocity and pressure are the maximum values reached anywhere in the LV at E-wave peak. Rotation, shear and strain are the maxima reached in figure8. Thus the first snapshots (figure6a, c and e) show the triple decomposition at a point in time when the flow is dominated by the strong inflow through the MV in early diastole, and the second snapshots capture the components after the initial E-wave jet has broken down into smaller flow structures in the ventricle.

	Simulation case	untreated -0.4	-0.2	0.0	0.2	0.4
	Area at mid diastole (cm 2 )	2.54	1.08 0.93 0.86 0.93 1.08
	E-wave peak velocity (m/s)	1.86	2.84 3.53 3.76 3.51 2.84
	E-wave peak pressure (kPa)	0.667	2.396 3.631 4.640 3.621 2.993
	Peak average rotation (s -1 )	21	35.4 37.0 39.8 39.9 34.1
	Peak average shear (s -1 )	174.6	344.8 397.5 439.1 409.7 368.7
	Peak average strain (s -1 )	51.4	82.5 86.7 95.3 93.6 81.3
	end of diastasis.

At t = 45, it is clear that different flow modalities are prevalent in different parts of the ventricle. Rotation (figure

This is a provisional file, not the final typeset article

Frontiers

jets have diminished, the rotation structures are fairly well distributed in the ventricle, with some areas in

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.

AUTHOR CONTRIBUTIONS

FUNDING AND ACKNOWLEDGEMENTS

The authors acknowledge the financial support from the Swedish Research Council (grant number 2018-04854), and the access to computing infrastructure by the Swedish National Infrastructure for Computing (grant number SNIC 2018/3-296).


Bernard Gibaud
email: bernard.gibaud@irisa.fr

THE DICOM STANDARD: A BRIEF OVERVIEW

Keywords: Medical Imaging, PACS, DICOM 1

The DICOM standard has now become the uncontested standard for the exchange and management of biomedical images. Everyone acknowledges its prominent role in the emergence of multi-vendor Picture Archiving and Communication Systems (PACS), and their successful integration with Hospital Information Systems and Radiology Information Systems, thanks to the Integrating the Healthcare Enterprise (IHE) initiative. We introduce here the basic concepts retained for the definition of objects and services in DICOM, with the hope that it will help the reader to find his or her way in the vast DICOM documentation available on the web.

INTRODUCTION

Exchanging images between various kinds of equipment has been an issue since the very beginning of digital medical imaging. The ACR/NEMA standard, issued in 1985 and 1988, did not provide a satisfactory solution and the medical imaging community had to wait until 1993 to have a real usable standard available. The publication of this standard was followed by several demonstrations in key international congresses, especially the RSNA in 1993, which convinced the whole community that the standard actually led to successful implementations, and that the manufacturers were willing to use it, and turn the page of proprietary formats and solutions. Since that time, the position of DICOM was significantly reinforced, with the creation of the DICOM Committee, and the collaboration with other standards development organizations such as the Comité Européen de Normalisation (CEN) in the mid-nineties, and ISO TC 215 "Health Informatics" in 1999.

DICOM has now become the international standard in the field of biomedical imaging. Its influence was critical in the emergence of multivendor technical solutions for Picture Archiving and Communication Systems (PACS), and in providing appropriate solutions for the integration with the other information systems involved, especially the Hospital Information Systems and the Radiology Information Systems. These integration issues were addressed by the Integrating the Healthcare Enterprise Initiative (IHE) -initially in the USA in 1998, under the auspices of the RSNA and HIMSS, then world-wide, with the launching of IHE-Europe in 2001, and similar initiatives in Asia more recently.

THE DICOM COMMITTEE

The DICOM committee gathers various kinds of contributors. All the big manufacturers of the field of biomedical imaging are members of the DICOM Committee. There are also many professional societies involved, in radiology, but also from other fields such as cardiology, dentistry, pathology, ophthalmology, as well as stakeholder institutions such as the American National Cancer Institute.

The scope is quite broad with currently 26 working groups addressing the various aspects involved (Table 1.) Among those working groups, one has to distinguish Working Group 6 ("Base standard"), whose function is to guarantee the overall technical consistency of the standard, and Working Group 10 ("Strategic advisory"), in charge of advising the DICOM Standard Committee on long term orientations that should be followed.

ORGANIZATION OF THE STANDARD

A multi-part document

The DICOM standard is organized as a multi-part document, made of 18 independent parts (Figure 1) [1]. The most important ones are Parts 3, 4, 16, 5 and 2. Part 3 "Information object definitions" provides the specification of the information objects to be exchanged (more than 1000 pages of text), as well as the definition of the semantics of each data element. The main reason why this part is so long and complex is related to the many existing imaging modalities (such as computed tomography, ultrasound, magnetic resonance, positron emission tomography etc.), that require many technical parameters. Part 4 "Service Class specifications" defines the services for exchanging information, either images or information that is useful to manage images. Part 16 "Content Mapping Resource" addresses the question of terminology, i.e. on the one hand, it defines how existing terminological resources can be used in DICOM (e.g. SNOMED, LOINC, UCUM), and on the other hand, how content items can be grouped together and re-used in DICOM Structured Reporting documents (notion of Templates). Part 5 "Data structure and encoding" specifies how the information objects specified in Part 3 can be organized into a linear bit stream, in order to be sent over a network connection or stored in a file. All aspects related to image compression are addressed in Part 5. Finally, Part 2 "Conformance" specifies how a manufacturer can claim conformance to the DICOM standard for a particular product or implementation, by writing a document called a "conformance statement". Part 2 explains in detail how this document must be written and the information it must contain.

Extensions and maintenance

The standard is being continuously updated, either to address new fields of biomedical imaging (for example, surgery applications), or to revisit existing fields (e.g. Magnetic resonance or PET images), in order to face new needs resulting from the evolution of imaging techniques.

Such extensions are made through the publication of "Supplements" to the standard. Several such Supplements are issued every year. Each Supplement includes text modifying the existing parts, e.g. to add a new information object in Part 3, or a new service class in Part 4, or new codes in Part 16, or new data elements in Part 6 "Dictionary" etc. All those changes are put together once a year to deliver a new issue of the standard, e.g. "DICOM 2008". However, this yearly issue should not be seen as a new "version" of the standard. There is only one single version of the standard, since all changes that are made are supposed to be "backward compatible": it means that new specifications bring new functionalities, but they never change the way previous services work. The only exception to this rule is the case of retirement, i.e. when service classes are retired from the standard (this happens only to services that are almost never used). The general motivation for this kind of approach is that imaging equipment can be in operation during quite long periods of time (6, 7, up to 10 years) and that it would be non-acceptable to change the standard in ways that would "break" existing implementations, and cause troubles to installed and working equipment.

A second way of updating the standard exists, through the mechanism of Correction Proposals (CPs). CPs provide the means to raise and fix minor issues and imprecision that cause (or may cause) interoperability problems. CPs are submitted, accepted, worked out and balloted as any other parts of the standard.

4.

BASIC CONCEPTS

Data Syntax

The syntax used in the DICOM standard is a binary one, based on the transfer of (Tag-Length-Value) triplets. It means that any elementary Data Element is given a binary Tag, e.g. Data element (0028,0010) Rows represents the number of rows in an image. This Tag, composed of 4 bytes is sent first, followed by the corresponding data type, represented in this particular case by ASCII characters 'U' and 'S' (for Unsigned Short), followed by the length of the value field, in this case "2" because the Data element (0028,0010) Rows is represented by an Unsigned Short, represented by 2 bytes, and finally comes the value itself, e.g. "256", represented in binary ('00' in the least significant byte, and '01' in the most significant byte). So the binary bit stream corresponding to this Data element would be in this example: … '28' '00' '10' '00' '55' '53' '02' '00' '00' '01' … As can be seen from this example the values are binary ones, which explains that images in DICOM cannot be edited like, e.g. XML documents. One must use a parser that decodes the TLV triplets in order to put the information in human-readable form.

Information objects

The specification of an information object (called in DICOM an Information Object Definition or IOD) consists in defining the set of Data elements that shall or may be transmitted. They are organized in groups called 'modules'. For example, the CT image IOD is depicted by a number of modules focusing on: (1) the general context of image acquisition (essential information on the patient, the study, the series), (2) the acquisition procedures (particularly the physical acquisition methods, the reconstruction algorithm etc.), (3) the image's characteristics (size, resolution etc.), and (4) the pixel data themselves. The 'module' concept, by gathering together data elements relating to one same information entity (for example 'Patient Module', 'General Study Module', 'General Image Module', 'Image Plane Module', 'Image Pixel Module'), facilitates their reuse in different IODs. These information entities are defined using information models, following the "entity-relationship" formalism. The basic DICOM model is a hierarchy in which a Patient can have one or more Studies, each containing one or more Series, each containing one or more Composite Objects (such as Images, Presentation States, etc.)

A distinction is introduced between Composite objects and Normalized objects. Composite objects gather data elements describing several entities of the real world, whereas Normalized objects focus on one single real world entity.

All sorts of images, whatever the modality (Computed Tomography, Magnetic Resonance, Ultrasound, XRay angio, etc.) are represented by Composite objects, as well as Presentation States, Fiducials, Registration objects, Structured Reports etc.

At the beginning of the standard, image objects were essentially 2D images. Multi-frame objects were introduced to represent in a single Dataset several related images, e.g. in ultrasound (Multi-Frame Ultrasound IOD) and nuclear medicine (Nuclear Medicine IOD). However, the new image objects developed recently for CT, MR, PET etc. (Enhanced CT, Enhanced MR, Enhanced PET, respectively) make use of a new extended multi-frame mechanism. This leads to a better and more efficient organization of common attributes (using so-called functional groups), as well as to faster transmission times.

Services

Two sets of services exist. A first set of services (called Composite services), applies to Composite objects, whereas a second set (called Normalized services) concerns Normalized objects. The first provide facilities to 'push' an object (service C-STORE), or to query a set of objects based on some search criteria (services C-FIND, C-GET and C-MOVE). A last service exist (C-ECHO) to simply test connectivity between two peer application entities.

Normalized services provide facilities to act upon a single object managed by a remote application entity, to create an instance (N-CREATE), delete it (N-DELETE), get information about it (N-GET), modify it (N-SET), launch an action (N-ACTION), report an event (N-EVENT-REPORT).

Service-Object Pair

Objects, i.e. IODs, can be associated with services to form what is called a 'Service Object Pair'. For example, a CT image IOD is associated with the Storage service to form the CT Image Storage SOP Class. This abstraction defines a service that can be provided by an application entity (this entity is called the SCP, for Service Class Provider), or used by another application entity (this entity is called the SCU, for Service Class User). These notions of SOP, SCU and SCP are extensively used in DICOM Part 2 to claim the conformance to DICOM of an implementation. A SOP Instance corresponds to the instantiation of the use of a service, e.g. for a particular image. SOP Classes and SOP instances have unique identifiers based upon the OSI Object Identification (numeric form) as defined by the ISO 8824 standard. The SOP instance UIDs provide identifiers for persistent objects such as images.

DICOM MAJOR ACHIEVEMENTS

Storage and query / retrieve of Composites objects

Storage and Query & Retrieve are the services classes used to communicate any sorts of Composite objects. In fact the exchange and query mechanisms are totally generic with respect to the objects to be exchanged. Requests are based on the use of unique identifiers for Studies, Series and Composite objects (using the SOP Class and SOP instance UIDs).

Composites objects may be exchanged in two ways, either in connected mode, or using removable media. In the first case, composite objects are transmitted as one message per Dataset (a Dataset corresponds to a Composite object instance), after a preliminary stage called Association establishment, during which the two peer application entities negotiate the SOP Classes and Transfer Syntaxes they want to use in the exchange. In the second case, the composite objects are copied on removable media, together with a file called DICOMDIR, providing a structured description of the Patient / Study / Series / Composite Objects relationships, with links to corresponding files on the media. DICOM does not specify any particular convention for the naming of files. Many kinds of media are supported like Cederoms, DVD, Flash memory etc.

Image management services

After providing the basic services for exchanging images, DICOM developed services for image management, in close collaboration with CEN TC251/WG4 "Medical Imaging", a standardization body that was quite active at that time in the area of image management [START_REF] Gibaud | Standardisation in the field of medical image management : the contribution of the MIMOSA model[END_REF]. Such services were acutely needed, both in the US and in Europe, for the development of multivendor PACS systems.

Modality Worklist was one of those services whose impact was prominent for successful PACS deployment. Modality Worklist provides the means for a modality to query an information system (e.g. a RIS) about the tasks that have been scheduled for this particular piece of equipment for a certain period of time. More precisely, the modality queries the Modality Scheduled Procedure Steps that were registered in the information system, and this mechanism provides the modality with essential identifying information about the Patient and Study concerned. This is really essential because it is the key for having consistent identifiers in the RIS and the PACS.

A second important service was developed at the same time, it is called "Storage Commitment". The basic idea is for the modality (or any kind of equipment producing images) to get an explicit confirmation that all the images that have been sent to the image manager have been well received and stored in a safe place, so that the modality can delete those images and free the corresponding storage spaces.

A third image management service is "Performed Procedure Step", allowing to "close the loop" with the information system. This service allows the information system to be notified about events such as the beginning and completion of the examination and the creation and updating of Performed Procedure Step objects, documenting what was precisely done during a procedure step (which could eventually differ from what was scheduled).

Additional services were added, such as "General Purpose Worklist", to support the organization of any sort of task (e.g., quality control, image processing, Computer Assisted Detection or CAD, reporting, transcription), based on a local definition of these tasks.

All these management services have been extensively used in the IHE Integration Profiles specified for Radiology applications, especially those addressing the issue of workflow (for example Scheduled Workflow, Post-Processing Workflow, Reporting Workflow).

5.3

Structured reporting DICOM's major focus is on images. However, it became obvious that standards were needed to support the exchange of observations made on images, whatever the origin of such observations, i.e. the analysis of images by a human observer, or through the use of an image analysis algorithm.

DICOM tried (and succeeded !) in addressing such issues in a generic way. The solution was provided by the Structured Reporting concept (SR), introduced in DICOM with Supplement 23, published in 1999.

The basic idea is to represent the observations related to a Structured Reporting document as a tree, composed of content nodes (called content items) connected by relationships [START_REF] Bidgood Wd | Clinical importance of the DICOM structured reporting standard[END_REF]. Such observations may represent measurements made on images, refer to regions of interest in specific images, or document the context of observation. DICOM introduced 14 content item types, listed in Table 2. What is important to note is that each node has a concept name, that is necessarily a coded entity, and a value, which may be coded or not, depending on the kind of node. In terms of new DICOM IODs and SOP Classes, Supplement 23 introduced three new classes, addressing various levels of complexity in the use of available content items. The first, called Basic Text SR IOD, addresses the needs of applications exchanging very simple documents, limited to text referring to images. The second, called Enhanced SR IOD, provides the same features plus the possibility to represent measurements and to refer to spatial and temporal coordinates. Finally, the third, called Comprehensive SR IOD may involve any kind of content nodes and authorizes the use of 'by reference' relationships (ability to refer to a node anywhere in the SR tree).

Another important feature of DICOM SR was the introduction of Templates. A Template is the specification of a sub-tree that can be instantiated in an SR document. It groups and specifies content item nodes in roughly the same way as modules do, but in a way that is more formalized and dedicated to tree-like structures. Basically, a Template is a table of all content item nodes of a sub-tree, specifying things like requirement type (mandatory or optional), multiplicity, code values allowed for the Concept name attached to each content item node, and eventually, code values allowed for the value of content item node (when represented by a coded entry).

Interested readers can refer to the book written by David Clunie [START_REF] Clunie D | DICOM Structured Reporting[END_REF] for (far) more details about DICOM SR.

Since the publication of Supplement 23, a wide range of SR documents templates have been standardized to address specific needs in the domains of CAD (e.g. Chest CAD SR, Mammo CAD SR), ultrasound procedure reports (e.g. Ultrasound OB-GYN Procedure report, Vascular Ultrasound Procedure Report, Echocardiographic Procedure Report, Intravascular Ultrasound Procedure Report) and CT/MR reports in cardiology (CT/MR Cardiovascular Analysis Report), among others.

Web Access to DICOM Objects

Web access to DICOM persistent objects (WADO) was introduced in the standard to facilitate the referencing and retrieval -via the internet protocols http and https -of DICOM persistent objects (images, structured reports, etc.), using URL/URI (Uniform Resource Locator/Identifier). This extension was made in 2003 with Supplement 85, in the form of Part 18 of the standard. It was also recognised by the ISO TC 215 (ISO 17432) [START_REF]Health Informatics -Messages and communication -Web Access to DICOM persistent objects[END_REF].

CONCLUSION

The development of the DICOM standard allowed to turn the page of proprietary standards for representing biomedical images and associated metadata. Initial efforts concerned primarily the exchange and management of images acquired on imaging equipment. Image processing is now being considered as well, with extensions of the standard dedicated to radiation therapy, CAD structured reports, spatial registration and fiducials, segmentation objects and even more recently segmentation surfaces such as meshes.

Besides, DICOM is open to the web technology to communicate and refer to images (or any other composite objects), as demonstrated with the WADO extension. However, DICOM remains based on its original binary syntax, in spite of growing pressure to move to some sort of XML encoding. For example efforts are being devoted in Working Group 20 to facilitate the translation of DICOM SR documents into HL7 Clinical Document Architecture documents, based on XML.

Finally, DICOM has started studying how processing tools, such as plugins or web services, could be more easily shared in the future, based on standard Application Programming Interfaces (API). This initiative, led by Working Group 23 "Application Hosting" received a strong support from the CaBIG initiative (Cancer Biomedical Informatics Grid) and the NCI, because of the perspectives it opens for in vivo imaging in cancer research, especially in molecular imaging, and CAD.

7.

Figure 1 .

 1 Figure 1. The different Parts of the DICOM standard (numbered 1 to 18)

Table 1 .

 1 the various working groups of the DICOM Committee

	DICOM Committee Working Groups
	WG-01: Cardiac and Vascular	WG-14: Security
	Information
	WG-02: Projection Radiography	WG-15: Digital Mammography
	and Angiography	and CAD
	WG-03: Nuclear Medicine	WG-16: Magnetic Resonance
	WG-04: Compression	WG-17: 3D
	WG-05: Exchange Media	WG-18: Clinical Trials and
		Education
	WG-06: Base Standard	WG-19: Dermatologic Standards
	WG-07: Radiotherapy	WG-20: Integration of Imaging
		and Information Systems

Table 2 .

 2 The various types of content items of SR trees

	Similarly, DICOM introduced a limited number of relationships, with
	relatively precise semantics such as 'CONTAINS', 'HAS OBSERVATION
	CONTEXT', 'HAS ACQUISITION CONTEXT', 'HAS PROPERTIES',
	'INFERRED FROM', 'SELECTED FROM' and 'HAS CONCEPT
	MODIFIER'.

ACKNOWLEDGEMENTS

Many thanks to Joël Chabriais, main representative of the French Society of Radiology in the DICOM Standards Committee, and currently radiologist in the hospital of Aurillac (France), for his continuous and enthusiastic participation to the DICOM standard's evolution and promotion.


Bastien X Nony

Mélanie C Rochoux

Didier Lucor

Thomas Jaravel

COMPOUND PARAMETRIC METAMODELLING OF LARGE-EDDY SIMULATIONS FOR MICROSCALE ATMOSPHERIC DISPERSION

Keywords: Computational Fluid Dynamics, Large-Eddy Simulation, Atmospheric Dispersion, Microscale, Surrogate Models, Parametric Uncertainty

In pollutant dispersion problems, mapping concentrations in the first tens or hundreds of meters from the source still remains a modelling challenge. Large-eddy simulations (LES) are able to represent time and space variability of turbulent atmospheric flow, which is of prime importance to assess public short-term exposure. However, they remain far from real time and subject to uncertainties, in particular to parametric uncertainties associated with the large-scale atmospheric forcing and the emission source position. In this work, we show that an efficient and accurate metamodel of the tracer concentration information provided by LES and encapsulating their associated uncertainties can be built using appropriate statistical tools combining machine learning and principal component analysis. We present a proof-ofconcept study based on a simplified but representative flow configuration (two-dimensional flow around a surfacemounted cube) using the AVBP LES solver and testing a variety of metamodels (linear regression, Gaussian processes, random forest, gradient boosting, etc.). Results reinforce the idea that for sufficiently statistically-converged quantities of interest and for a sufficiently large LES data set, a compound surrogate model can succeed in synthesizing information from the LES in the whole computational domain (with a Q 2 predictivity coefficient above 90 %). Downstream of the obstacle, the Q 2 coefficient of all metamodels reaches excellent results over 90%. Upstream, the tracer concentration is subject to strong discontinuities; combining metamodels allows to guarantee a good predictivity coefficient over 75%.

INTRODUCTION

Industrial accidents often involve a pollutant plume dispersion potentially harmful to human health, economy and/or environment. A new generation of decision support tools at the interface between modeling and statistical learning could be designed to help monitoring an emergency by accurately simulating the plume dispersion and accounting for a range of possible scenarios at future time frames.

Near-field plume dispersion is controlled by complex turbulent flow dynamics enhanced by complex geometry, for example in urban areas where separation and recirculation zones are induced by the presence of buildings of different height and geometry. Computational Fluid Dynamics (CFD) approaches, based on Reynolds-averaged Navier Stokes (RANS) and Large-Eddy Simulation (LES), is now becoming popular to accurately represent these processes (Philips et al., 2013). However, CFD approaches remain computationally intensive and therefore far from real time. Replacing the CFD model by a metamodel (i.e. a statistical model that has learnt the relationship between uncertain inputs and some quantities of interest over a CFD data set) can partly overcome this limitation. García-Sánchez et al. (2017) have demonstrated the capacity to train a metamodel based on polynomial chaos expansion, which mimics the RANS model response for a wide range of parametric variations on three uncertain parameters: the upstream roughness height, the wind direction and magnitude. Once the training step is done using a RANS data set, the resulting metamodel can be used to estimate at no further cost, the metamodel model response for a new set of parameters. This was illustrated for the "Joint Urban 2003 field experiment" in Oklahoma City, which is fully representative of urban geometry complexity. The objective of this work is to provide a detailed comparison of metamodels to identify which are the most relevant for pollutant dispersion and in particular if polynomial multi-linear regression approaches as in García-Sánchez et al. (2017) are sufficient to reproduce the CFD response to parametric variations. We consider here a two-dimensional canonical case (i.e. a two-dimensional flow around a surface-mounted cube) to compare a variety of metamodels, evaluate their performance and assess their robustness with respect to the training set size and data noise (associated with the time-average process). The number of uncertain parameters is kept low: three scalar parameters are considered: the inlet wind speed and the emission source position. The quantity of interest is the time-averaged tracer concentration affected by unsteady flow features, such as vortex shedding. The numerical data set is built using the AVBP LES solver (http://www.cerfacs.fr/avbp7x/). In the present study, the objective is to design a robust and accurate metamodel to represent the AVBP response to variations in the three uncertain parameters. This paper is structured as follows: i) the numerical case study and setup are introduced; ii) the metamodelling strategy is presented; and iii) the metamodel results are analysed.

NUMERICAL CASE STUDY AND SETUP

The test case is representative of a passive tracer dispersion in the lower part of the atmospheric boundary layer, which is modelled as an open two-dimensional computational domain in the reference frame (x, z) (Figure 1). A single obstacle with side length H = 1 m is considered. The domain consists of a rectangle 31H long (inlet flow direction x) by 10H high (vertical direction z). The domain height of 10H was chosen to avoid upper boundary effects. AVBP is used to solve the compressible Navier-Stokes equations with an artificial compressibility approach and with a fine mesh grid made of 240,000 triangular elements. A uniform spatial resolution of 0.04H (4 cm) was adopted to solve the flow around the obstacle with 25 grid points in each direction. A third-order accurate advection of the Taylor-Galerkin family is used. NSCBC (Navier-Stokes Characteristic Boundary Conditions) are used to properly handle acoustics at the inlet and outlet. The upper boundary models an atmospheric free surface through a plane symmetry condition. The lower boundary (the ground surface and the obstacle) is modelled as an adiabatic wall with a zero velocity. A constant pressure condition is imposed at the outlet.

METAMODELLING STRATEGY

The objective in this work is to design a metamodel able to predict the mean (time-averaged) tracer concentration field (containing Nnodes mesh nodes) as a function of the three uncertain parameters:

(1)

M is referred to as the response surface where the quantities of interest are the scaled concentration coefficients at different mesh nodes obtained from the uncertain parameters. The scaled concentration is defined by K = c/c0 with c the time-averaged concentration and c0 the normalizing concentration:

(2)

Only the nodes where the ensemble mean tracer concentration is above 𝐾 ≥ 4.59 × 10 23 (10 -2 kg m -3 ) inside the box [-4, 10]H×[0, 4]H are considered in the metamodeling procedure, resulting in Nnodes≈16,000. Outside this box, concentrations are very small, falling below the solver numerical accuracy, and are thus discarded for metamodelling. To perform a detailed analysis of the metamodel performance, a very large ensemble of LES (n = 700) is generated. The response surface is sampled using a quasi Monte-Carlo approach, assuming uniform statistical distributions for the three parameters given the range of variations in Equation 1. A low-discrepancy Halton's sequence is adopted in this study to obtain a homogeneous sampling of the low-dimensional parameter space. AVBP is run for each sampled set of parameters. The ensemble is segmented into a training set (70% of the LES, i.e. 489) and a validation set (30% of the LES, i.e. 211). Each resulting mean tracer concentration field K is decomposed using Principal Component Analysis (PCA) to reduce the output dimension from Nnodes≈16,000 to 50 principal components (Figure 2). Tests have shown that this is enough to accurately represent the tracer concentration over the whole domain. To limit the computational cost related to the LES data set generation, the AVBP simulation time is adjusted according to the inlet flow velocity Uinlet. Higher inlet velocities imply faster asymptotic convergence of the mean tracer concentration. A LES lasts from 150 s to 1,600 s in physical time to reach approximately a fixed number of 30 vortex shedding based on an estimated Strouhal value of St = 0.02. Using this approach, running 100 LES requires about 40,000 CPU hours on CERFACS' supercomputer.

Different regression metamodels are trained to approximate the response surface in Equation 1 (via the SciKit Learn library) using the LES training data set: multi-linear regression with and without penalty (Ridge, LASSO/Least Absolute Shrinkage and Selection Operator, OMP/Orthogonal Matching Pursuit), decision trees (random forest and gradient boosting), Gaussian processes. All of them are put in competition to represent the relation between each PCA component and the three uncertain parameters. In the resulting compound metamodel, we only keep a unique "champion" metamodel per PCA axis, i.e. the one that achieves the best performance according to the Q 2 predictivity coefficient (Figure 3):

(3)

with pi(K) and pi(𝐾 7 ) the projected components on the ith PCA axis of the AVBP prediction K and a given metamodel prediction 𝐾 7 , respectively, over the validation data set. The Q 2 coefficient gives the percentage of the variance explained by the metamodel (a perfect metamodel features Q 2 = 1; in practice, a Q 2 coefficient above 90% is considered as excellent). In this work, the performance of the compound and standalone metamodels is quantified using the Q 2 metrics on tracer concentrations at the Nnodes nodes over the validation data set. It is therefore a spatial field. To help the analysis of the results, different Q 2 statistics are computed: the spatially-average value of the Q 2 coefficient over the computational domain (the global Q 2 ) but also over limited-area regions (upstream, around the obstacle and downstream).

RESULTS

Figure 4 presents the spatial field of the Q 2 coefficient obtained for the compound metamodel. Table 1 presents the associated Q 2 statistics per limited-area domain. Results show a good performance of the compound metamodel in the whole computational domain, except in some upstream areas where the Q 2 is locally very low. This is consistent with a global Q 2 coefficient higher than 94%. The best Q 2 performance is obtained in the area downstream (the downstream Q 2 statistics is above 97%), while the Q 2 statistics decreases when moving upstream (the upstream Q 2 statistics remains above 80%). The dispersion process occurring downstream makes the mean tracer concentration field smoother for all sets of uncertain parameters. This is easier to capture by the tested metamodels than upstream of the obstacle, where the mean tracer concentration fields locally feature high values in distinct areas that depend on the tracer emission source location. Because only few realizations are considered in this area in the validation data set, the associated variance is low, prediction errors are penalizing, resulting in a poor Q 2 performance. Table 1 highlights the added value of the compound metamodel compared to standalone metamodels (when all PCA components are predicted using the same metamodel). The compound metamodel obtains the best performance for all Q 2 statistics. There is no much difference in the downstream area (where all metamodels reach a Q 2 statistics above 97%). However, having a compound metamodel limits the drop in performance in the upstream area: linear regression without penalty and Gaussian process regression feature an upstream Q 2 statistics equal or lower than 60%, while this statistics remains above 80% for the compound metamodel (gradient boosting follows closely at 78%). For the 90% quantile of concentration maximum residual errors, the compound metamodel value is K = 3.94, which is better than gradient boosting (4.49), linear (5.64) and Gaussian process regressors (6.60). Table 1 also shows that the performance of the compound metamodel remains the best among all tested metamodels and relatively good when reducing the size of the training set (from 489 to 100 LES) or when introducing noise in the training set (by dividing by two the AVBP simulation time). There are now more differences between the performance of the compound metamodel and that of the gradient boosting metamodel, especially in the upstream area. The compound metamodel good behaviour is explained by its flexibility in adjusting to PCA components that are difficult to predict. The PCA components 1 to 4 are best predicted using Gaussian process regression, while the PCA components 5 to 11 are better predicted by gradient boosting. This implies that Gaussian processes better represent large-scale tracer concentration structures, while gradient boosting is more adapted to track sharp concentration gradients near the tracer emission source position. Figure 5 gives an example of the compound metamodel solution for a validation case. The predicted mean tracer concentration is compared to the AVBP solution. Downstream plume dispersion is well predicted. The upstream tracer concentration is slightly under-predicted in the wake of the emission source. There is also some noise due to accumulated prediction errors on the high-level PCA components upstream.

CONCLUSION

The compound metamodel approach has shown its ability to reconstruct time-averaged information from LES on a simplified pollutant dispersion case. Mixing different metamodels (Gaussian process and gradient boosting) allows more PCA components to be considered upstream of the obstacle to improve the metamodel prediction and robustness to the training set size and noise. Future work includes adding a temporal dimension to fully take advantage of the LES, and applying it to more realistic cases. Philips, D.A., Rossi, R. and G. Iaccarino, 2013: Large-eddy simulation of passive scalar dispersion in an urban-like geometry. J. Fluid Mech.,723,

Figure 1 .

 1 Figure 1. Scheme of the two-dimensional computational domain with boundary conditions. Example of a mean tracer concentration field is given in the background for a source centered at (-2.26, 0.98) and an inlet wind speed of 5.58 m s -1 . Associated instantaneous snapshots of the tracer concentration fields are also given.

Figure 2 .

 2 Figure 2. Scheme of the LES data set generation and PCA postprocessing before metamodeling.

Figure 3 .

 3 Figure 3. Scheme of the metamodeling procedure.

Figure 4 .

 4 Figure 4. Spatial field of the Q 2 predictivity coefficient (Equation3) using the compound metamodel.

Figure 5 .

 5 Figure 5. Mean tracer concentration field for Uinlet = 5.58, xsrc = -2.26, ysrc = 0.98 for (top panel) the AVBP solution and (bottom panel) the compound metamodel prediction.

Table 1 .

 1 Comparison of the metamodel Q 2 statistics for a validation data set made of 211 LES. In brackets are indicated the Q 2 statistics obtained when dividing by two the AVBP time window, in square brackets the statistics obtained when the training data set is reduced from 489 to 100 LES.

	Metamodel/Q 2	Global Q 2	Upstream Q 2
			x

∈ [-4, -0.5]H Obstacle Q 2 x ∈ [-0.5, 1.5]H Downstream Q 2 x ∈ [1.

5, 10]H Compound 94% (93%) [85%] 80% (78%) [62%] 94% (93%) [80%] 97% (96%) [91%] Linear regression (without penalty) 90% (89%) [66%] 60% (60%) [-43%] 83% (82%) [44%] 97% (96%) [93%] Gradient boosting 93% (92%) [80%] 78% (78%) [52%] 92% (92%) [66%] 97% (95%) [88%] Gaussian process regression 87% (86%) [2%] 41% (40%) [-320%] 81% (81%) [-37%] 98% (96%) [80%]