Abstract
Purpose Evaluation of hemodynamics is crucial to predict growth and rupture of cerebral aneurysms. Variational data assimilation (DA) is a powerful tool to characterize patient-specific intra-aneurysmal flows. The DA method inversely estimates a boundary condition in fluid equations using personalized flow data; however, its high computational cost in optimization problems makes its use impractical. This study proposes a practical DA approach to evaluate patient-specific intra-aneurysmal flows.
Methods To estimate personalized flows, a variational DA method was combined with computational fluid dynamics (CFD) analysis and observed intra-aneurysmal velocity data, and an inverse problem was solved to estimate the spatiotemporal velocity profile at a boundary of the aneurysm neck. To circumvent an ill-posed inverse problem, model order reduction based on a Fourier series expansion was used to describe temporal changes in state variables.
Results In numerical validation using synthetic data from a direct CFD analysis, the present DA method achieved excellent agreement with the ground truth, with a velocity mismatch of approximately 18%. In flow estimations for three patient-specific datasets, the velocity mismatch for the present DA method was markedly lower than that for the direct CFD analysis and would mitigate unphysical velocity distributions in flow data from phase contract magnetic resonance imaging.
Conclusion By focusing only on the intra-aneurysmal region, the present DA approach provides an attractive way to evaluate personalized flows in aneurysms with greater reliability than conventional CFD and better efficiency than existing DA approaches.
Introduction
A cerebral aneurysm is a cerebrovascular disorder in which part of vessel wall develops an outward bulge. Rupture of a cerebral aneurysm is a major factor in subarachnoid hemorrhage and has a mortality rate of up to 50% [1], although the rupture rate is only approximately 2% per year [2]. Since there is a trade-off between the risks of rupture and clinical intervention, a high-level evaluation to predict the rupture risk is highly desirable.
It has been argued that growth and rupture of an aneurysm are generated by a three-way relationship among the pathobiology of the wall, aneurysmal geometry, and intra-aneurysmal flow [3,4]. To quantify the hemodynamics of cerebral aneurysms, phase-contrast magnetic resonance imaging (MRI), also termed four-dimensional (4D) flow MRI, and computational fluid dynamics (CFD) have been conducted and compared [5,6,7]. However, these approaches cannot adequately quantify hemodynamic features because the spatiotemporal resolution of 4D flow MRI is limited and it is difficult to set patient-specific boundary conditions for CFD [8].
Data assimilation (DA) is a powerful tool to quantify time-dependent intra-aneurysmal flows in each patient. It predicts a more statistically reliable velocity field based on the observed and modelled velocities from 4D flow MRI and CFD, respectively. In particular, an inverse approach is preferably used to strictly satisfy physical constraints [9]. Several types of inverse DA algorithms have been developed for time-dependent problems; these algorithms estimate boundary conditions, such as the inlet/outlet velocity, using a feedback control method [10], variational method [11,12], and sequential (Kalman filter-based) method [13]. However, these DA approaches analyze the main vessels in addition to an aneurysm, resulting in a high computational cost with respect to the computation time and memory requirement for such time-dependent problems. Moreover, since solutions are often sensitive to material parameters such as blood viscosity and vessel shape reconstructed from image-segmentation, uncertainty estimates under many different analysis conditions are required for each patient, which further increases the computational cost. These difficulties should be mitigated so that DA can be used to quantify patient-specific intra-aneurysmal flows.
This study aimed to develop a practical DA approach to evaluate patient-specific intra-aneurysmal flows. Using a variational DA method coupled with CFD and measured velocity data, and limiting the domain of analysis to the inside of the aneurysm, an inverse problem was solved to estimate the spatiotemporal velocity profile at the aneurysm neck. Furthermore, to circumvent an ill-posed inverse problem, model order reduction based on a Fourier series expansion was employed to describe temporal changes in state variables. The model was validated numerically using synthetic data obtained from CFD with given boundary conditions, and then the model feasibility was investigated by applying it to three patient-specific MRI datasets.
Materials and Methods
Patient-specific data
Time-of-flight magnetic resonance angiography (TOF-MRA) and 4D flow MRI were conducted using three patients, each of whom had an aneurysm at an internal carotid artery. The study was approved by the ethics committees for human research (IRB Number: R2019-227). The TOF-MRA images (pixel size, 0.3906 mm; slice thickness, 0.6 mm) were used for reconstruction of vessel geometries. The spatiotemporal velocity fields of the blood stream were reconstructed from 4D flow MRI data (pixel size, 0.7031 mm; slice thickness, 1.0 mm; temporal resolution, 12 frames per heartbeat) acquired using velocity encoding (VENC) at 40 cm/s (VENC40) and 120 cm/s (VECN120).
Extraction and smoothing of vessel shapes were performed using Mimics (Materialise, Belgium) and Meshmixer (Autodesk, USA). For further usage, surface shapes of the vessels were implicitly described using a discrete level-set function with a signed distance. The 4D flow MRI and TOF-MRA datasets were co-registered using a rigid transformation based on global geometry information stored in the headers of their corresponding DICOM files. The velocity field was extracted from the 4D flow MRI data using VENC, and the velocity inside the vessel was then extracted using the level-set function for the vessel reconstructed from the TOF-MRA image (Fig. 1). Supplementary Fig. S1 shows the vessel shapes for each of the three patients. Only velocity vectors at positions greater than one pixel from the vessel wall were extracted to eliminate artifacts located near the wall.
Application of DA to cerebral aneurysms
4D variational (4D-Var) method
In this study, a 4D-Var method was applied for DA. In this method, the interior of an aneurysm without branching vessels is introduced as an analysis domain, and an optimization problem is solved to minimize velocity mismatch between the mathematical model and experimental observation in that domain. An optimal control concept is used to solve the problem by imposing a spatiotemporal velocity profile on the aneurysm neck as a boundary condition for the incompressible Navier–Stokes equations, and this velocity profile is inversely estimated as a design variable. The mathematical description can be given as
where j(u, g) is the objective function consisting of two cost functions; the first represents the total data mismatch between the model velocity u(x, t) and observed velocity 
in the spatiotemporal domain comprising the analysis space x ∈ Ω, observation space 
, and time 
, where T is the period of pulsatile flow; and the second represents a zeroth-order Tikhonov regularization function with a parameter ε for the design variable g(x, t) (i.e., the boundary velocity on the aneurysm neck Γbc) to mitigate sensitivity to data noise and an ill-posed inverse problem. Parameter Q is the operator mapping the solution field of the mathematical model to the observation. In Eq. (1), the incompressible Navier– Stokes equations act as constraints for the model velocity u by imposing a no-slip condition at the fixed wall Γwall and an inlet/outlet velocity g at the aneurysm neck Γbc. Here, p is the fluid pressure, ρ is the fluid density, η is the fluid viscosity, and n is the unit normal vector on the surface Γbc.
To simplify the numerical setup, the extracted aneurysms were rotated so that the aneurysm extrusion pointed along the z direction and the inlet/outlet surface Γbc was defined to be a cross-section through the neck in the xy plane (Fig. 2).
By applying the Lagrange multiplier technique, the minimization problem (1) can be rewritten
where λ, α, φ, and ξ are the Lagrange multipliers. The stationary conditions result in sets of partial differential equations (PDEs) for the original (state) system
and the adjoint system
where Q* is the inverse operator mapping the observation to the model field. Analogously, the functional derivative of j in terms of g is given by
where ∂n = n ⋅ ∇ is the normal derivative on the surface Γbc. In these derivations, a pulsatile periodicity ⋅|t=0 = ⋅|t=T is imposed. More details of the derivation are provided in the Supplementary Information (SI_text1).
Note that the equality H = 0 gives first-order optimality conditions of the optimization problem, and in general the above nonlinear systems are iteratively solved for the design variable g so that H → 0.
Optimization calculation
A gradient descent method was used to solve the optimization problem. From Eq. (5), the evolutional equation for g is then given by
where k enumerates the optimization step and Δτ is the hyperparameter in the gradient descent to control corrections in the optimization interval [k, k + 1]. The model and adjoint variables are obtained by solving systems (3) and (4) as F(uk, pk; gk) and G(λk, αk; uk, Uobs), where the pulsatile periodicity in flow should be sufficiently satisfied. System (4) was solved numerically in a time-reversed manner from t = T to t = 0.
Model order reduction
Focusing on the periodicity of the problem, model order reduction was implemented to express temporal changes in variables, including the design variable, in terms of Fourier series expansions. Let us denote the Fourier series expansion of g by
where Nl is the maximum expansion order, ω = 2π/T is the angular frequency, and al (x) and bl (x) are the Fourier coefficient distributions on Γbc. Denoting 
and 
, the mappings to evaluate coefficients in the Fourier cosine and sine series in Eq. (6), namely al (x) = Cl (g(x, t)) and bl (x) = Sl (g(x, t)), can be rewritten as
where sk = −αkn + η∂n λk. Thus, once the Fourier coefficients for the term sk are evaluated, al and bl can be updated using Eq. (8). In this study, Cl (sk) and Sl (sk) were calculated using trapezoidal numerical integration at the same time as solving the adjoint system (4).
The present method of model order reduction has three potential advantages. First, it mitigates an ill-posed inverse problem resulting from large degrees of freedom when directly updating g, which is discretely introduced in every time step of the fluid simulations. Second, the continuous profiles of the Fourier series smooth temporal changes. Third, the method dramatically reduces the memory required to store temporal changes in g and u when solving the adjoint system (4).
Mapping operators between modelled and observed data
The mapping operator Q was defined as
where the kernel D is a function of a parameter h, δ is the Dirac delta function, and ΔΩh is the size of small control domain xh ∈ Ω (in general, ΔΩh = h3). Note that the observation position it is not guaranteed to equal that of the model, and thus smoothing is required by applying the kernel at the resolution of the CFD mesh. Since a Cartesian mesh was used to solve the fluid equations, as described later, a three-dimensional smooth Dirac delta function was used for D [14,10]. Analogously, the inverse mapping operator Q* was defined as
where 
is the velocity mismatch between the model and observation, and 
and 
are the discrete position and time, respectively, associated with the observation data, namely 
and 
. Here, M and N are the number of observation data for space and time, respectively. The first term of the objective function j in Eq. (1) was approximated by
where 
is the observation domain size.
CFD solver
The original and adjoint systems in Eqs. (3) and (4), respectively, were solved using a Cartesian mesh-based CFD solver [10,15] based on the boundary data immersion method [16]. This method solves a unified PDE model for fluid and solid dynamics using a smooth characteristic (or density) function, so a simple fixed-mesh system such as the Cartesian coordinate system can be adopted. Analogous to the existing solver [10,15], a pressure projection method was used to couple the velocity and pressure (and the related adjoint variables) on a Cartesian mesh.
The following material parameters were used: ρ = 1050 kg/m3, η = 0.0035 Pa, and T = 480 ms. The dimensions of the mesh elements were set to Δx = Δy = Δz = 0.2 mm, and the time interval was set to Δt = 0.05 ms, yielding a Courant number of approximately 0.1 for sample B and 0.2 for samples A and C with respect to their corresponding maximum velocities. Simulations were executed on single CPU (A64FX, Fujitsu, Japan) of either the Fugaku (RIKEN Center for Computational Science, Kobe, Japan) or Flow Type I (Information Technology Center, Nagoya University, Japan) supercomputers.
Analysis cases
Direct simulation using main vessel branches and the observed inlet flow condition
To validate the present DA method, direct simulations were compared with numerical solutions. In the simulations, the aneurysm and the main vessel branches were included in the analysis domain. A Dirichlet boundary condition for the velocity was imposed on a cross section of an internal carotid artery inlet and prescribed a uniform spatial profile and a Fourier series-based temporal profile evaluated from 4D flow MRI data at VENC120 around the inlet region. The simulations were performed for five pulsatile periods, and the result at the final period was used for evaluation.
The observed velocities in the aneurysm region were smaller than those in the main branches and had peak magnitudes of less than 40 cm/s. Moreover, the velocity fluctuations resulting from VENC40 were also lower than those from VENC120 (Fig. S2). This implies that the intra-aneurysmal flow data at VENC40 are more reliable than those at VENC120 in terms of the signal-to-noise (S/N) ratio. Therefore, the inlet velocities used in the direct simulations were adjusted so that the velocity in the aneurysm was close to that measured using 4D flow MRI at VENC40.
Comparison of DA analyses using synthetic data and 4D flow MRI data
Eqs. (3) and (4) were evaluated using two pulsatile times that almost satisfied the flow periodicity, and the solutions at the second time were used for evaluation of the Fourier coefficients in Eq. (8). From preliminary trials, Nl = 5 was set for the Fourier-series expansion, Δτ = 1000 was set in the gradient descent, and the Tikhonov parameter was set to ε = 0 (i.e., no regularization was considered). The number of iterations of the optimization was set to 300. The parameter h for the kernel D mapping the velocity between the model and observation was set to h = Δx so that the kernel was smoothly distributed across neighboring mesh elements, and observation data at positions two or more mesh elements away from the wall and inlet/outlet boundaries was used.
For validation, DA analysis was performed in the aneurysm region using synthetic data created from a direct simulation by down-sampling its spatiotemporal resolution to match that of the 4D flow MRI data. The numerical results of the direct simulation were then used as the ground truth.
For practical applications, DA analyses were conducted using 4D flow MRI data acquired at VENC40 for three patients.
Metrics for velocity mismatch
Using the L2 norms
for spatial and spatiotemporal quantities, where qm = q(xm)|t and 
, normalized metrics for the velocity mismatch between the model and observation (Eu) were defined in space and time:
Other metrics for the velocity mismatch between the model and observation were also introduced with respect to the vector (Ru), magnitude (R|u|), and angle (R∠u):
Results
Validation using synthetic data
Fig. 3 shows estimated results of the DA analysis using synthetic data and the ground truth. The DA analysis well captured the velocity distributions in two cross sections. The velocity mismatch in each frame fell within the 16%–27% range. In the longitudinal cross section, some discrepancies between the DA analysis and the ground truth were observed around the aneurysm neck. The direct simulation as the ground truth resulted in high velocities at the aneurysm neck and a corresponding inlet jet along the side wall. However, the DA analysis did not adequately capture these features. This might explain the maximum mismatch at frame 10, at which time the jet-like inlet flow increased from a low-flow state as the systolic phase was approached.
Fig. 4 compares the velocity in the ground truth with that obtained from DA analysis, both of which are plotted at the observation resolution. Excellent agreement was achieved for velocity vectors in the systolic phase. Although the Bland–Altman plot shows there was some additive bias (mean) and fluctuation (standard deviation) in mismatches of the velocity magnitude, the sizes of these mismatches were much smaller than the velocity values themselves.
Comparisons with 4D flow MRI data
The intra-aneurysmal velocity vectors obtained from 4D flow MRI, the proposed DA analysis, and direct simulation at systolic phases for three patients are shown in Fig. 5. Although both the DA analysis and direct simulation could capture the overall behavior of the flow fields obtained from 4D flow MRI, there were several large discrepancies in the velocities between patients B and C. The mismatches (Eq. (12)) between the velocity obtained from 4D flow MRI and those obtained from the direct simulation and DA analysis are summarized in Table 1. For all three patients, the velocity mismatches between the 4D flow MRI and DA analyses were 38%–48% lower than those between the 4D flow MRI and the direct simulations. However, the velocity mismatches between the DA and 4D flow MRI analyses for patients B and C were approximately twice as large as the corresponding mismatch for patient A. The velocity magnitudes (i.e., the spatiotemporally averaged velocities using the L2 norm) for the 4D flow MRI were 0.132, 0.095, and 0.129 m/s for patients A, B, and C, respectively. Notably, the size of these values was in the order A > C > B, which followed the opposite trend of the velocity mismatches (A < C < B).
Bland–Altman plots of the agreement between the velocity magnitude obtained from 4D flow MRI and velocities obtained from the DA analysis and direct simulation are shown in Fig. 6. For all patients, velocity differences between the mathematical model and 4D flow MRI were much smaller for the DA analysis than for the direct simulation. Moreover, variation in the velocity differences among the patients was smaller for the DA analysis.
To further investigate velocity mismatch between the DA analysis and 4D flow MRI, metrics for vector, magnitude, and angle in Eq. (14) were evaluated (Fig. 7). The metric for the magnitude was always lower than that for the vector, and the angle metric was approximately 0.1–0.3 (i.e., 10°–54°). Again, the size of the metric mismatch for patients A–C was in the order A < C < B, which followed the opposite trend of the velocity magnitudes (A > C > B).
Effects of VENC values on the DA analysis
To investigate the influence of data noise on the velocity estimation, DA analysis using VENC120 data was performed for patient A. The velocity mismatches for both VENC40 and VENC120 are shown in Fig. 8. The velocity mismatch at VENC40 was clearly much smaller than that at VENC120.
Discussion
A novel numerical approach to quantify patient-specific intra-aneurysmal flows was developed using a 4D-Var DA method. The growth and rupture of cerebral aneurysms are influenced by hemodynamic factors and follow biophysical changes in the vessel wall [3], highlighting the importance of personalized predictions. The DA analysis is a key tool for such predictions because it can provide detailed hemodynamic factors, and the 4D-Var formulation is well-suited to predicting flow fields under unknown boundary and initial conditions. However, the computational cost of solving a constrained optimization problem is large and generally prohibits the clinical application 4D-Var DA. The proposed DA approach, which focuses on an intra-aneurysmal analysis domain, dramatically reduces the calculation time because the computational cost of forward simulations for the original and adjoint systems in every step of the optimization is comparatively low. Using model order reduction to describe temporal changes in variables offers further advantages by avoiding an ill-posed inverse problem and reducing the corresponding memory requirement. A similar idea to parametrize the spatiotemporal profile of the inlet velocity with a few control parameters was previously proposed [13]. However, the method was restricted to simple velocity profiles in both space and time, limiting its application to real-world datasets. Another advantage of the present Fourier series-based approach is that it does not require estimation of any initial conditions in time-dependent systems. In existing DA methods, main vessels and several branching vessels are included in the analysis domain, and thus the computational cost is large for every iteration. Moreover, an image processing-based extraction procedure for main vessels is sensitive to the CFD simulation [17,18,19]. The proposed DA approach bypasses this problem, at least for the main vessels including branches, which could reduce the computational burden of pre-processing procedures. For these reasons, the proposed DA approach is anticipated to become a practical clinical tool for estimating patient-specific inter-aneurysmal flows.
Using the synthetic data obtained from the direct simulation of flow velocity in an aneurysm and the main vessel branches with a boundary condition imposed on vessel geometry, the present DA method could reasonably reproduce the original (ground truth) flow field. The spatially averaged velocity mismatch was within approximately 16% in the systolic phase (until frame 7) but increased to approximately 27% between the end of diastole (frame 8) and the beginning of the next systole (frame 10). During these frames, a jet-like inflow from the aneurysm neck resulted in high flow instability and sudden changes in the flow magnitude and direction, which would lead to a large velocity mismatch. The objective function for the velocity mismatch did not converge to zero even in the DA analysis using the ground truth data. There are two possible reasons for this. First, the inlet/outlet boundary of the aneurysm was vertically extended from the neck to stabilize the numerical simulations, and this change in aneurysm shape may have caused a mismatch with the ground truth. A second possibility is that the solution reached a local minimum, and a different set of spatiotemporal profiles of the inlet/outlet velocity was obtained compared with the original profile. This is consistent with the observation that the estimation accuracy of the present DA method was lower around the aneurysm neck than inside the aneurysm. Notably, although both of the above scenarios may have occurred, reasonable flow estimations inside the aneurysm were nevertheless obtained.
In image-based blood flow simulations, hemodynamic factors inside an aneurysm are dependent on inflow boundary conditions [8]; nevertheless, setting complete spatiotemporal profiles on the inlet of main vessel is difficult owing to resolution and noise limitations in 4D flow MRI. Moreover, the flow field inside the aneurysm is strongly dependent on the accuracy of shape reconstruction of main vessels because the magnitude of the intra-aneurysmal flow is much smaller than that through the main vessels. This could explain why the direct simulation using main vessels reproduced the 4D flow MRI velocity with low accuracy. In contrast, the present DA method, where the spatiotemporal profiles on the inlet/outlet boundary of the aneurysm neck are automatically set by solving the optimization problem, achieved 38%–48% less velocity mismatch than the direct simulation. Nevertheless, both the direct simulation and present DA method have uncertainty in the vessel and aneurysm shapes, and thus uncertainty analyses are required to accurately evaluate the flow fields [20]. In this regard, the present DA method has the advantage of focusing only on the aneurysm.
The velocity mismatch between the mathematical model and 4D flow MRI varied among patients. The order of the relative sizes of the mismatches for patients A–C was opposite that of the relative magnitudes of the spatiotemporally averaged velocities obtained from 4D flow MRI; however, the difference in the average velocity values between patients A and C was small compared with the corresponding difference in the velocity mismatches (i.e., A ≅ C > B), whereas the velocity mismatch between patient B and C was small (i.e., A < C ≅ B). When focusing on flow fields, all patients had an organized circulating flow, but several high-velocity vectors were observed in the 4D flow MRI of patients B and C, and some of them appeared to have artifacts near the aneurysm wall. 4D flow MRI is affected by temporal averaging of multiple heartbeats and motion artifacts [21], which result in an incomplete reproduction of actual flow profiles [6,22]. Since DA analysis of patient C using the synthetic data provided an accurate estimation, low-accuracy estimation for the DA analysis using the real-world dataset may be caused by artifacts in the 4D flow MRI. Importantly, low-accuracy velocity estimates in particular patients may result if the surface shape of the aneurysm is not reconstructed well in low-velocity regions using TOF-MRA measurement [23]. The source of the inconsistencies between patients B and C remains unclear, and further investigation is needed.
In 4D flow MRI, the choice of VENC is crucial to adequately quantify flow velocities. The value should be set to cover the maximum flow velocity to avoid aliasing artifacts; however, the S/N ratio is low in the low-velocity region because of phase errors in the 4D flow MRI measurements [24]. The flow velocity inside an aneurysm is generally lower than that in the main vessels, and thus VENC40 could more accurately evaluate the intra-aneurysmal flow velocity than VENC120 because it results in a larger S/N ratio. This explains why DA analysis using VENC40 data resulted in smaller velocity errors than DA analysis using VENC120. Since the proposed DA approach only requires observation of the intra-aneurysmal velocity, it is possible to use a low VENC value in the 4D flow MRI measurement to improve the S/N ratio.
This study had several limitations. First, as mentioned above, reconstruction accuracy of the aneurysm shape limited the estimation accuracy. In this regard, a notable technique that incorporates a shape factor of the vessel wall into an optimization problem as a design variable has been developed [25], and its implementation may improve the estimation accuracy of the present DA method. Also, an approach based on the Darcy equation to describe the vessel wall can be trialed [26]. Second, since this study imposes a point-wise approximation to the kernel, more reasonable approximations to estimate the values may be considered (e.g., volume integral estimation using a reasonable kernel) [12, 27]. Third, although flow fields were considered as a hemodynamic factor, several other hemodynamic factors related to wall shear stress affect the growth and rupture of aneurysms [28]. Incorporation of these indices into the DA model will enhance its applicability. Fourth, although the present DA method has a lower computational cost than other DA approaches, calculations still require nearly 1 d to sufficiently converge the objective function during optimization. More efficient gradient-based optimization techniques, such as the Broyden–Fletcher– Goldfarb–Shanno method [11], will mitigate this issue. Fifth, standardized blood viscosity has a large and unpredictable impact on hemodynamic factors in CFD analyses [29]. Uncertainty analyses or extension of the optimization problem to include blood viscosity will be needed to correctly quantify patient-specific hemodynamic features in DA analyses. Sixth, the aneurysm analysis regions were manually extracted. However, hemodynamic factors around the aneurysm neck play an important role in pathobiological changes, and thus careful extraction procedures are required to accurately reproduce the neck shape. More sophisticated reconstruction techniques [30] will improve the determination of the analysis domain.
In summary, a novel and practical tool to estimate patient-specific intra-aneurysmal flows based on the 4D-Var DA method using 4D flow MRI data was introduced and validated. The proposed DA method may better reproduce intra-aneurysmal flows in vivo than conventional CFD simulations. Moreover, the present DA approach could complete simulations in at most one day on a small-scale parallel computer system. The proposed approach is expected to become a key tool for early, patient-specific prediction of aneurysm growth and rupture.
Data Availability
All data produced in the present study are available upon reasonable request to the authors.
Statements and Declarations
We declare that there are no competing interests to be reported for this publication, and all authors were fully involved in the study and preparation of the manuscript, and that the material has not been and will not be submitted for publication elsewhere.
Acknowledgements
This work was supported by the MEXT Program for Promoting Researches on the Supercomputer Fugaku (Development of human digital twins for cerebral circulation using Fugaku, JPMXP1020230118) and used computational resources of the supercomputer Fugaku provided by the RIKEN Center for Computational Science (project ID: hp230208, hp240220). Some computations were also carried out using the supercomputer “Flow” at the Information Technology Center, Nagoya University. This work was also supported by JSPS KAKENHI Grant Number JP22H00190. The authors thank Prof. Shigeo Wada and Prof. Tomohiro Otani (Osaka University) for fruitful discussion; Mr. Daiki Kurokawa, Mr. Yu Taniguchi (Osaka University), and Mr. Kengo Tanaka (Tokyo Metropolitan University) for preliminary studies; and Edanz (https://jp.edanz.com/ac) for editing a draft of this manuscript.
References
