{'fraction_non_alphanumeric': 0.028038484193991752, 'fraction_numerical': 0.010602788140585117, 'mean_word_length': 4.68691380080393, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Blood flow analysis in carotid artery bifurcation by two-dimensional ultrasonic-measurement-integrated simulation

If highly precise elucidation of the blood flow characteristics in a carotid bifurcation was possible, it would be widely applicable to diagnosis of circulatory diseases such as arteriosclerosis and cerebrovascular disease. This study was conducted to establish a new flow-dividing ratio estimation method applicable to an unsteady flow on a two-dimensional ultrasonic-measurement-integrated simulation of a carotid artery bifurcation for which it has been previously difficult to obtain a stable solution. In this new method, the flow-dividing ratio was directly adjusted by specifying the flow rate in a branch so that the difference of the Doppler velocities in the external carotid artery was decreased. The effectiveness of the proposed method was confirmed by a numerical experiment using the actual shape of a carotid artery bifurcation, and the superiority of the two-dimensional ultrasonic-measurement-integrated simulation over the ordinary simulation in terms of the reproducibility of the blood flow structure was clarified by analysis using clinical ultrasound data.


Introduction
Extensive research on arteriosclerosis has shown that it is closely related to hemodynamics (Schmid-Schönbein and Granger, 2003). In particular, it is known that a carotid bifurcation is the most common site of arteriosclerosis and stenosis. Since the stress acting on the carotid sinus in the internal carotid artery just before the bifurcation changes spatiotemporally and complexly (Hoi, et al., 2010;Lee, et al., 2009), this stress has a serious influence on progression of arteriosclerosis. It has been said that atherosclerosis is likely to develop in a part subjected to the low wall shear stress with reduction of substance transport between blood and the vessel wall by low velocity or stagnation of the blood flow (Chatzizisis and Giannoglou, 2006). However, it turns out that such development cannot be explained by such a simple mechanism due to advances in methodologies of experiment and numerical analysis, and it is not clear whether the wall shear stress acting on a vascular endothelial cell or the wall shear rate related to the substance permeability is the most critical factor. Elucidation of the pathologic condition of circulatory disease and its cause will possibly lead to early treatment or development of a new method of treatment. Also, identification of the blood flow structure in a carotid bifurcation is possibly applicable to diagnosis of cerebrovascular disease. That is to say, an internal carotid artery in a carotid bifurcation is connected to the cerebral artery, and it is suggested that the ratio of the flow rate between the external carotid artery and the internal carotid artery changes when blood vessel resistance arises due to cerebrovascular diseases (Kawai, 1981). For these reasons, if highly precise elucidation of the blood flow characteristics in the carotid bifurcation was possible, it would be widely applicable to diagnosis of circulatory diseases.
On the other hand, in current medical treatment, measurement of a blood flow in a carotid artery by an ultrasonic diagnostic device is widely used for diagnosis of circulatory diseases. Although blood vessel shape and Doppler velocity, the component of the blood flow velocity along an ultrasound beam, are obtained in real time by this device, the detailed structure of the blood flow is not. Funamoto, et al. proposed a two-dimensional ultrasonic-measurement-integrated simulation (hereinafter referred to as "two-dimensional UMI simulation") (Funamoto, et al., 2005) as a method to determine the detailed structure of an intravascular blood flow. This is a numerical simulation method, in which the difference of the Doppler velocities between the ultrasonic measurement and the numerical simulation is fed back to the simulation to reproduce the information on intravascular blood flow field correctly and in detail. The common carotid artery without branches was analyzed using this method, and the possibility of clinical application was shown (Funamoto, et al., 2011).
As for the two-dimensional UMI simulation of a carotid bifurcation blood flow, obtaining the ratio of the blood flow towards the branched blood vessel stably has been difficult, and blood flow has only been analyzed under ideal conditions. Although preliminary analysis of the two-dimensional UMI simulation of a steady flow in a simple-shaped branch model has been conducted to estimate the flow-dividing ratio (ratio between the flow rate in a branch of external carotid artery and that in a common carotid artery) by correcting the resistance of the downstream boundary in a previous study, this method is difficult to apply to an unsteady flow (Funamoto, et al., 2010).
Therefore, in this study, we attempted to establish a new flow-dividing ratio estimation method applicable to an unsteady flow. In this method, the flow-dividing ratio is directly adjusted by specifying the flow rate in a branch so that the difference of the Doppler velocities in a blood vessel behind the branch is decreased. The effectiveness of the proposed method was confirmed by a numerical experiment using the actual shape of a carotid bifurcation. The superiority of the two-dimensional UMI simulation over the ordinary simulation was considered in terms of the reproducibility of the blood flow structure by analysis using clinical ultrasound data.

Two-dimensional UMI blood flow simulation system
An outline of the two-dimensional UMI blood flow simulation system is herein presented. A block diagram of the system is shown in Fig. 1. Color Doppler measurement of a carotid bifurcation is performed with an ultrasonic diagnostic device (LOGIQ7, GE Healthcare Japan, Japan) with a linear type ultrasonic probe (10L, GE Healthcare Japan, Japan), and the data are saved in DICOM format. The measurement data are transmitted to a server (Altix UV1000, SGI, USA) through a memory card. A blood vessel shape is extracted by binarization between the blood Fig. 1 Schematic of a two-dimensional UMI blood flow simulation system. domain and the blood vessel domain from the time-averaged color Doppler image in the server. Numerical simulation of the blood flow is performed based on this vessel shape by correcting the velocity vectors through a feedback signal proportional to the difference between the calculation results and the measurement data. Iterative calculations of the flow analysis are performed until the inflow rate and the flow-dividing ratio are estimated. After the convergence of the solution, the velocity vectors, which cannot be acquired in measurement, are obtained and visualized on the graphics workstation.

Analysis object
This study dealt with blood flow in a human carotid bifurcation. The clinical data used for the analysis were color Doppler images of the longitudinal section in a carotid bifurcation of a healthy male volunteer in his forties (Fig. 2). The acquisition of the ultrasonic color Doppler images was conducted in accordance with the guideline of GE Healthcare Japan Corporation on clinical studies. In Fig. 2, the monochrome domain is a tissue image displayed in B mode, and the domain colored red and blue is the blood flow part. The blood is flowing from left to right. The straight blood vessel on the left side is the common carotid artery (CCA), and the two branches on the right side are the external carotid artery (ECA) and the internal carotid artery (ICA). As for the measurement conditions, the number of ultrasound beams was 288, the beam interval was 135 m, the spatial resolution in the depth direction was 96 m, and the time resolution was 86 ms. Measurement was conducted at a center frequency of 5 MHz, a pulse repetitive frequency (PRF) of 3.3 kHz, and a wall filter frequency of the minimum setup of the measurement device. The incidence angle of the ultrasound beam was 0 deg to the y-axis, the number of measurement time steps was 32 (in four cardiac cycles), and the measurement time interval t was 0.086 s.

Computational grid generation
For computation, the x-y coordinate shown in Fig. 2 was defined. The blood vessel shape was extracted from the time series data of Fig. 2, and a computational grid system was generated. The blood vessel shape and the orthogonal grid system are shown in Fig. 3. It is noted that the blood vessel shape extracted from Fig. 2 was rotated in counter clockwise direction with an angle of 8 deg so that the common carotid artery (CCA) was parallel to the x-axis. Resultant practical incidence angle of the ultrasound beam is 8 deg (Fig. 3). As for the blood vessel shape extraction, the blood domain and the blood vessel domain were binarized from the time-averaged color Doppler image. The number of grid points of a computational grid system was N x ×N y = 188×100, and the grid point intervals were set to x = 205 m and y = 194 m, respectively. Assuming human blood for the fluid, the density was set to =1.00×10 3 kg/m 3 and the kinetic viscosity was set to  = 4.0×10 -6 m 2 /s. As to the characteristic value of nondimensionalization in the flow computation, the characteristic length was set to the vessel diameter 4.86×10 -3 m at the upstream boundary, and the characteristic velocity was arbitrarily set to 0.1 m/s. The fluid density was also used as the characteristic value.

Numerical analysis method
For the purpose of practical clinical applications, this study dealt with the two-dimensional flow on the two-dimensional longitudinal section though the actual intravascular blood structure is three dimensional. Also, neglecting the elastic deformation of the blood vessel, the calculation was performed for the rigid blood vessel wall. In this case, the governing equations of the UMI simulation are two-dimensional unsteady incompressible Navier-Stokes equations (1) and the pressure equation (2).
Equation (2) is obtained by calculating the divergence of Eq. (1) and substituting the continuity equation into that equation. In the above expressions, u = (u, v) is the velocity, t is the time, is the viscosity, is the density, and p is the pressure. f is the external force term of the feedback signal given in the following equation proportional to the difference between the Doppler velocities, V m of the ultrasonic measurement and V c of the numerical calculation, and is applied in the feedback domain within the calculation domain.
where K V * is the feedback gain (non-dimensional), U is the characteristic velocity, L is the characteristic length and b is the unit vector along the ultrasound beam. Positive Doppler velocity and the orientation of the unit vector b correspond to the direction with increasing y coordinate. The case with K V * = 0 corresponds to an ordinary simulation without feedback.
The above governing equations are discretized based on a finite volume method using a staggered grid system, and the resultant difference equations are solved with a method similar to the SIMPLER method (Hayase, et al., 1992). The feedback domain was domain C from 1/8 to 7/8 from the upstream boundary within the calculation domain shown in Fig. 3.

Fig. 3
Extracted shape of longitudinal section of the carotid bifurcation and computational grid.

Inflow rate and flow-dividing ratio estimation method
The estimation methods of the inflow rate and the flow-dividing ratio of the boundary condition of the flow analysis are now described. The flow-dividing ratio is defined in Eq. (4) as the rate between the flow rate in ECA (q 1 ) and the inflow rate (q 0 ). Note that "flow rate" stands for "flow rate on the two-dimensional longitudinal section" (the unit is [m 2 /s]) in this paper intended for a two-dimensional blood flow field.
The inflow rate estimation in the upstream CCA part is now explained. The golden section search (GSS) (Press, 1986) is employed to search for the upstream inflow rate q 0 so that the sum of the absolute values of Doppler velocity errors for the measurement value (Eq. (5)) is minimal in the domain A from 1/8 to 3/8 from the upstream boundary in the calculation domain shown in Fig. 3 (Yamagata and Hayase, 2009) as follows: is the inflow rate of the upstream boundary. N A is the number of grid points in domain A. The flow-dividing ratio estimation of the downstream boundary is explained below. The fixed-point iteration method (FIM) (Burden and Faires, 1985) was used to search for the flow-dividing ratio so that the sum of the differences for the measured and calculated Doppler velocities (Eq. (6)) is 0 in domain B from 6/8 to 7/8 from the upstream boundary within the calculation domain shown in Fig. 3.
The index i for the flow-dividing ratio r is the number of the iteration of FIM, N B is the number of grid points in domain B, and is the relaxation coefficient.
The downstream boundary condition was given by the flow-dividing ratio r i obtained by Eq. (6) as follows: where u e is the velocity at each grid point at the downstream boundary, u e-1 is the velocity at the upstream grid next to the downstream boundary, q 2 is the flow rate in the longitudinal section in ICA, and  is the relaxation coefficient.
(q 1 +q 2 )/q 0 is the term to prevent divergence in an early stage of iterations.

Calculation flow chart
The above calculation process is shown by the flow chart in Fig. 4. First, the estimation of the inflow rate q 0 is performed concerning the velocity field at the first time step. The flow rate q 0 is updated based on GSS, and the velocity and the pressure are calculated by the flow analysis described in the section 2.4. After the residual at convergence for the flow rate is satisfied by iterative calculation, the estimation of the flow-dividing ratio r is performed with q 0 fixed. The velocity and the pressure are calculated by flow analysis, updating flow-dividing ratio r based on FIM, and the solution satisfying the convergence criteria is regarded as the final solution of the one-time step. Note that these calculations are performed applying the feedback signals. Furthermore, the blood flow fields after the next time step are calculated in the same process. When the acquisition of desired blood flow fields in a time series is completed, the analysis is finished.

Calculation conditions
In this study, both a numerical experiment using the Doppler velocity of a numerical solution calculated in advance, and analysis using clinical data were conducted. First, the calculation conditions of the numerical experiment are described. The object of analysis is the numerical solution calculated in advance as a blood flow model (called as the standard solution). The calculation conditions of the standard solution, the UMI simulation and the ordinary simulation are summarized in Table 1. Note that the ordinary simulation corresponds to the UMI simulation with the null feedback signal (K V * = 0). As for the standard solution, its flow-dividing ratio was set to r = 0.5, and its time-variation of the mean flow velocity at the upstream boundary was the sinusoidally oscillating flow (1 Hz) with a steady flow component. Its velocity field in the systolic phase, whose mean flow velocity at the upstream boundary has a maximum value 0.45 m/s, is shown in Fig. 5. The upstream velocity profile was assumed to be parabolic. The convergence criteria for flow analysis as well as GSS and FIM are shown in Table 1. In calculation of UMI and ordinary simulations, the relaxation coefficients in Eq. (6) and Eq. (7) were determined by preliminary calculation. The initial flow-dividing ratio of UMI and ordinary simulations was set to 0.2, close to that of a healthy adult (0.3) (Swillens, et al., 2009). Assuming that the actual velocity profile is unknown, the upstream velocity profile of UMI and ordinary simulations was given as the uniform one different from that of the standard solution.
Next, the analysis using the clinical data is described. The object of analysis is in vivo blood flow field for which the ultrasonic measurement was performed. The UMI and ordinary simulations were performed using the Doppler velocity measurement. Calculation conditions were identical to those of the numerical experiment given in Table 1.

Evaluation of calculation accuracy
In order to evaluate the calculation accuracy of the two-dimensional UMI and ordinary simulations, Doppler velocity error norm e V (non-dimensional) was defined as follows: where N is the total grid number of estimation domain, n is the grid point index, and V type is the literature data of the maximum mean Doppler velocity in CCA (0.39 m/s (Kitaoka, et al., 2009)). The velocity vector error norm eu (non-dimensional) was also defined as follows:  (Funamoto et al., 2011) was observed. The errors of the Doppler velocity, the velocity vector and the flow-dividing ratio at K v * = 200 decreased to 26%, 43% and 11% of those of the ordinary simulation at K v * = 0, respectively.
The numerical experiment clarified that the errors of the velocity vector and the flow-dividing ratio could be reduced by reduction of the Doppler velocity error, and the effectiveness of the UMI simulation method in a carotid bifurcation was confirmed.

Analysis using clinical data
The calculation result of the analysis using actual clinical data is secondly described. Variation of the time averaged value of the Doppler velocity error norm e _ V with feedback gain is shown in Fig. 9. Similar to the result of the numerical experiment, the Doppler velocity error norms at K v * = 200 decreased to 33% of that of the ordinary simulation. The Doppler velocity error becomes smaller as the gain is larger in the result, but the noise component is reproduced, and the flow field is unnatural in accordance with an increase of gain, in general. Therefore, the result of the UMI simulation at K v * = 200 is investigated in comparison with that of the ordinary simulation below.
The time-averaged spatial distribution of the Doppler velocity error norm is shown in Figs. 10 for UMI and ordinary simulations. This result is the time-averaged value for three cardiac cycles except for the first two cardiac cycles at the start of calculation. In the UMI simulation, the error in the feedback domain is reduced in comparison with that of the ordinary simulation, but the reduction is relatively insignificant in ICA. As the cause of this, it is considered that the ICA is not contained in the estimation domain of the flow-dividing ratio and that the clinical data reflect the three-dimensional actual blood flow whose flow rate on the two-dimensional longitudinal section does not satisfy q 0 = q 1 +q 2 .
Time variation of the upstream mean flow velocity u ave , the flow-dividing ratio r and the Doppler velocity error norm e V are shown in Figs. 11. The upstream mean flow velocity u ave corresponds to the inflow rate q 0 as it is expressed by u ave = q 0 /D 0 , where D 0 is the blood vessel width. The time variation of u ave of the UMI simulation is mostly periodic, while that of the ordinary simulation has larger amplitude than that of the UMI simulation and is poor in periodicity.  (c) show the results at the systolic phase (11th step, t = 0.95 s), the end-systolic phase (13th step, t = 1.12 s) and the diastolic phase (17th step, t = 1.46 s), respectively. At any phase, the Doppler velocity distribution of the UMI simulation at the middle is close to that of the clinical data at the top in the feedback domain implying the in vivo complex velocity field is properly reproduced, while that of the ordinary simulation at the bottom is different from that of the clinical data.
Looking at each time in detail, at the systolic phase in Fig. 12(a), a vortex is generated on the part where the vessel diameter expands towards the bifurcation area from CCA by x = -0.005 in the UMI simulation. In the ordinary simulation, the Doppler velocity distribution is almost uniform from the upstream to the branch, and the vortex shown in the UMI simulation does not exist. At the end-systolic phase in Fig. 12(b), in the UMI simulation, the velocity in the lower internal carotid sinus is small and stagnation is generated. In the ordinary simulation, a vortex is generated in the lower internal carotid sinus and the bifurcation area in front of ECA. At the diastolic phase in Fig. 12(c), although the total velocity is small and little difference between the UMI simulation and the ordinary simulation is found, the velocity in the lower internal carotid sinus in the UMI simulation is small and stagnation is generated similar to the result at the end-systolic phase. This result shows the same tendency as the blood flow measurement result of the carotid bifurcation of MRI (Markl, et al., 2010). In the ordinary simulation, stagnation is small in the lower internal carotid sinus. According to the above results using clinical data, spatiotemporally different blood flow fields between the UMI simulation and the ordinary simulation were obtained, and the UMI simulation reproduced the stagnation generated in the internal carotid sinus, a common site of arteriosclerosis.
Limitations of this paper are discussed below. A real blood flow field is three-dimensional, containing the secondary flow on a cross section perpendicular to the vessel axis, and generally it cannot be exactly reproduced by two-dimensional analysis though the flow field on a two-dimensional longitudinal section was well reproduced by numerical analysis with feedback from ultrasonic measurement. Additionally, the flow-dividing ratio of a real blood vessel calculated from three-dimensional velocity vectors in a cross section perpendicular to the vessel axis is probably different from the present one evaluated from the two-dimensional velocity vectors on the two-dimensional longitudinal section. However, it is considered that the present flow-dividing ratio has a value relatively close to flow by choosing the domain of flow rate estimation in the main and branch vessels so as to contain the vessel axes as much as possible. In future work, a more precise estimation of the flow-dividing ratio will be possible in the framework of a three-dimensional UMI simulation. Moreover, comparison of the result of two-dimensional UMI simulation with those of MRI measurement, a model experiment, and a numerical experiment will be necessary to verify the validity of the present method. An appropriate correlation between the flow-dividing ratio of the present two-dimensional analysis and that of the real three-dimensional blood flow field is crucial for usage in clinical practice.
In order to perform blood flow analysis easily in a clinical setting, the two-dimensional UMI simulation system in this study deals with the blood vessel configuration reconstructed on the orthogonal grid system by binarizing the ultrasonic measurement data. Therefore, the shape of the blood vessel wall is stepwise, as shown in Fig. 3, and there is a region with relatively large error near the wall, as shown in the analysis result for clinical data in Fig. 10. Additionally, (a) Systolic phase (t = 0.95s) (b) End-systolic phase (t =1.12s) (c) Diastolic phase (t =1.46s) Fig. 12 Doppler velocity distribution and velocity vector field in the analysis using the clinical data.
an appropriate smoothing process is essential in calculation of the wall shear stress. However, it is considered that the effect of the stepwise wall shape on the result of analysis of the flow field is limited near the vessel wall since the spatial scale of the error occurring near the wall is small. In this work, the blood vessel was treated as a rigid body, but a real blood vessel has elasticity, and it deforms in response to pulsation. According to the literature (Robert et al., 2004), the change of the vessel diameter by pulsation is about 5%. In the present UMI simulation, the velocity field of the blood flow was properly reproduced, but the error occurred in the blood vessel diameter as the time variation of the blood vessel shape was ignored. As a result, errors in the flow rate and the flow-dividing ratio were introduced.
Finally, the present method is applicable for blood flow in a bifurcation for which ultrasonic measurement is possible, such as a renal artery and a portal vein as well as a carotid bifurcation.

Conclusions
This study aimed to establish a new flow-dividing ratio estimation method applicable to an unsteady flow on a two-dimensional UMI simulation of a carotid artery bifurcation. In this method, a flow-dividing ratio was directly adjusted by specifying the flow rate in a branch so that the difference of the Doppler velocities in a blood vessel behind a branch was decreased. The effectiveness of the proposed method was confirmed by a numerical experiment using the actual shape of a carotid bifurcation and the superiority of the two-dimensional UMI simulation over the ordinary simulation in terms of the reproducibility of the blood flow structure was clarified by analysis using clinical ultrasound data. In future work, the effects of the three-dimensionality of a blood vessel, the stepwise vessel shape on the orthogonal grid, and the elasticity of the blood vessel will be considered. It is noted that the present method is applicable for blood flow in a bifurcation for which ultrasonic measurement is possible.