\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}  {\large\bf}{\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}  {\large\bf}}
\makeatother
\usepackage{secdot}

\makeatletter
\def\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}{\small\bf}{\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}{\small\bf}}
\makeatother

\usepackage{amsmath, amssymb}

\usepackage[dvipdfmx]{graphicx}
\usepackage{float}
\graphicspath{{./}}
\usepackage[labelsep=period, figurename=Fig. ,tablename=Table]{caption}
\usepackage{subcaption}
\renewcommand{\thesubfigure}{\roman{subfigure}}
\renewcommand{\thetable}{\Roman{table}}
\usepackage{booktabs}
\usepackage{multirow}
\usepackage{tabularx}
\newcolumntype{C}{>{\centering\arraybackslash}X}

\usepackage{cite}
\bibliographystyle{jpsj}
\usepackage{url}

\usepackage{ulem}
\usepackage{comment}
\pagestyle{plain}
\usepackage{color}
\usepackage{ascmac}
\usepackage{fancybox}
\usepackage{seqsplit}

\def\vector#1{\mbox{\boldmath $#1$}}
\newcommand{\tc}{\textcolor{red}}

\title{Does soccer players' sprint follow the Fujimura-Sugihara model?: validation using soccer tracking data}

\usepackage{authblk}
\author[1]{Takuma Narizuka\thanks{{\it E-mail address}: pararel@gmail.com (T. Narizuka).\\
}}
\author[2]{Kenta Takizawa}
\author[3]{Yoshihiro Yamazaki}
\affil[1]{Faculty of Data Science, Rissho University, Kumagaya, Saitama 360-0194, Japan}
\affil[2]{Department of Physics, Faculty of Science and Engineering, Chuo University, Bunkyo, Tokyo 112-8551, Japan}
\affil[3]{Department of Physics, School of Advanced Science and Engineering, Waseda University, Shinjuku, Tokyo 169-8555, Japan}
\date{}

\providecommand{\keywords}[1]
{
  \small
  \textbf{\textit{Keywords---}} #1
}

\begin{document}
	\maketitle
	\baselineskip 20pt

\begin{abstract}
\baselineskip 20pt
In soccer game analysis, the widespread availability of play-by-play and tracking data has made it possible to test mathematical models that have been discussed mainly theoretically.
One of the essential models in soccer game analysis is a motion model that predicts the arrival point of a player in $ t $ s.
Although many space evaluation and pass prediction methods rely on motion models, the validity of each has not been fully clarified.
This study focuses on the motion model proposed by Fujimura and Sugihara (Fujimura-Sugihara model) under sprint conditions based on the equation of motion.
A previous study indicated that the Fujimura-Sugihara model is ineffective for soccer games because it generates a circular arrival region.
In this study, we quantitatively examined the validity and limitations of the Fujimura-Sugihara model using soccer tracking data.
We show that the boundary of the player's arrival region is circular rather than elliptical.
Furthermore, we show that the initial speed dependence of the arrival region satisfies the solution of the Fujimura-Sugihara model.
We also proposed a method to estimate kinetic parameters directly from tracking data and discussed the limitations of the Fujimura-Sugihara model for soccer games.
\end{abstract}

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}  {\large\bf}{Introduction}
Over the past decade, various types of data have become available for sports data analysis.
In particular, the prevalence of play-by-play \cite{Pappalardo2019} and tracking data \cite{Pettersen2014} for soccer games has enabled new analyses that were previously impossible\cite{Gudmundsson2017}.
Some examples include ball-passing network analysis \cite{Duch2010, Buldu2019}, formation analysis \cite{Bialkowski2014, Narizuka2019}, and space evaluation \cite{Fernandez2018, Spearman2018, Narizuka2021}.
The widespread use of these new data has also made it possible to test mathematical models that have mainly been discussed theoretically.
This study aimed to validate a motion model based on soccer tracking data.

The motion model is an essential mathematical model in soccer games, which calculates the arrival point of a player in $ t $ s based on the current location and velocity.
The two practical applications of motion models in soccer games are space evaluation and pass prediction.
In space evaluation, a fundamental concept is the ``dominant region,'' defined as the region where a specific player can arrive before other players \cite{Taki2000, Gudmundsson2017}.
In general, we can estimate the dominant region of each player by comparing the arrival times of all players to each location in the field.
In addition, the outcome of a given pass can be estimated by calculating the arrival time of each player on the ball trajectory \cite{Fujimura2005, Spearman2017, Alguacil2020, Anzer2022}.
Because the motion model can calculate the arrival times of a player to a specific location, it is essential for soccer game analysis.

Thus far, various motion models have been proposed.
The simplest motion model assumes uniform linear motion for all players, resulting in the Voronoi region as the dominant region.
More realistic models have also been proposed; for example, Taki and Hasegawa assumed uniform accelerated motion of players \cite{Taki1996, Taki2000}, and Fujimura and Sugihara considered viscous resistance \cite{Fujimura2005}.
These motion models are based on the equation of motion and are often referred to as the ``physics-based motion model.''
They provide the arrival points of players under specific conditions such that each player moves to all locations by sprinting.

Another type of motion model calculates the arrival points of players based on machine learning \cite{Gudmundsson2014, Brefeld2019, Caetano2021}.
This ``probabilistic motion model'' can predict realistic arrival points based on previous tracking data, though it is costly for learning.
We note that because the tracking data used for learning include various running patterns other than a sprint, the predicted arrival point does not necessarily mean the position arrived by sprinting.

Both physics-based motion models and probabilistic motion models aim to predict the arrival point of players.
Physics-based motion models also play a role in elucidating the principles of players' movement laws.
This study focuses on the motion model proposed by Fujimura and Sugihara (hereinafter Fujimura-Sugihara model) in sprint conditions to elucidate soccer players' movement laws during sprinting.
Specifically, we investigated the validity and limitations of the Fujimura-Sugihara model based on soccer tracking data.
We stress that the motion model in the sprint condition is significant in situations where we evaluate which player can receive the ball.

In our investigation, we first focused on the shape of the arrival region of players predicted by the motion model.
The Fujimura-Sugihara model generates a circular arrival region \cite{Fujimura2005}.
However, it has been pointed out that the arrival region should be elliptical considering the difference in the acceleration ability of players depending on the direction of motion \cite{Brefeld2019, Fernandez2018}.
We show that the circular arrival region is obtained from soccer tracking data, and the arrival region's initial speed dependence also satisfies the Fujimura-Sugihara model's solution.
Next, we propose a method to estimate the kinetic parameters of the Fujimura-Sugihara model.
Contrary to previous experiment-based methods, this method can estimate valid parameters directly from tracking data \cite{Fujimura2005}.
Finally, we discuss the limitations of the Fujimura-Sugihara model for soccer games based on the time dependence of the estimated kinetic parameters.

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}  {\large\bf}{Method}

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}{\small\bf}{Data}
The datasets used in this study were from 54 soccer games in the top division of the Japan Professional Football League (J1 League).
Each game was played by 18 teams in 2016.
The primary data of our dataset is absolute positional coordinates $ (x, y) $ of all players every 0.04 s, collected by using the TRACAB system \cite{tracab}.
The $ x $ and $ y $ coordinates are considered to contain an error of $ \pm 1 $ m by assessing the accuracy of the TRACAB system.
These datasets were provided by DataStadium Inc., Japan, which was authorized to collect and sell these data under a contract with the J League \cite{DataStadium}.
This contract also ensures that the use of relevant datasets does not infringe on the rights of players and clubs belonging to the J League.
Although the datasets were proprietary, we received explicit permission from DataStadium Inc. for use in this study.
This study's Data analyses and visualizations were performed using Python packages on an iMac Pro system with a 3-GHz 10-Core Intel Xeon W processor and 128 GB of memory.

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}{\small\bf}{Fujimura-Sugihara model}
We summarize the motion model proposed by Fujimura and Sugihara \cite{Fujimura2005}.
For the position $ \vec{x}(t) $ of a player at time $ t $, the Fujimura-Sugihara model is given by the following equation of motion:
\begin{align}
	m \frac{d^{2} \vec{x}(t)}{d t^{2}} &= F \vec{n} - k \frac{d \vec{x}(t)}{d t},
	\label{eq:fujimura}
\end{align}
where $ m $ is the mass of the player, $ F $ and $ \vec{n} $ are the magnitude and direction of the driving force, respectively, and $ k $ is the coefficient of viscous resistance.
In other words, the player accelerates in the direction of $ \vec{n} $ with a magnitude $ F $ and it becomes harder to accelerate in proportion to its velocity $ d\vec{x}(t)/dt $.
Given an initial position $ \vec{x}_{0}=(0, 0) $ and an initial velocity $ \vec{v}_{0} $, the solution is given as follows:
\begin{align}
	\vec{x}(t) &= \frac{1 - \exp(-\alpha t)}{\alpha} \vec{v}_{0}
	+ V_{\mathrm{max}} \left(t - \frac{1 - \exp(-\alpha t)}{\alpha}\right)\vec{n},
	\label{eq:solution}
\end{align}
where $ V_{\mathrm{max}} = F/k $ and $ \alpha = k/m $ are arbitrary constants referred to as ``kinetic parameters.''
Here, we define the first and second coefficients in Eq. \eqref{eq:solution} as
\begin{align}
	A(\alpha, t) &= \frac{1 - \exp(-\alpha t)}{\alpha}, \label{eq:A}\\[10pt]
	B(\alpha, V_{\mathrm{max}}, t) &= V_{\mathrm{max}} \left(t - \frac{1 - \exp(-\alpha t)}{\alpha}\right) = V_{\mathrm{max}} (t - A(\alpha, t)). \label{eq:B}
\end{align}
When the direction of the driving force $ \vec{n} $ changes arbitrary, the arrival points in $ t $ [s] of the player are distributed on a circle with center $ A(\alpha, t) \vec{v}_{0} $ and radius $ B(\alpha, V_{\mathrm{max}}, t) $; this circle is referred to as ``arrival circle'' hereinafter.
Fujimura and Sugihara conducted a sprint experiment with three amateur hockey players to estimate the kinetic parameters.
By fitting the obtained speed curve with the solution of the Fujimura-Sugihara model, they empirically obtained the kinematic parameters $ V_{\mathrm{max}}=7.8 $ m/s and $ \alpha=1.3 $ 1/s as the typical values during sprinting \cite{Fujimura2005}.

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}{\small\bf}{Investigation method}
To characterize the arrival points of players in $ \Delta t $ [s] using tracking data, we adopt the method used in previous studies \cite{Fujimura2005, Brefeld2019}.
First, the velocity $ \vec{v}(t) $ of each player at time $ t $ is converted to start at the origin $ (0, 0) $ and point in the positive direction of the $ x $ axis (Fig. \ref{coord}).
We then plot the location of the same player at $ t+\Delta t $ on the same coordinates.
By repeating this plot for various players and $ t $, we obtain a heat map formed by the arrival points in $ \Delta t $ [s].

The shape of the heat map reflects the movement pattern of each player during $ \Delta t $.
For example, when players lose speed during $ \Delta t $, they arrive at a point within the heat map.
However, they arrived at a point around the boundary in the case of sprinting.
Because we investigated the validity of the Fujimura-Sugihara model under sprint conditions, we focused on the shape of the heat map's boundary and compared it with the solution \eqref{eq:solution}.
We set $ \vec{v}_{0} = (v_{0}, 0) $, where $ v_{0} > 0 $, in Eq. \eqref{eq:solution} for comparison with the heat map.

In the following analyses, we calculated heat maps using all data during the playing time of 54 games for all players except the goalkeepers.
The controlling parameters were the initial speed $ v_{0} $ [m/s] and time interval $ \Delta t $ [s].
Specifically, for a fixed $ \Delta t $, we calculated the heat maps for the initial speed $ [v_{0}, v_{0} + \Delta v_{0}) $, where $ \Delta v_{0} = 0.3 $ m/s ($ \simeq 1 $ km/h).
We set one cell side in the heat map to $ 0.2\times \Delta t $ [m], and only cells adjacent to more than $ c $ nonblank cells were used for the analyses.
To exclude isolated cells as outliers, we manually set the value of $ c $ to eight for $ \Delta t < 1 $, six for $ 1 \leq \Delta t \leq 3 $, and four for $ 3 < \Delta t $.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=6cm]{./figure/coord.pdf}
	\caption{Coordinate system for the calculation of heat maps.}
	\label{coord}
\end{figure}

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}  {\large\bf}{Result}

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}{\small\bf}{$ v_{0} $ dependence of heat map}
Figures \ref{heatmap}(a) and \ref{heatmap}(b) present the heat maps of the players' arrival points for various initial speeds, $ v_{0}=1,\ 3,\ 5,\ 7 $ m/s, where $ \Delta t = 1 $ and $ 2 $ s.
We focused only on the shape of each heat map's boundary, although the color gradation of the heat map is proportional to the number of data points.
Remarkably, Figs. \ref{heatmap}(a) and (b) show that the boundary of the heat map is not elliptical but circular for any initial $ v_{0} $.
This result is consistent with the solution \eqref{eq:solution} of the Fujimura-Sugihara model.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=16cm]{./figure/hmap.pdf} \\
	\caption{Heat maps of players' arrival points for various initial speed $ v_{0}=1, 3, 5, 7 $ m/s, where (a) $ \Delta t = 1 $ s and (b) $ \Delta t = 2 $ s. The dotted line in each panel represents the arrival circle \eqref{eq:solution} with kinetic parameters estimated from the heat maps.}
	\label{heatmap}
\end{figure}

Next, we approximated the boundary of each heat map as a circle and estimated its center coordinates $ (x_{\mathrm{c}}, y_{\mathrm{c}}) $ and radius $ r_{\mathrm{c}} $.
After excluding the isolated heatmap cells as outliers, we calculated the maximum and minimum values of the heat maps $ x $ and $ y $ coordinates, $ x_{\mathrm{max}} $, $ x_{\mathrm{min}}$, $ y_{\mathrm{max}} $, and $ y_{\mathrm{min}} $.
Then, $ (x_{\mathrm{c}}, y_{\mathrm{c}}) $ and $ r_{\mathrm{c}} $ were estimated as follows:
\begin{align*}
	(x_{\mathrm{c}}, y_{\mathrm{c}}) &= \left(\frac{x_{\mathrm{max}}+x_{\mathrm{min}}}{2}, \frac{y_{\mathrm{max}}+y_{\mathrm{min}}}{2}\right), \\
	r_{\mathrm{c}} &= \frac{x_{\mathrm{max}} - x_{\mathrm{min}} + y_{\mathrm{max}}-y_{\mathrm{min}}}{4}.
\end{align*}
Figure \ref{xyc-v0} shows the $ v_{0} $ dependence of $ x_{\mathrm{c}} $ and $ y_{\mathrm{c}} $.
We find that $ y_{\mathrm{c}} $ is independent of $ v_{0} $; namely, $ y_{\mathrm{c}} \simeq 0 $.
However, $ x_{\mathrm{c}} $ is proportional to $ v_{0} $, particularly for $ v_{0} \lesssim 6 $ m/s.
The dotted lines in each panel of Fig. \ref{xyc-v0} represents the regression line for the points where $ v_{0} \leq 6 $ m/s.
These results are consistent with the solution \eqref{eq:solution} of the Fujimura-Sugihara model.

We also examined the $ v_{0} $ dependence of the radius $ r_{\mathrm{c}} $ of the estimated circle.
Figure \ref{rc-v0} shows that $ r_{\mathrm{c}} $ becomes almost constant, particularly for $ v_{0} \lesssim 6 $ m/s.
The dotted lines in each panel of Fig. \ref{rc-v0} represents the average value of $ r_{\mathrm{c}} $ where $ v_{0} \leq 6 $ m/s.
As $ B(\alpha, V_{\mathrm{max}}, t) $ in Eq. \eqref{eq:solution} is independent of $ v_{0} $, this result is also consistent with the solution \eqref{eq:solution} of Fujimura-Sugihara model.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=14cm]{./figure/xyc-v0.pdf}
	\caption{$ v_{0} $ dependence of the center coordinates $ x_{\mathrm{c}} $ and $ y_{\mathrm{c}} $ of the estimated circle of each heat map, where (a) $ \Delta t = 1 $ s and (b) $ \Delta t = 2 $ s. The dotted line in each panel represents the regression line for points where $ v_{0} \leq 6 $ m/s.}
	\label{xyc-v0}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=14cm]{./figure/rc-v0.pdf}
	\caption{$ v_{0} $ dependence of the radius $ r_{\mathrm{c}} $ of the estimated circle of each heat map. The dotted line represents $ r_{\mathrm{c}} = 6.16 $ m for (a) $ \Delta t = 1 $ s, and $ r_{\mathrm{c}} = 12.94 $ m for (b) $ \Delta t = 2 $ s. These values are obtained by averaging the points where $ v_{0} \leq 6 $ m/s.}
	\label{rc-v0}
\end{figure}

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}{\small\bf}{Estimation of kinetic parameters using heat maps}
According to solution \eqref{eq:solution}, the $ x $ coordinate of the arrival circle is proportional to $ v_{0} $, and the proportionality coefficient is given by Eq. \eqref{eq:A}.
Because we obtained the result that $ x_{\mathrm{c}} $ calculated from the heat map is proportional to $ v_{0} $ (refer to Fig. \ref{xyc-v0}), we can estimate the kinetic parameter $ \alpha $ using Eq. \eqref{eq:A}.
Specifically, the proportionality coefficients obtained from the regression line shown in Fig. \ref{xyc-v0}(a) and (b) are 0.58 and 0.84; then, we obtain $ \alpha=1.23 $ and $ 1.04 $ 1/s for $ \Delta t = 1 $ and $ 2 $ s, respectively.

The radius of the arrival circle is given by Eq. \eqref{eq:B}.
As shown in Fig. \ref{rc-v0}(a) and (b), the radius of the estimated circle of each heat map has become $ r_{\mathrm{c}}\simeq 6.16 $ and 12.94 m for $ \Delta t = 1 $ and $ 2 $ s.
Thus, the estimated values of $ \alpha $ and Eq. \eqref{eq:B} can yield another kinetic parameter $ V_{\mathrm{max}} $; we obtain $ V_{\mathrm{max}}=14.53 $ and $ 11.19 $ m/s, respectively.

The dotted lines in each panel of Fig. \ref{heatmap} shows the arrival circle calculated by the solution \eqref{eq:solution} with the above estimated kinetic parameters $ \alpha $ and $ V_{\mathrm{max}} $.
We found that the solution \eqref{eq:solution} of the Fujimura-Sugihara model can correctly predict the boundary of the arrival region.

\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}{\small\bf}{$ \Delta t $ dependence of kinetic parameters}
We also analyze the dependence of the kinetic parameters on $ \Delta t $.
As a result of the same analyses in the previous section for different $ \Delta t $, we confirm the characteristics shown in Fig. \ref{heatmap} to Fig. \ref{rc-v0} for all $ \Delta t $.
We present the $ \Delta t $ dependence of the estimated kinetic parameters in Fig. \ref{delta_t}.
Although $ \alpha $ and $ V_{\mathrm{max}} $ are assumed to be constant in the Fujimura-Sugihara model, we found that the estimated kinetic parameters vary with $ \Delta t $, particularly in the range of $ \Delta t \lesssim 1 $ s.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=8cm]{./figure/kp-dt.pdf}
	\caption{$ \Delta t $ dependence of kinetic parameters, $ \alpha $ and $ V_{\mathrm{max}} $.}
	\label{delta_t}
\end{figure}

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}  {\large\bf}{Discussion and Conclusion}
We investigated the validity of the Fujimura-Sugihara model based on heat maps of the players' arrival points obtained from soccer tracking data.
Our results can be summarized as follows:
First, the boundary of the heat map became a circle rather than an ellipse.
Second, $ x_{\mathrm{c}} $ was proportional to $ v_{0} $ and $ y_{\mathrm{c}} $ was independent of $ v_{0} $.
Third, $ r_{\mathrm{c}} $ was independent of $ v_{0} $.
These results are consistent with the solution \eqref{eq:solution} of the Fujimura-Sugihara model.
In addition, we proposed a method for estimating the kinetic parameters included in the Fujimura-Sugihara model.
This estimation method is more reasonable than the previous method based on experiments \cite{Fujimura2005} in the sense that the parameters are estimated directly from the tracking data.
Meanwhile, the estimated kinetic parameters varied with $ \Delta t $, particularly in $ \Delta t \lesssim 1 $ s; this result is inconsistent with the assumption of the Fujimura-Sugihara model.

Note that as shown in Fig. \ref{heatmap} to Fig. \ref{rc-v0}, the plot in large $ v_{0} $ region is different from the others.
This discrepancy is because of the few data points of the heat maps in the large $ v_{0} $ region.
For example, a player with large $ v_{0} $ is unlikely to move in the opposite direction during $ \Delta t $.
Namely, the larger the $ v_{0} $, the fewer the data points in the region $ x < 0 $ in the heat maps.
The estimated values of $ (x_{c}, y_{c}) $ and $ r_{c} $ for large $ v_{0} $ values appear to be incorrect because they are calculated using the maximum and minimum values of the heat maps $ x $ and $ y $ coordinates.
In the future, the use of large-scale tracking data could eliminate such problems.

We comment on the results of the $ \Delta t $ dependence of the kinetic parameters shown in Fig \ref{delta_t}.
This result requires careful consideration from two perspectives.
First, the $ x $ and $ y $ coordinates of our tracking data contain an error of $ \pm 1 $ m by assessing the accuracy of the TRACAB system.
This error cannot be ignored in the analysis of the heat maps, especially when $ \Delta t \lesssim 1 $ s.
Therefore, we must check the validity of the rapid increase in $ \alpha $ and $ V_{\mathrm{max}} $ in Fig. \ref{delta_t} using more accurate data.
The second possibility is that the results in Fig. \ref{delta_t} indicate the limitation of the Fujimura-Sugihara model.
The Fujimura-Sugihara model must be extended to reproduce the behavior shown in Fig. \ref{delta_t}.
One possible extension is a model in which the magnitudes of the driving force $ F $ and viscous resistance $ k $ in Eq. \eqref{eq:fujimura} are time-dependent:
\begin{align}
	m \frac{d^{2} \vec{x}(t)}{d t^{2}} &= F(t) \vec{n} - k(t) \frac{d \vec{x}(t)}{d t}.
	\label{eq:fujimura2}
\end{align}
This equation is known as the variable coefficient second-order linear ODE.
As the time dependence of kinetic parameters is conceivable in soccer, the analysis of this extended motion model can be a challenging and significant future topic in soccer game analysis.

Previously, Brefeld et al. demonstrated that the arrival region of soccer players becomes elliptical \cite{Brefeld2019}.
Fern\'{a}ndez et al. modeled the player influence area as elliptical when the initial speed is large \cite{Fernandez2018}.
However, Anzer et al. pointed out that the smaller the speed interval $ \Delta v_{0} $, the closer the shape of the arrival region is to being circular rather than elliptical \cite{Anzer2022}.
This study supports the results of Anzer et al., with a more detailed and comprehensive analysis.
Furthermore, we have recently shown that soccer players' sprints satisfy the characteristics of the Fujimura-Sugihara model, that is, $ v_{0} $ dependence of the arrival circle for fixed $ \Delta t $.
Our results suggest that a relatively simple model can describe soccer players' sprints, although a slight discrepancy exists between actual observations and model predictions.
Furthermore, the motion model under the sprint condition was utilized for pass prediction \cite{Fujimura2005, Spearman2017, Alguacil2020} and for modeling the dominance and influence of soccer players at each location in the field \cite{Narizuka2021}.
Our observations are fundamental characteristics that the motion model should satisfy and provide direction for modeling soccer players' motions.

\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus -.2ex}{2.3 ex plus .2ex}  {\large\bf}*{Acknowledgements}
The authors are grateful to DataStadium Inc., Japan, for providing the player tracking data for this study.
This work was partially supported by the Data-Centric Science Research Commons Project of the Research Organization of Information and Systems, Japan, and by a Grant-in-Aid for Young Scientists (18K18013) from the Japan Society for the Promotion of Science (JSPS).


\section{Introduction}
Coarctation of the aorta (CoA) refers to a local narrowing of the aortic arch. It makes up 6 - 8$\%$ of all congenital heart diseases \cite{Rafieianzab2021}, and is often associated with other cardiovascular diseases, such as aortic arch hypoplasia, subaortic stenosis, ventricular and atrial septal defects \cite{Hoffman2002,Reller2008,Alkashkari2019}. The coarctation leads to high blood pressure and thus heart damage. Stent intervention, which is generally performed based on clinical experience, is a recommended therapy to reduce the pressure gradient and restore blood flow.

With increasing computational power, in silico modeling is emerging as a
promising tool to help clinicians with intervention planning and to evaluate the outcome of therapies, such as stenting for intracranial aneurysm \cite{Zhong2016,Berg2018}, abdominal aortic aneurysm \cite{ Pionteck2021}, and type-B aortic dissection \cite{Chen2018,Kan2021}. By taking personalized information as input, modeling also supports the design of patient-specific medical implants \cite{McCloy2001}.

For in silicio modeling of personalized stent intervention in CoA, a protocol for virtual geometry deformation \cite{Neugebauer2016} and a validated numerical method to accurately predict the blood flow in the aorta are required. Regardless of the erythrocytes, leukocytes, and platelets in blood, the flow in the aorta is normally modelled as Newtonian fluid \cite{pedley1980} considering the relatively large Reynolds number $Re$, which is proportional to the flow velocity and aorta diameter and inversely proportional to the blood viscosity. Computational fluid dynamics (CFD) plays an important role in resolving hemodynamics \cite{Taha2005}. Due to the personalized and complex 3D geometry and jet flows induced by heart contraction and local narrowing, laminar flow, turbulent flow and transition between them may coexist spatiotemporally \cite{Stein1976, Ku1997}. Thus, to accurately resolve such aortic flow, both turbulence and complex geometry should be considered in CFD simulations. Three approaches, including Reynolds-averaged Navier–Stokes equations (RANS), large eddy simulation (LES), and direct numerical simulation (DNS), are typically used for turbulence modeling. From RANS to DNS, both the accuracy and computational demand increases due to more and more degrees of freedoms that need to be resolved.

So far, mainly RANS and LES were used to study aortic flow in literature \cite{Caballero2013}. With a transitional model, RANS was adopted to resolve flows in patient-specific thoracic aortic aneurysm by Tan et al. \cite{Tan2009} and aortic dissection by Cheng et al. \cite{Cheng2013,Cheng2015}. Simulations were carried out using ANSYS CFX, a commercial finite volume-based solver. Kouseral et al. studied flow stability in a normal aorta using the same numerical method, and compared their numerical results with experimental data from in vivo magnetic resonance imaging (MRI) \cite{Kousera2012}. They concluded that the RANS based shear stress transport transitional model was capable of capturing the correct flow state when low inflow turbulence intensity (1.0\%) was specified. Miyazaki et al. validated three CFD models for aortic flows in the aorta of a healthy adult and a child with double aortic arch \cite{Miyazaki2017}. Laminar, LES and the renormalization group (RNG) \textit{k}-$\varepsilon$ model were considered and compared. Simulations were performed using another finite volume-based solver, ANSYS Fluent. Their results show that the RNG \textit{k}-$\varepsilon$ model has the highest correlation with data from 4D flow MRI. Recently, Manchester et al. used LES to study the blood flow in patient-specific aorta with aortic valve stenosis \cite{Manchester2021}. Here, the finite volume based open-source library OpenFOAM, was used. After investigating the fluctuating kinetic energy, wall shear stress (WSS) and energy loss, they concluded that turbulence played an important role in aortic hydrodynamics.

It should be noted that severe turbulence will be encountered in CoA due to a more complex geometry and larger $Re$, which might lead to higher requirements on the CFD method. The aforementioned conventional CFD methods are based on discretizations of macroscopic governing equations, such as the Navier-Stokes (NS) equations. Alternatively, the lattice Boltzmann method (LBM) is based on the mesoscopic Boltzmann equation and has multiple advantages, including simply handling of complex geometric shapes, ease of programming, and suitability for parallelization \cite{kruger2017lattice, he2009lattice, chen1998lattice}. Therefore, the LBM is increasingly used for the simulation of turbulent flow \cite{Chen2003} and biological fluid flows \cite{Wang2014}. Hennt et al. simulated the unsteady blood flow in a patient-specific geometry with a moderate thoracic aortic coarctation, and demonstrated that the LBM based DNS was capable of resolving such complex flow \cite{Henn2013}. Recently, Mirzaee et al. studied aortic flows for 12 patients with CoA using the LBM based LES, particularly with the Smagorinsky turbulence model \cite{Mirzaee2017}. A reasonable agreement for pressure drop between the numerical results and the catheter measurements was achieved. Nevertheless, to guide in silico stent intervention for CoA, a comprehensive validation for the LBM based LES for complex flow is still missing.

Since the 1970s and 1980s, MRI has become an important clinical and scientific tool that is widely used for diagnosis, monitoring of treatment procedures, and for biomedical research \cite{Markl2012}. Compared with X-ray and computed tomography, one of the advantages of MRI is the use of non-ionizing radiation \cite{McRobbie2007, Landheer2020}. In addition to obtaining anatomical information, MRI can also be used for quantitative flow measurements using phase-contrast imaging \cite{Moran1982, DonnellM1985,Soulat2020} including measurement of aortic blood flow \cite{Miyazaki2017,Saitta2019}.

In this study, the LBM based LES and DNS were adopted to resolve blood flow in
the human aorta. Geometries for a patient-specific aorta with CoA, before and
after stent intervention, were considered and physical phantoms were created
using 3D printing for use in MRI flow experiments. Flow measurements obtained
with MRI scans were then used as boundary conditions for simulations. Obtained
numerical results using LES and DNS were then compared with experimental
4D-Flow data. To further validate the LBM based LES, we also compared with
in vivo data. We demonstrated that LES is capable of accurately simulating
complex aortic flow and further applied it for in silico stent implantation.
Details of the methodology are given in Section 2. Numerical results and
experimental measurements for aortic flows are presented and comapred in Section
3. The application for stent selection is provided in Section 4.
Discussion and conclusions can be found in Section 5.

\section{Methodology}

\subsection{MRI Experiments}

The anatomical structure of the aorta and the flows therein were acquired by
MRI, which provided realistic geometries and boundary conditions for CFD
simulations. Comparison between CFD and phase-contrast flow MRI for flows in
3D-printed phantoms and in vivo aorta were performed respectively.

For the phantoms, the 3D anatomies of the heart and aorta of a 14-year-old
patient with CoA, before and after stent intervention, were reconstructed from
images obtained by a Magnetom Skyra 3T (Siemens Healthineers, Erlangen,
Germany). The stent (diameter = 12 mm) used in this patient was a covered
Cheatham-platinum (CP) stent made of platinum-iridium (NuMed, Orlando, USA).
The sequence parameters
are listed in Table 1. Using ITK-SNAP \cite{YUSHKEVICH20061116}, the geometry
that starts from the aortic root and ends above the diaphragm was segmented
based on the grey values and exported as STL file. The main branches, such as
the right subclavian, the left subclavian, the right carotid artery and the
left carotid artery, were included. To have an uniform surface mesh, the
generated geometries were then remeshed using Autodesk Meshmixer
\cite{SchmidtR2010}. The schematic diagram of the experiment is shown in Fig.
\ref{fig:FlowPhatom}. Two aortic models, including the pre-interventional and
post-interventional geometries, were printed using the Stratasys’ high-end 3D
laser printer Connex 3 using biocompatible MED610 as material. The phantoms
were connected to a pump. Forced water flows therein were then measured using
4D phase-contrast Flow MRI \cite{Markl2012, Dyverfeldt2015, Potters2014}.
The sequence parameters can be found in Table 1. Every case was measured three
times with about 30 mins per measurement. The averaged flow fields were then
used for comparison.

For the in vivo validation, the aortic blood flow of a 3-year-old patient,
was obtained using a 2D flow sequence instead of a 4D flow sequence, to reduce the
duration of measurement. In the 2D measurement, the through-plane velocity of
the flow  was measured in two planes located in the ascending aorta and
descending aorta respectively. All in vivo measurements were made with the use
of ECG triggering and respiratory gating. More details can be found in Table 2.
For further CFD simulation, the aortic geometry was segmented and reconstructed
using the same procedure as mentioned above.

\begin{table*}
\small
\caption{MRI sequences used for the flow-phantom study:
In vivo MRI scans were performed pre and post intervention
to obtain the anatomical structure. These geometries were
used for 3D-printing a phantom then used
for 4D-flow MRI measurements.}

\begin{center}
\begin{tabular}{c|cccc}
\hline\hline
          & Geometry (pre) & Geometry (post) & 4D flow MRI (pre) & 4D flow MRI (post) \\
\hline
Sequence type &{3D FLASH (TWIST)} &{3D T1 weighted FLASH} &3D Cartesian FLASH & 3D Cartesian FLASH \\
Acceleration & 3 $\times$ 2 & 2 $\times$ 2 & - & - \\
Matrix size &352 $\times$ 246 &448$\times$252  &384 $\times$504  &416$\times$364\\
Number of slices  &80 &88 &144 &144\\
Slice thickness (mm)  &1.30 &1.20 &0.77 &0.77\\
Pixel size (mm$^2$) &1.02$\times$1.02 &0.89$\times$0.89 &0.77$\times$0.77 &0.77$\times$0.77\\
Repetition time (ms)  &2.75 &3.70 &36.40 &70.40\\
Echo time (ms)  &1.00 &1.31 &4.61 &7.46 \\
Flip angle ($^\circ$)  &20 &25 &7 &7\\
Velocity encoding (cm/s)  &- &- &50 &40\\

\hline\hline
\end{tabular}
\label{tab:sequences}
\end{center}
\end{table*}

\begin{table*}
\small
\caption{MRI sequences for geometry and 2D flow
for in vivo validation.}
\begin{center}
\begin{tabular}{c|cc}
\hline\hline
          & Geometry & 2D flow MRI\\
\hline
Sequence type &{3D T1 weighted FLASH} &{2D T1 weighted FLASH}\\
Acceleration & 3$\times$2 & 2 $\times$ 2  \\
Matrix size &320$\times$260  & 192$\times$119\\
Number of slices  &88 &30 \\
Slice thickness (mm)  &1.00 &5.00\\
Pixel size (mm$^2$) & 1.00$\times$1.00 & 1.56$\times$1.56\\
Repetition time (ms)  & 312.01&39.44\\
Echo time (ms)  &1.64 &2.67\\
Flip angle ($^\circ$)  &20 &20\\
Velocity encoding (cm/s)  &- &200\\

\hline\hline
\end{tabular}
\label{tab:sequences}
\end{center}
\end{table*}

\begin{figure}
\centering
\center{\includegraphics[trim={0 0 0 0}, clip, width=0.9\linewidth]{FlowPhantom.pdf}}
\caption{Schematic diagram of the phantom measurement. Left: anatomical 3D
structures of a heart and great vessels (with CoA) were segmented and
reconstructed based on MRI images to obtain the geometry of the aorta.
Middle: 4D-Flow MRI of 3D printed phantoms
connected to a pump were performd for pre-interventional  and
post-interventional geometries. Right: 4D-Flow MRI yields three
components (U, V, W) of the instantaneous flow velocity at each voxel.}
\label{fig:FlowPhatom}
\end{figure}

\subsection{Numerical Modeling}

The CFD simulations in this study were performed using the LBM, which is based
on the kinetic theory, particularly the Boltzmann equation which describes the
movements of fluid particles \cite{kruger2017lattice, he2009lattice,
chen1998lattice,He1997}. For simulations, the space and time are discretized into
finite nodes and time steps. Starting from an initial state, the configuration
of the fluid particles at each time step evolves in two sub-steps, streaming
and collision. During streaming, fluid particles at a node move to the
neighbouring nodes along specified discrete directions as defined by
the lattice. The streamed particles at a node collide with each other and
change their velocity distribution functions \cite{Roberto1992}. For 3D flows,
the most popular lattice is the D3Q19, which is used in this work.

Different operators, such as the single-relaxation time BGK \cite{Qian1992} and
the multi-relaxation-time (MRT) operators \cite{DHumieres2002}, can be used to
approximate the particle collision. We chose the MRT operator due to its better
numerical stability. The governing equation for the LBM with MRT operator reads
\begin{equation}
f_{i}\left(\mathbf{x}+\mathbf{e}_{i} \Delta t, t+\Delta t\right)-f_{i}(\mathbf{x}, t)=\Lambda_{ij} ( f_j^{eq}(\mathbf{x}, t)-f_j(\mathbf{x}, t)),
\label{eq:lbm}
\end{equation}
in which $f_{i}$ is the particle velocity distribution function along the $i$th
direction; $\mathbf{x}$ and $t$ are the spatial coordinate and time
respectively; $\Delta t$ is time step; $\mathbf{e}_{i}$ is the discrete
velocity of the lattice along the $i$th direction. The right-hand side of Eq.
(\ref{eq:lbm}) represents the collision process in momentum space.
$\Lambda_{ij}=\mathbf{M}^{-1}\mathbf{S}\mathbf{M}$; $\mathbf{M}$ is a given
transformation matrix for the lattice; $\mathbf{S}$ is a diagonal matrix.
Macroscopic parameters, such as the fluid density, pressure and velocity, are
moments of $f_i$.

The left-hand side and right-hand side of Eq. (\ref{eq:lbm}) represent the
streaming and the collision processes respectively. The simplicity of this
equation implies that the LBM is readily parallelizable as the non-local
streaming is linear while the non-linear collision is local
\cite{kruger2017lattice, he2009lattice, chen1998lattice}. Thus, the LBM is
increasingly used for turbulence modeling, especially DNS with high performance
modern computers. Additionally, due to its particle feature, even with a simple
Cartesian grid the LBM can resolve flow with complex geometry, such as the
patient-specific aortas considered in this study.

The LBM based DNS and LES were investigated in this study, based on the
open-source library Palabos \cite{LATT2021334}. DNS resolves the flow at all
scales without empirical model, thus is regarded as a kind of numerical
experiment. As mentioned before, the calculation cost of DNS is very high
especially for flows with large Reynolds number \cite{Moin1998,Wang2014}.
Alternatively, LES explicitly solves large eddy current and implicitly
calculates small eddies by using a sub-grid scale (SGS) model, thus balancing
accuracy and computational cost. The Smagorinsky SGS model \cite{Joseph1963}
was incorporated into the LBM in this study.

For all simulations, the inlet velocity with Poiseuille profile was specified
at the ascending aorta. The flow rates were given based on MRI measurements and
can be found in the following sections. Outlet boundary condition with a
reference pressure was applied to the descending aorta. The curved aortic wall
was assumed to be no-slip and treated with an extrapolation scheme
\cite{Guo2002}.

\section{Validation and Comparison}

\subsection{Phantom Experiments}

The two phantoms filled with water were used in the MRI
experiments and compared to CFD simulations of the same geometries.
The main branch blood vessels were closed, to
reduce their influence. Inlet and outlet of the geometries were extended
artificially for the connection of the water pipe. Water is incompressible
and Newtonian. Its density and kinematic viscosity are 1.0 kg/m$^{3}$ and 1.0
$\times 10^{-6} \text{m}^2\text{/s}$, respectively. The averaged velocity at
the inlet is 0.1 m/s.  Different spatial resolutions were considered for every
single case. After mesh independence tests, 12.30 million (DNS) and 3.45
million (LES) lattice nodes were selected for the pre-interventional geometry
and 5.12 million (DNS) and 1.25 million (LES) lattice nodes for the
post-interventional one.

Due to the complex geometry, \textit{e.g.} multiple plane curvatures and
branches, blood flow in the patient-specific aorta is unsteady and complicated.
Instantaneous velocity contours on a sagittal plane and a coronal plane in the
pre-interventional geometry are given in Fig.
\ref{fig:phantom-instantFields.pdf}. Due to the relatively low temporal
resolution of 4D flow MRI, here only results from DNS and LES are presented. It
can be seen that the flow therein is turbulent. Because of the local narrowing
in the stenosis, flow is accelerated in the pre-interventional geometry and the
local Reynolds number on the stenosis plane is more than 2500. Jet flow, which
leads to high blood pressure and high wall shear stress (WSS), is observed. DNS
provides more flow details due to higher spatial resolution.

\begin{figure}
    \centering
    \center{\includegraphics[width=0.9\linewidth]  {DNSLES-pre-instantFlowField.pdf}}
    \caption{Instantaneous velocity contours on a sagittal plane and a coronal plane
	in the pre-interventional geometry of an aorta of a patient with CoA
	simulated by CFD. Left: DNS; right: LES.}
 \label{fig:phantom-instantFields.pdf}
\end{figure}

As the flow is unsteady, temporal averaging was performed for both CFD and MRI results
and the following comparison is based on the time-averaged flow fields.
The main flow features can be found in Fig. \ref{fig:phantom-streamlines.pdf}.
The visualization of 3D streamlines and velocity vectors on the sagittal plane
(insets) shows the complexity of the flow within the patient’s aorta,
especially where the stenosis occurs. Again, jet flow and recirculation are
observed in the pre-interventional geometry in the streamlines and highlighted
in the zoomed-in insets. Helical streamlines can also be found in all cases.
For the MRI results, some streamlines start from the vessel wall as no-slip
boundary condition is not guaranteed in MRI data. Nevertheless, all methods,
including DNS, LES and MRI, resolved the main flow features. Moreover, as the
aorta is deformed and flattened after stent implantation, the flow resistance
in the post-interventional geometry is reduced. The pressure drop is reduced
from 790 Pa (DNS) and 778 Pa (LES) to 9 Pa (DNS) and 8 Pa (LES), respectively.
Those results indicate that stent implantation restored the aortic flow
effectively.

\begin{figure}
    \centering
    \center{\includegraphics[width=0.9\linewidth]  {streamline-post-pre.pdf}}
    \caption{Streamlines of the flow computed with CFD simulations of the aorta
	compared to MRI measurements of the flow in the 3D-printed phantoms.
	The time-averaged flow is shown. The insets show the velocity vectors on
	the sagittal planes. Top: pre-interventional
	geometry; bottom: post-interventional geometry.}
	\label{fig:phantom-streamlines.pdf}
\end{figure}

Figure \ref{fig:MeanVelocityofSixPlanes} presents quantitative comparison of
mean velocity magnitude on six specified cross-sectional planes. These six
planes, as shown in the left panel of Fig. \ref{fig:MeanVelocityofSixPlanes},
represent the ascending, arch, pre-stenosis, on-stenosis, post-stenosis and
descending of the aorta respectively. The mean velocity magnitude was
calculated according to $ \sum_{i=1}^{N}\sqrt{U_{i}^{2}+V_{i}^{2}+W_{i}^{2}}$
$/N $, with $N$  the number of points on a cross plane. It can be seen from the
right panel of Fig. \ref{fig:MeanVelocityofSixPlanes} that the MRI results are
a little larger than the numerical ones on planes Pre 2, Pre 3 and Pre 5. Using
MRI results as reference, the relative deviations for LES and DNS are $4-28\%$
and $7-27\%$ respectively. Largest deviation can be found on Pre 5, mainly due
to the difficulty for the MRI measurement induced by the recirculation after
the stenosis. On Pre 4, MRI data is between the DNS ($-4\%$) and LES ($7\%$)
ones. For the post-interventional geometry, CFD results agree with the
experimental data well on all planes, with relative derivations less than
10$\%$.

\begin{figure}
    \center{\includegraphics[trim={0 0 0 0}, clip, width=0.9\linewidth]
    {MeanVelocityofSixPlanes.pdf}}
    \caption{Comparison of mean velocity magnitude on specified planes
	for the pre-interventional and post-interventional aorta geometry.
	Left: six specified cross planes, including ascending, arch, pre-stenosis,
	on-stenosis, post-stenosis and descending of aorta. Right: The
	time-averaged mean velocities from 4D flow MRI and the CFD simulations.} \label{fig:MeanVelocityofSixPlanes}
\end{figure}

Velocity contours on those specified planes, in pre-interventional and
post-interventional geometries, are given in Figs.
\ref{fig:pre-planes-VelocityMagnitude.pdf} and
\ref{fig:post-planes-VelocityMagnitude.pdf} respectively. Generally, flow
patterns are qualitatively well-matched among LES, DNS and MRI for both
geometries. Both LES and DNS provide more details and smoother fields due to
higher spatial resolution compared with MRI. Discrepancy between the CFD
results and MRI one is observed on Pre 4 and 5. We believe this is because of
the complexity induced by the jet flow and recirculation after the stenosis,
which is visualized in Figs. \ref{fig:phantom-instantFields.pdf} and
\ref{fig:phantom-streamlines.pdf}. Nevertheless, the results from LES agree
very well with those from DNS on all specified planes.

\begin{figure}
    \centering
    \center{\includegraphics[width=0.9\linewidth]  {pre-planes-VelocityMagnitude.pdf}}
    \caption{Flows in pre-interventional aorta. Velocity magnitude distributions were compared among CFD simulations and 4D Flow MRI on six planes as shown in Fig. \ref{fig:MeanVelocityofSixPlanes}. Time-averaged results.}
 \label{fig:pre-planes-VelocityMagnitude.pdf}
\end{figure}

\begin{figure}
    \centering
    \center{\includegraphics[width=0.9\linewidth]  {post-planes-VelocityMagnitude.pdf}}
    \caption{Flows in post-interventional aorta. Velocity magnitude distributions were compared among CFD simulations and 4D Flow MRI on six planes as shown in Fig. \ref{fig:MeanVelocityofSixPlanes}. Time-averaged results.}
 \label{fig:post-planes-VelocityMagnitude.pdf}
\end{figure}

Velocities along the centerlines of the pre-interventional and
post-interventional geometries were also compared. The centerline was extracted
using the VMTK extension for 3D Slicer \cite{PINTER201919}. Since velocity in
the pre-interventional geometry changes more acutely because of the stenosis,
here only the results for pre-interventional geometry are presented. Discrete
points along the centerline, starting from the ascending aorta after the
artificial extension, are considered and can be found in the left panel of Fig.
\ref{fig:Pre-centerline-U_V_W.pdf}. Velocity components and magnitude are given
in the right panel. It can be seen that velocity components change sharply in
the stenosis region which leads to a peak in the profiles of velocity
magnitudes. Results from LES, DNS, and MRI agree with each other in most region
of the aorta, except the coarctation. Similarly to Fig.
\ref{fig:pre-planes-VelocityMagnitude.pdf}, discrepancy mainly happens in the
MRI result around the stenosis, while good agreement between DNS and LES is
always achieved.

\begin{figure}
    \centering
    \center{\includegraphics[width=0.9\linewidth] {Pre-centerline-U_V_W.pdf}}
    \caption{Comparison of flow velocity along the centerline of the pre-interventional aorta. Left: locations of the points starting from the ascending aorta. Right: from top to bottom, velocity components and magnitude. Time-averaged results.}
 \label{fig:Pre-centerline-U_V_W.pdf}
\end{figure}

Based on the above comparisons, it can be concluded that CFD (LES and DNS)
results agree with the data from 4D-Flow MRI. More information, such as
pressure drop and WSS can be easily generated from the CFD results. Since LES
needs less computational resource but provides acceptable accuracy, we conclude
that LES is capable in resolving aortic flow and adopt it for the following
research.

\subsection{In vivo Validation}

To further validate the numerical modeling, we also used in vivo measurements
in addition to the phantom study. The patient-specific aortic flow was compared
between LES and flow MRI. Considering the long measurement duration needed by
4D-Flow MRI \cite{Diagnosis2014}, 2D flow MRI was used for the in vivo scans.
Velocities, thus flow rates, on specified planes perpendicular to the main flow
direction were measured with 2D phase-contrast flow MRI. The experimentally
obtained flow rate on a plane in the acceding aorta is shown in the left panel
of Fig. \ref{fig:invivo.pdf}. The two time instants, $t_1$ and $t_2$ during
systole, were considered and their flow rates were used as inlet boundary
condition for the LES simulations. The experimentally recorded velocity
distributions in the descending aorta at these two instants were used for
comparison.

In numerical simulations, the aortic geometry was segmented based on the MRI
slices. We didn't consider the deformation of the aorta, or fluid-solid
interaction, and assumed that the wall was not moving. A parabolic velocity
profile based on the experimentally recorded flow rate was specified as inlet.
The opening in the descending aorta was defined as outlet with a reference
pressure. Since the vessel branches were open in this test, the difference of
flow rates between the inlet and outlet were assigned to the branches according
to their cross areas \cite{Mirzaee2017, Pirola2017}.

The time-averaged flow fields are given in the right panel of Fig.
\ref{fig:invivo.pdf}. At $t_1$ (the instant with peak flow rate), the through-plane component of mean velocity from LES is 0.69 m/s, with relative deviation 8$\%$ in reference to the MRI result, 0.75 m/s. Similarly, at $t_2$, the through-plane component of mean velocities are 0.55 m/s (LES) and 0.49 m/s (MRI). Thus, good agreement between LES and MRI was achieved
also in vivo. Moreover, LES with proper boundary condition could provide more
flow details due to higher spatial resolution, as shown in both in vivo and
phantom tests.

\begin{figure}
    \centering
    \center{\includegraphics[trim={0 0 0 0}, clip, width=0.9\linewidth] {patientT1T2.pdf}}
    \caption{Comparison of in vivo aortic flow between LES and 2D flow MRI.
	Left: the experimentally recorded flow rate crossing the ascending
	plane; locations of the specified ascending plane and descending plane
	are given in the inset. Right: he through-plane component of velocity distributions on the descending plane at two different time instants during the cardiac cycle. }
	\label{fig:invivo.pdf}
\end{figure}

\section{Application for Stent Selection}

The main motivation for silico stent implantation is to help clinicians to
evaluate the surgical plans based on predicted results and to be able to select
an optimal stent already before surgery. A fast virtual stenting approach
proposed by Neugebauer et al. \cite{Neugebauer2016} was implemented in this
work to generate virtually deformed geometries. Several parameters, such as the
aorta bending resistance, aorta stiffness, stent stiffness, stent position and
diameter were considered in this approach to represent the interaction between
aorta and stent. The deformed geometry deformation was then obtained based on
the deformed centerline and deformed surface vertices. The original
pre-interventional geometry and its deformed versions can be found in the left
panel of Fig. \ref{fig:Deformed.pdf}. The deformed geometries were exported as
STL files, which were remeshed and further imported into the CFD solver for
flow simulation. It should be noted that in this study the physician chose a
stent with diameter 12 mm independently and without input from the in slicio
model. We compared the virtually deformed geometry based on this stent with the
post-interventional one reconstructed from MRI images in the right panel of
Fig. \ref{fig:Deformed.pdf}. Centerlines for both geometries were obtained and
the distance ( $\leq 1.60$ mm) between two centerlines is presented. Similarly,
the spatial dependent diameters are also given, with deviation less than 0.65
mm. It can be concluded that the virtually deformed geometry agrees with the
clinically deformed one quantitatively.

\begin{figure} \centering
	\center{\includegraphics[width=0.9\linewidth]  {Deformed.pdf}}
    \caption{Virtual stent implantation and validation of deformation. Left:
	pre-interventional geometry, virtually deformed geometries with stent
	diameter 8 mm, 12 mm, and 16 mm respectively. Right: validation of
	virtual stent implantation, with inset showing the realistic
	post-interventional geometry.} \label{fig:Deformed.pdf}
\end{figure}

The above validated LES was used to resolve the
flows in the pre-interventional geometry and the virtually deformed ones. We used
the same boundary conditions as mentioned in subsection 3.2, with flow rate
13.2 ml/s at the inlet. The numerically obtained results are presented in Figs.
\ref{fig:deformedAorta-streamlines.pdf}-\ref{fig:deformedAorta-WSS.pdf}. The
color coded streamlines in Fig. \ref{fig:deformedAorta-streamlines.pdf} show
that a stent with diameter 8 mm is inadequate to reduce the stenosis and jet
flow with a large local flow still can be observed in the narrowing region. For
stent diameter 12 mm or 16 mm, the jet flow disappears, with substantially
reduced maximum velocity compared to the pre-interventional geometry. The WSS
distributions are given in Fig. \ref{fig:deformedAorta-WSS.pdf}. WSS describes
the mechanical force generated by blood flow on the vessel wall, thus plays an
important role in chronic adaption and remodelling \cite{Humphrey2021}. It is
defined as $\tau_{y=0}= \mu \frac{\partial u }{\partial y}\mid  _{y=0} $, where
$\mu$ is the dynamic viscosity of the flow, $u$ is the flow velocity along the
wall and $y$ is the height above the wall.  As shown in Fig.
\ref{fig:deformedAorta-WSS.pdf}, high WSS is observed in the stenosis region in
the pre-interventional geometry and the deformed one with stent diameter 8 mm.
After the stenosis, high WSS is also found in a part of the descending aorta
due to the impact of high-speed jet flow (see Fig.
\ref{fig:deformedAorta-streamlines.pdf} for reference). A stent with diameter
16 mm enlarges the stenosis most and therefor leads to the smallest WSS in the
same region. However, as this stent is larger than the size of the aorta, it
also leads to a relatively large WSS before the stent, compared to the case
with stent diameter of 12 mm.

\begin{figure} \centering
    \center{\includegraphics[width=0.9\linewidth]  {streamline-of-deformed-aorta.pdf}}
    \caption{Streamlines in the pre-interventional and virtually deformed
	aortas. Time-averaged results. From left to right: pre-interventional
	geometry, virtually deformed geometries with stent diameter 8 mm, 12
	mm, and 16 mm respectively.} \label{fig:deformedAorta-streamlines.pdf}
\end{figure}

\begin{figure}
    \centering
    \center{\includegraphics[width=0.9\linewidth]  {WSS-deformed.pdf}}
    \caption{WSS on the pre-interventional and virtually deformed walls.
	Time-averaged results. From left to right: pre-interventional geometry,
	virtually deformed geometries with stent diameter 8 mm, 12 mm, and 16
	mm respectively. } \label{fig:deformedAorta-WSS.pdf}
\end{figure}

To quantitatively compare the four cases with different stent diameter, the
pressure drop and maximum WSS are given in Fig.
\ref{fig:deformedAorta-comparsion.pdf}. It can be seen that the pressure drop
is 119 Pa in the pre-surgical geometry, and is reduced to 34 Pa and 28 Pa in
the geometry with stent diameter 12 mm and 16 mm respectively. For the maximum
WSS on the aortic wall, the geometry with stent diameter 12 mm provides the
smallest value, 1.07 Pa. It is understandable that a larger diameter stent
results in less flow resistance thus smaller pressure drop, assuming that the
aortic wall is always deformable. On the other hand, size of the aorta wall and
its nonlinear response to possible strain should also be considered. If a stent
is too large for the aorta, it is conceivable that in addition to flow there
will be external mechanical force from the stent acting on the aortic wall.
Thus a stent with diameter 12 mm should be the optimal solution for the current
patient-specific aorta, which agrees with the physicians' independent choice in
this case.

\begin{figure}
    \centering
    \center{\includegraphics[width=0.55\linewidth]  {virtualDeformation-comparsion.pdf}}
    \caption{Quantitative comparison of pressure drop and maximum WSS in the aortas before and after deformation.}
 \label{fig:deformedAorta-comparsion.pdf}
\end{figure}

\section{Discussion and Conclusion}

Image-based in silico stent implantation
\cite{Pionteck2021,Chen2018,Kan2021,Neugebauer2016} and CFD
\cite{Zhong2016,Berg2018} together provide a new framework for stent planning
and interventional procedure evaluation. Besides a protocol for virtual
geometry deformation, a CFD method is needed to accurately resolve the flows in
the aortic geometries. However, blood flow in patient-specific aorta is
complicated \cite{Stein1976, Ku1997}. Laminar flow, turbulent flow and
transition between them may coexist spatiotemporally. In our study, we firstly
evaluated the accuracy of LES to predict such complicated flow. Two CFD
methods, the LBM based LES and DNS, in cooperation with flow MRI, were
considered. Both phantom and in vivo validations show that the LBM based LES,
which keeps a balance between numerical accuracy and computational requirement,
is a reasonable choice for resolving aortic flow. The validated LES was then
used to predict the flows in virtually deformed geometries with different stent
diameters. By comparing the flow fields, pressure drop, and maximum WSS, it was
found that the optimal stent was the one with diameter 12 mm, which agrees with
the physicians' independent choice.

To restore blood flow, in addition to numerical methods, accurate geometry and
boundary conditions are also important. Based on MRI scans, aortic geometry can
be segmented and reconstructed from high-quality image slices
\cite{M.A.Hadhoud2012,Avendi2016,Kan2021}. Furthermore, flow MRI is used for
visualization and quantification of aortic flow. 2D flow MRI with through-plane
velocity encoding is usually performed in  clinical applications
\cite{Diagnosis2014,Tanaka2010,Wentland2013}. But the 2D flow measurement is
affected by the selection of the cross-sectional plane. 4D flow MRI,
alternatively, is able to obtain time-dependent 3D blood flow, which is
resolved in all three dimensions of space and the dimension of time during the
cardiac cycle. It can used for the estimation of the flow pathways and the WSS.
But 4D flow imaging takes a significant amount of time, which prevents wide
clinical application. Thus for patient-specific in silico stent implantation,
MRI and LES need to work together and both are indispensable. Particularly, MRI
provides data for geometry and boundary conditions, while LES predicts aortic
flows for further evaluation.

In this work, we validated the LBM based LES with both phantom and in vivo
experiments, and provided a realistic example of the in silico stent implantation.
There are still many improvements possible which could be considered in the
future. Firstly, the aortic flow is unsteady, but we didn't consider time
dependent boundary condition in our simulations. We argue that current boundary
conditions are enough for us to compare different methods and provide results
for stent evaluation. By modeling the realistic cardiac cycle, one may get more
instantaneous flow information at the cost of longer computation time.
Secondly, we assigned the flow rates to the vessel branches according to their
cross areas. An alternative is the Windkessel model which considers resistance
and capacitance of the vessel network. Thirdly, we used a fast geometric method
to mimic the complex interaction between the aortic wall and the stent.
Ideally, one would take into account the mechanical properties of the aortic
wall and the stent, and model their interaction using finite element method.
Unfortunately it is still a big challenge to get accurate orthotropic
properties of the fiber reinforced aortic wall and to model the contact problem
numerically. Thus a simplified geometric method is a reasonable start. Finally,
this methodology could also be extended to other stenosis, such as cerebral
artery stenosis.

\section*{Conflict of interest}
The authors have no conflict to disclose.
\section*{Acknowledgments}
Computational resources from HPC at GWDG and MPCDF are appreciated.
Scientific comments and suggestions from the reviewers are also gratefully acknowledged.
\section*{Ethics Declarations}
This article does not describe studies with human or animal subjects.
The anonymous images used in this article were obtained from
a clinical database. Informed consent of the patients that permits
research use of their data has been obtained beforehand.

\bibliographystyle{apalike}