**3. Modeling of coil deformations during embolization**

In this section, we describe the deformation model that is used to capture the mechanical properties of the coil. Moreover, as the behavior of the coil during embolization strongly depends on the contacts with the vessel wall, we presents a constraint-based collision response.

#### **3.1. Coil model**

There are different types of detachable coils but most of them have a core made of platinum, and some are coated with another material or a biologically active agent. All types are made of a soft platinum wire of less than a millimeter diameter and therefore are very soft. The softness of the platinum allows the coil to conform to the arbitrary shape of an aneurysm.

The deformation model of the coil is based on the recent work of Dequidt et. al. [6]. Coil dynamics is modeled using serially linked beam elements:

$$\mathbf{M}\dot{\mathbf{v}} = \mathbf{p} - \mathbb{F}\left(\mathbf{q}, \mathbf{v}, \mathbf{q}\_0\right) + \mathbf{H}\mathbf{f},\tag{1}$$

where **<sup>M</sup>** <sup>∈</sup> **<sup>R</sup>**(*n*×*n*) gathers the mass and inertia matrices of beams. **<sup>q</sup>** <sup>∈</sup> **<sup>R</sup>***<sup>n</sup>* is the vector of generalized coordinates (each node at the extremity of a beam has six degrees of freedom: three of which correspond to the spatial position, and three to the angular position in a global reference frame). The rest position **q**<sup>0</sup> represents the reference position. Tuning the rest position allows to simulate different families of coil: for instance, helical shape or more complex 3D shape (see Figure 4). **<sup>v</sup>** <sup>∈</sup> **<sup>R</sup>***<sup>n</sup>* is the vector of velocity. **<sup>F</sup>** represents internal visco-elastic forces of the beams, and **p** gathers external forces. **f** is the vector of the contact forces with the aneurysm wall, and **H** gathers the contact directions. To integrate this model we use backward Euler scheme with a unique linearization of **F** per time step. Moreover the linear solver takes advantage of the nature of our model. All beam elements being serially connected, **F** is a tridiagonal matrix with a band size of 12. nd we solve the linear system using the algorithm proposed by Kumar *et al.* [20]. The solution can thus be obtained in *O*(*n*) operations instead of *O*(*n*3).

<sup>1</sup> http://www.cgal.org/

**Figure 4.** Example of coils used in our simulations: (left) Boston Scientific helical coil GDC 10; (right) 3D GDC built with omega loops [2].

#### **3.2. Simulation of coil deployment**

10 Will-be-set-by-IN-TECH

1. Manually place initial blobs inside the aneurysm and related blood vessel. Only an

Once we get a smooth surface model, meshing the domain is done using the interleaved optimization algorithm based on Delaunay refinement and Lloyd optimization [41], and implemented using the CGAL library 1. Each refinement step acts on the size of elements, while each optimization step acts on the shape of elements. Using this method, we can also define a size field to obtain anisotropic and non-uniform mesh. We control the fidelity of the generated mesh to the previously obtained surface model, by specifying the distance tolerance

In this section, we describe the deformation model that is used to capture the mechanical properties of the coil. Moreover, as the behavior of the coil during embolization strongly depends on the contacts with the vessel wall, we presents a constraint-based collision

There are different types of detachable coils but most of them have a core made of platinum, and some are coated with another material or a biologically active agent. All types are made of a soft platinum wire of less than a millimeter diameter and therefore are very soft. The softness of the platinum allows the coil to conform to the arbitrary shape of an aneurysm. The deformation model of the coil is based on the recent work of Dequidt et. al. [6]. Coil

where **<sup>M</sup>** <sup>∈</sup> **<sup>R</sup>**(*n*×*n*) gathers the mass and inertia matrices of beams. **<sup>q</sup>** <sup>∈</sup> **<sup>R</sup>***<sup>n</sup>* is the vector of generalized coordinates (each node at the extremity of a beam has six degrees of freedom: three of which correspond to the spatial position, and three to the angular position in a global reference frame). The rest position **q**<sup>0</sup> represents the reference position. Tuning the rest position allows to simulate different families of coil: for instance, helical shape or more complex 3D shape (see Figure 4). **<sup>v</sup>** <sup>∈</sup> **<sup>R</sup>***<sup>n</sup>* is the vector of velocity. **<sup>F</sup>** represents internal visco-elastic forces of the beams, and **p** gathers external forces. **f** is the vector of the contact forces with the aneurysm wall, and **H** gathers the contact directions. To integrate this model we use backward Euler scheme with a unique linearization of **F** per time step. Moreover the linear solver takes advantage of the nature of our model. All beam elements being serially connected, **F** is a tridiagonal matrix with a band size of 12. nd we solve the linear system using the algorithm proposed by Kumar *et al.* [20]. The solution can thus be obtained in *O*(*n*)

**Mv**˙ = **p** − **F** (**q**, **v**, **q**0) + **Hf**, (1)

approximate radius is needed for these blobs.

3. Run the above implicit modeling algorithm

between the two surfaces.

response.

**3.1. Coil model**

operations instead of *O*(*n*3).

<sup>1</sup> http://www.cgal.org/

2. Extract data points by casting rays from the blob centers

**3. Modeling of coil deformations during embolization**

dynamics is modeled using serially linked beam elements:

#### *3.2.1. Modeling contacts with aneurysm walls*

Simulating coil embolization requires to accurately model contacts that occur between the coil and the aneurysm wall. The contact model must provide the following features: first, account for the stick and slip transitions that take place during the coil deployment, second include a compliant behavior of vessel wall (we choose the one that is close to Boussinesq model [30]) and finally the friction motion of the coil along the aneurysm wall. For modeling contacts with friction, we use two different laws, based on the contact force and on the relative motion between the coil and the aneurysm walls. The contact law is defined along the normal **n** and the friction law, along the tangential (**t**, **s**) space of the contact.

The contact model, based on Signorini's law, indicates that there is complementarity between the gaps *δ<sup>n</sup>* and the contact forces f*<sup>n</sup>* along the normal direction, that is,

$$0 \le \delta^n \perp \mathbf{f}^n \ge 0.$$

With the Coulomb's friction law, the contact force lies within a spacial conical region whose height and direction are given by the normal force, giving two complementarity conditions for stick and slip motion:

$$\begin{aligned} \left[\delta^{t}\;\delta^{s}\right] &= \mathbf{0} \Rightarrow \left| \left[\mathbf{f}^{t}\;\mathbf{f}^{s}\right] \right| < \mu\left|\left[\mathbf{f}^{\mu}\right] \right| \quad \text{ (stick condition)}\\ \left[\delta^{t}\;\delta^{s}\right] &\neq \mathbf{0} \Rightarrow \left[\mathbf{f}^{t}\;\mathbf{f}^{s}\right] = -\mu\left|\left[\mathbf{f}^{\mu}\right] \right| \frac{\left[\delta^{t}\;\delta^{s}\right]}{\left\|\left[\delta^{t}\;\delta^{s}\right] \right\|} \text{ (slip condition)} \end{aligned}$$

Where the vector [*δ<sup>t</sup> δs*] provides the relative motion in the tangential space and *μ* represents the friction coefficient.

The obtained complementarity relations could create *singular* events when it changes from one state to another: For instance, when a collision occurs at instant *t �*, the velocity **v**(*t �*) of the coil, at that point, changes instantaneously. The acceleration could then be ill-defined and we can observe some quick changes in the dynamics. Each friction contact creates three non-holonomic constraints along the normal and tangential directions. Our approach allows for processing simultaneously multiple friction contacts, including self-contacts on the coil.

#### *3.2.2. Simulation steps*

The processing of one simulation step begins with solving Equation 1 for all forces except contact forces (**f** = 0). This *free motion* corresponds essentially to the deformation of the beam elements under gravity and user force input. Once the *free motion* has been computed, collision detection computes the contact points between the coil model and the aneurysm

**Figure 5.** Examples of our simulation results: (left) real coil embolization; (right) our simulated coil embolization with 3D coils.

surface and the points of self-collision. When collisions are detected, the *contact response* is computed. This is a complex aspect that influences greatly the overall behavior of the coil model. To describe the mechanical behavior during contact, the mechanical coupling between the different contact points is modeled. This information is provided by evaluating the compliance matrix in the contact space, called **W**, for both the coil and the aneurysm. Let's consider a contact *α* on the node *i* of the coil (with one constraint along the contact normal **n** and two along the tangential friction directions **t**, **s**). **H***<sup>α</sup>* is the matrix of the frame [**nts**]. The mechanical coupling of this contact with a contact *β* (with frame **H***β*) on node *j* can be evaluated with the following 3 × 3 matrix:

$$\mathbf{W}\_{(\boldsymbol{\alpha},\boldsymbol{\beta})} = \mathbf{H}\_{\boldsymbol{\alpha}}^{T} \left( \frac{\mathbf{M}}{h^{2}} + \frac{d\mathbf{F}}{hd\mathbf{v}} + \frac{d\mathbf{F}}{d\mathbf{q}} \right)^{-1}\_{(i,j)} \mathbf{H}\_{\boldsymbol{\beta}} = \mathbf{H}\_{\boldsymbol{\alpha}}^{T} \mathbf{C}\_{(i,j)} \mathbf{H}\_{\boldsymbol{\beta}\boldsymbol{\gamma}}$$

where **C**(*i*,*j*) is the 3 × 3 sub-matrix of global compliance matrix **C** (inverse of tangent matrix) at the rows of node *i* and the columns of node *j*. For the aneurysm wall, the formulation of the coupling is simpler:

$$\mathbf{W}\_{(\boldsymbol{\alpha},\boldsymbol{\beta})} = \frac{\mathbf{g}(d\_{ij})}{e} \mathbf{H}\_{\boldsymbol{\alpha}}^T \mathbf{H}\_{\boldsymbol{\beta}\boldsymbol{\epsilon}}$$

where *e* is an elasticity parameter that is homogeneous to young modulus and *g*(*dij*) is a Gaussian function of the distance, defined on the surface, between contact point *i* and *j*. The Gaussian function allows a fall-off of the coupling with increasing distance between the contact points. This model is close to the Boussinesq approximation which provides a distribution of the normal contact stress from the elasticity of the surface, around a point of contact [30].

The result of the *contact response* involves finding the friction contact forces that respect Signorini and Coulomb laws. Several works ([18] and [8]) present Gauss-Seidel iterative approaches that solve this problem. The solver needs an evaluation of a global compliance matrix **W**, which is the sum of the compliances of the coil and the aneurysm wall. It also needs the value of the relative displacement of the contacting points during the free motion δfree. When the contact forces are found, during the last step, called *contact correction* we compute the motion associated to the contact forces.

In this section, we presented an efficient dynamic coil model, *i.e.*, FEM beams based model for the intrinsic mechanical behavior of the coil, combined with a process for computing collision response (with the aneurysm wall) and auto-collision response. However, even if the results are quite encouraging, as illustrated in Figure 5, the mechanical modeling is not fully complete. Indeed, the behavior of the coil is also greatly influenced by the blood flow. Moreover, we also aim at predicting the influence of the coils on the hemodynamic status near the aneurysm. As a result, it is important to simulate the flow within the aneurysm, and to propose a method for blood-coil coupling.
