Breaking Blood Flow with Wires in Aneurysm Coiling Treatment Simulations

By Fabian Holzberger, Jan Kirschke, Markus Muhr, Natalia Nebulishvili, Julian Schwarting, and Barbara Wohlmuth

Cerebral aneurysms are abnormal focal outpouchings of large intracranial arteries that result from structural weakening in the arterial wall (see Figure 1). Aneurysm rupture causes severe intracranial hemorrhage and has substantial rates of morbidity and mortality. Roughly 500,000 people around the world die from aneurysm rupture annually, and 30 percent of survivors suffer from severe neurological disorders. The pathogenesis of aneurysm formation involves multiple factors, including hemodynamic stress, wear and tear on the arterial wall, and high-flow states. Hypertension, smoking, and connective tissue disease may contribute to this process [4].

To mathematically describe the blood flow velocity field \(\underline{u}\) and associated pressure distribution \(p\) within the aneurysm sac and adjacent vessels, we employ the incompressible Navier-Stokes equations as a fluid model:

\[\begin{eqnarray} \frac{\partial\underline{u}}{\partial t}+(\underline{u}\cdot\nabla)\underline{u}+\frac{1}{\rho}(\nabla p-\nabla\cdot(2\mu(|\dot{\gamma}|)\epsilon(\underline{u})))&=&\underline{f}\\
\nabla\cdot\underline{u}&=&0. \end{eqnarray}\]

Here, \(\rho\) is the density of blood, \(\epsilon(\underline{u})\) is the viscous strain tensor, \(\mu\) is the blood’s viscosity depending on the shear rate \(|\dot{\gamma}|\), and \(\underline{f}\) is a force density. In their simplest form, these equations model a mass- and momentum-conservative fluid flow within a given (vascular) geometry that is subjected to boundary conditions at vessel walls, the inlet, and the outlet. For our numerical experiments, we obtained and adapted vascular geometries from freely available datasets that provide segmented and reconstructed three-dimensional (3D) geometries based on medical imaging (see Figure 2). Regarding boundary conditions, we consider the vessel walls to be rigid in accordance with the comparably small pulsation amplitude in the cerebral vessels of interest. We derive pulsatile inflow and outflow velocity and pressure conditions at the truncated domain ends from a 0D/1D model of the whole/global circulatory system [5], hence rendering a detailed, 3D “zoomed in” local simulation in the vicinity of the aneurysm. Figure 2a illustrates streamlines from a fluid simulation of an untreated aneurysm, where a clear circulation pattern is visible within the sack. This and subsequent images are snapshots from a heartbeat’s diastole phase.

To numerically solve this fluid mechanical system, we employ a lattice Boltzmann method (LBM) that is available in open-source software such as the waLBerla framework [1] in a highly parallelized fashion.

Coiling Treatment and Coil Formation Simulation 

Practitioners utilize two approaches to treat intracranial aneurysms: surgical clipping and the insertion of endovascular devices. While surgical clipping has been the standard of care in recent decades, the continual evolution of endovascular techniques is equally promising. In these minimally invasive procedures, doctors use a microcatheter to insert medical devices (i.e., platinum wires, stents, or flow diverters) into the cerebral vessels and/or aneurysm and trigger the occlusion of the aneurysm sac. During endovascular coiling, they fill the aneurysm cavity with platinum coils that trigger the occlusion in just a few minutes; Figure 3 displays angiographic images of the coiling process and Figure 4 depicts a simulation result. After surgical intervention, the endothelium may regrow to complete a successful aneurysm occlusion within months; however, inflammatory processes, high stresses in the aneurysm walls, and/or compaction of the thrombus can also trigger a recanalization that causes the patient to relapse. Due to its impact on the thrombus’ homogeneity and stability, the choice of coil highly influences the treatment’s overall success. Yet despite specialized endovascular procedures and advanced imaging techniques, device selection is mostly governed by the neuroradiologist’s personal experience [3, 7, 9].

We believe that numerical simulations of coil placement can support coil choice by estimating the outcome in silico with respect to coil shape, volumetric occlusion, and—by a subsequent hemodynamical fluid simulation—the coil’s influence on blood flow, aneurysm perfusion rate, and stresses on the aneurysm walls. We use discrete elastic rods [2] as a comparably lightweight approach to model the mechanics of a coil insertion procedure. We model the coiling wire as a one-dimensional (1D) curve that is represented by a set of discrete material points \(\underline{x}_i\) with mass \(m_i\). A connection (edge) forms between every two consecutive points \(\underline{x}_i\) and \(\underline{x}_{i+1}\); we then define a material coordinate frame on this connection to represent the orientation of the curves’ cross sections. Knowing the orientation and relative positions of the material points allows us to define measures of mechanical strain, which give rise to forces that act on the curve’s material points. In this case, strains stem from the influence of an external insertion force, contact with the aneurysm wall, the wire itself, and deviation between the wire’s natural shape (see Figure 4b) and its shape after insertion. The forces include stretching \(\underline{F}_{\textrm{s}},\) bending \(\underline{F}_{\textrm{b}},\) and torsion \(\underline{F}_{\textrm{t}}\). Finally, we can summarize the wire dynamics with a system of ordinary differential equations

\[m_i\ddot{\underline{x}}_i+k\dot{\underline{x}}_i=\underline{F}_{\textrm{f}}(\underline{F}_{\textrm{s},i}+\underline{F}_{\textrm{b},i}+\underline{F}_{\textbf{t},i})\qquad \textrm{for} \qquad i=1,2,\dots,n_{\textrm{masses}}\]

that is integrated in time via standard methods, with a friction law \(\underline{F}_{\textrm{f}}(\cdot)\) for the contacts. The parameter \(k\) introduces further damping that reflects the presence of blood around the coil during the insertion process. To reduce the risk of perforation during the procedure and fill the aneurysm’s volume more evenly, the wire’s natural shape is a helix (see Figure 4b). It is thus able to form loops of a specific radius when emerging from the microcatheter (see Figure 4a).

The 1D coiling structure can use the wire diameter of real devices to “inflate” to three dimensions, allowing us to place the fully resolved 3D coil as an obstacle for blood flow within the aneurysm geometry. By redoing the previous hemodynamics simulation, we can then make statements about the occlusion quality of a given coil and aneurysm. For example, Figures 2c and 2d demonstrate the flow pattern within an aneurysm with a fully inserted coil, which reduces the overall perfusion and suppresses the vortex pattern from the untreated case in Figures 2a and 2b. Such in silico operations yield a tool with which clinicians can assess the applicability of a treatment method before a real in vivo operation.

A Porous Medium Surrogate Model 

The final coil shape from both a numerical simulation and a real surgical operation is very sensitive to small parameter perturbations, like the insertion point or angle. Furthermore, a major drawback of the fully resolved blood flow analysis is the need for a fine-grained numerical grid resolution that captures all details of the coil obstacle. Instead of pursuing such a direct numerical simulation, we now aim to simplify the model to relevant, coarser features of a coil while maintaining a comparable overall influence on the hemodynamics.

The main idea is to replace the fully inserted, tightly packed coil with a porous volume of arbitrarily diffuse resolution that fills the aneurysm sac. This porous volume has a similar suppressing effect on the blood flow and also still incorporates the coil’s spatio-structural features with spatially heterogeneous porosity \(\phi=\phi(\underline{x})\) and permeability fields \(\mathbf{K}=\mathbf{K}(\underline{x})\), where \(\mathbf{K}\) is a symmetric, positive-definite permeability tensor. The porous medium is analytically and numerically reflected in our fluid model via a Darcy-Forchheimer term on the right side \(\underline{f}\) of the Navier-Stokes momentum equation [6]:

\[\underline{f}=-{\phi\nu}\mathbf{K}^{-1}\underline{u}-{\phi F_e}{\mathbf{K}^{-1/2}}|\underline{u}|\underline{u},\]

where \(\nu\) is an effective viscosity and \(F_e\) is some geometric scaling function. This representation slows and eventually blocks blood circulation within the aneurysm.

To obtain sensible porosity and permeability fields that resemble a coil’s structure, we employ a homogenization procedure that considers a local, volumetric neighborhood of sufficient size for each spatial (grid) point \(\underline{x}\) within the aneurysm geometry and computes the volume fraction of “in neighborhood and within the inflated coil volume” divided by “total neighborhood volume.” This fraction is called packing density \(c(\underline{x})\) and yields porosity via \(\phi(\underline{x})=1-c(\underline{x})\), which is merely a “diffuse/coarse shadow” of the original coil geometry (see Figure 5). We can influence the level of “diffusivity” by enlarging or reducing the size of the averaging neighborhood, thus allowing interpolation between a computationally demanding, extremely detailed geometry representation and (in extreme) a nearly homogeneous, constant porosity situation that could have much coarser resolution. We can then infer permeability \(\mathbf{K}(\underline{x})\) from porosity via Kozeny-Carman relations [8], which account for the local wire orientation to arrive at an anisotropic permeability tensor.

Simulation Results and Implications 

After creating both the fully resolved coil geometry and the porous medium surrogate, we compare the results of the corresponding hemodynamic simulations under two circumstances: (i) The untreated situation from before and (ii) among each other. Here, a common quantity of interest is the wall shear stress (WSS)—defined as \(\textrm{WSS}=\boldsymbol{\sigma}\underline{n}-(\underline{n}^{\top}\boldsymbol{\sigma}\underline{n})\underline{n}\), where \(\boldsymbol{\sigma}\) is the viscous stress tensor and \(\underline{n}\) is the outward normal of the aneurysm geometry—which is potentially a primary influential factor for aneurysm (re)growth or rupture risk. Figure 6 compares the WSS distribution on the aneurysm surface in three simulation cases. A clear reduction of WSS in the coiled case when compared to the untreated aneurysm indicates the treatment’s success. Although the porous medium surrogate model tends to slightly overestimate the WSS, both the fully-resolved coiled model and the porous medium model yield similar results — illuminating the porous medium’s applicability as a surrogate approach.

Beginning with a real-life and patient-specific aneurysm geometry, we have mechanically simulated a coil insertion procedure and obtained a 3D coil shape. Using a fully resolved hemodynamical LBM simulation, we are able to analyze this shape based on its influence on blood circulation reduction. We also derived a porous medium surrogate for coiling devices that—by a steerable degree of spatial heterogeneity—enables a balance between a rich(er) resolution of geometrical features and a coarser, more computationally efficient representation.

---

Markus Muhr delivered a minisymposium presentation on this research at the 2023 SIAM Conference on Computational Science and Engineering, which took place in Amsterdam, the Netherlands, earlier this year. More details about this project are available on its webpage.

Acknowledgments: This research was partially funded by the Deutsche Forschungsgemeinschaft (WO 671/11-1, WO 671/20-1) within the SPP2311, project number 465242983. 

References
[1] Bauer, M., Eibl, S., Godenschwager, C., Kohl, N., Kuron, M., Rettinger, C., … Rüde, U. (2021). waLBerla: A block-structured high-performance framework for multiphysics simulations. Comput. Math. Appl., 81, 478-501.
[2] Bergou, M., Wardetzky, M., Robinson, S., Audoly, B., & Grinspun, E. (2008). Discrete elastic rods. In ACM SIGGRAPH 2008 papers (pp. 1-12). Los Angeles, CA: Association for Computing Machinery.
[3] Campos, J.K., Lien, B.V., Wang, A.S., & Lin, L.-M. (2020). Advances in endovascular aneurysm management: Coiling and adjunctive devices. Stroke Vasc. Neurol., 5(1), 14-21.
[4] Etminan, N., & Rinkel, G.J. (2016). Unruptured intracranial aneurysms: Development, rupture and preventive management. Nat. Rev. Neurol., 12, 699-713. 
[5] Fritz, M., Köppl, T., Oden, J.T., Wagner, A., Wohlmuth, B., & Wu, C. (2022). A 1D–0D–3D coupled model for simulating blood flow and transport processes in breast tissue. Int. J. Numer. Method. Biomed. Eng., 38(7), e3612.
[6] Guo, Z., & Zhao, T.S. (2002). Lattice Boltzmann model for incompressible flows through porous media. Phys. Rev. E, 66(3), 036304.
[7] Pierot, L., & Wakhloo, A.K. (2013). Endovascular treatment of intracranial aneurysms: Current status. Stroke, 44(7), 2046-2054.
[8] Schulz, R., Ray, N., Zech, S., Rupp, A., & Knabner, P. (2019). Beyond Kozeny-Carman: Predicting the permeability in porous media. Transp. Porous Media., 130, 487-512. 
[9] Sindeev, S., Prothmann, S., Frolov, S., Zimmer, C., Liepsch, D., Berg, P., … Friedrich, B. (2019). Intimal hyperplasia after aneurysm treatment by flow diversion. World Neurosurg., 122, e577-e583.

	Fabian Holzberger is a Ph.D. student in the Numerical Mathematics research group at the Technical University of Munich (TUM). He studies the mechanical simulation of endovascular coiling devices in patient-specific geometries and works on models for blood clotting within coiled aneurysms that are suitable for the lattice Boltzmann method.
	Jan Kirschke is an experienced neurointerventionalist and a managing consultant in the Department of Neuroradiology at TUM. His research focuses on computer-aided quantitative imaging in neuroradiology and investigates the reliable, reproducible, and automatic extraction of clinically important image parameters.
	Markus Muhr is a postdoctoral researcher in the Numerical Mathematics research group at TUM. He earned his doctoral degree from TUM in 2022, and his research interests lie in numerical simulations of seismic wave phenomena and medical applications such as nonlinear acoustics, ultrasound lithotripsy, and lattice Boltzmann methods for hemodynamics.
	Natalia Nebulishvili is a Ph.D. student in the Numerical Mathematics research group at TUM, where she studies porous medium surrogate models for endovascular treatment devices and non-Newtonian models for blood flow that are suitable for a lattice Boltzmann method framework.
	Julian Schwarting is a clinician scientist and resident in the Department of Neuroradiology at TUM. His research focuses on the translational investigation of subarachnoid hemorrhage and treatment strategies for patients with cerebral aneurysms.
	Barbara Wohlmuth is Chair of Numerical Mathematics and director of the International Graduate School of Science and Engineering at TUM. She is a member of the German National Academy of Sciences Leopoldina. Her research interests include mathematical modeling, numerical analysis of partial differential equation systems, high-performance algorithms, and geophysical and biomedical applications.