## Abstract

Cardiac catheterism is important because it offers many advantages in comparison to open surgery, for example, fewer injuries, lower risk of infections, and shorter recovery times. Simulators play a fundamental role in training packages, and virtual learning environments are less stressful. Moreover, they can also be used in certification boards and in performance assessments. A realistic and interactive simulator must be fast. In this work, the physical model of the guidewire used in catheter simulations has been improved. In particular, we determined a simple analytic expression to calculate the direction of a guidewire segment, which minimizes the total energy. The surface energy resulting from the guidewire--artery interaction and the bending energy of the guidewire is approximated up to the second order, which gives rise to interactions between segments. Furthermore, the multiple segment relaxations are introduced, enhancing the convergence especially at the beginning of the relaxation cycle. The formulas are written in matrix form of dimension $4M\xd74M$, where $M$ represents the number of segments varied in the update step. The method results in a more stable static solution.

## 1 Introduction

Over the last decades, minimally invasive surgery has revolutionized many surgical procedures (Basdogan, Kim, Muniyandi, & Srinivasan, 2004). Endovascular interventions include a variety of techniques that give access to the vascular system through small incisions. The cardiac catheterization is a procedure commonly used to diagnose and treat heart conditions (Balaji & Shah, 2011).

The catheter and the guidewire (GW) are manipulated within the vascular anatomy in order to navigate to the desired place in a blood vessel. This treatment is delivered using image guidance so that a skillful instrument navigation and a thorough understanding of the anatomy are critical in order to avoid complications. The physician must learn how to manipulate the catheter to the correct position (Wu, Chen, Chang, & Lin, 2016), and it is imperative to develop technically and conceptually the competences (Schneider, 2003). Moreover, it is mandatory to have wide experience in catheter-based subspecialties (Pelletier, Kaneko, Peterson, & Thourani, 2017). Medical companies have a vested interest in ensuring that their devices are appropriately and properly used.

The training methods include live observation of procedures, practicing on anatomic phantoms, and hands-on training using human cadavers or live animals (Cheng et al., 2017; Lenoir, Cotin, Duriez, & Neumann, 2006). Yet the phantom lacks realistic tissue and a complex human vasculature is expensive (Li et al., 2012). In the past, hands-on training was considered the best available method (Lunderquist et al., 1995; Mori, Hatano, Maruyama, & Atomi, 1998). However, it has ethical issues and it is also expensive, owing to the costs associated with the use of animals and because the instruments can be used only once (Coles, Meglan, & John, 2011). Moreover, animals are not a good substitute for humans because their vascular networks differ from humans' networks. In addition, real x-ray equipment can imperil the health of trainees.

The combination of traditional learning methods and technology enhances trainee satisfaction and skill acquisition levels (Engum, Jeffries, & Fisher, 2003; Tsang et al., 2008). Chui et al. (1996) developed an incremental finite element method (FEM) for simulation and analysis of catheter navigation. However, it was computationally intensive and not suitable for training in scenarios involving complex vascular models and catheter behavior. The first interventional-radiology simulator for catheterization training was DaVinci (Anderson, Raghavan, Wang, Mullick, & Kong, 1997). One year latter, DaVinci's techniques were extended to develop the interventional-cardiology simulator ICard (Lim, Shetty, Chui, Wang, & Cai, 1998; Wang, Chui, Lim, Cai, & Mak, 1998).

The simulation can be used as an initial training step to develop skills before further proceeding with training in animals or humans (Alderliesten, Konings, & Niessen, 2004). Several training sessions are performed with the aid of simulation techniques, providing certain levels of proficiency for the physician. Furthermore, simulators offer advantages: no radiation is required, specific cases can be handled with multiple difficulty levels, and there are no additional costs per training session. Simple devices can teach hand-eye coordination, while more sophisticated simulators can teach complex tasks and sequencing (Pelletier et al., 2017). Not to mention the fact that the virtual learning environment offers a less stressful and more controlled situation.

Simulators are able to differentiate advanced from novice operators, suggesting that it is a valid tool in the assessment of performances (Tsang et al., 2008). The virtual reality simulators are not used exclusively for training purposes, since they can be easily customized to provide both medical programs and certification boards with an objective tool for assessing physician skill and knowledge (Willaert, Aggarwal, Herzeele, Cheshire, & Vermassen, 2012). For instance, an important parameter which must be evaluated is the completion time for a specific task, because it is related to the total radiation exposure during the procedure. The physician must make prompt decisions, avoiding unnecessary advancements or rotations to minimize the completion time. Additionally, it is possible to measure the fluoroscopy time, the volume of contrast used to achieve a satisfactory result, and the placement accuracy of stents and valves (Pelletier et al., 2017). Surgical simulations are useful not only for training, but also helpful for pre-operative planning and intra-operative guidance (Bui, Tomar, & Bordas, 2019).

Computer models are increasingly being used to predict the behavior of instruments, making it possible to optimize their design (Cai et al., 2004). Moreover, the performance of instruments can be tested before the manufacturing phase. Until now, the selection of the GW has been often based on specialist's experience, which does not always result in a successful procedure (Al-Moghairi & Al-Amri, 2013). Simulations can be used to evaluate performances so that the GW with the proper mechanical properties to access the target location is chosen (Cardoso & Furuie, 2016).

Realism and computation time are two important requirements in GW modeling. Thus, it is crucial to evaluate and control the simulator errors so that the computational resources can be focused more on large instead of small error areas, particularly where clinically relevant (Bui, Tomar, Courtecuisse, Audette, et al., 2018). Moreover, it is impossible to validate a complete surgical simulation approach and to understand the sources of error without evaluating both the modeling and the discretization error (Bui, Tomar, Courtecuisse, Cotin, & Bordas, 2018). The best way to validate a model is by letting a specialist try it out and judge the outcome based on his or her real experiences (face validity) (Carter et al., 2005).

The GW relaxation is very fast, so it can be represented as a sequence of equilibrium states which depend on the *instantaneous* boundary conditions. The modeling is usually based on quasi-static mechanics, which is acceptable as the loading of the GW is slow and inertial effects can be ignored. The static equilibrium of the GW is determined by finding the minimum energy of a system with a very large number of degrees of freedom. Depending on the boundary conditions, there can be multiple minima. For example, a GW folded back inside an artery will not unfold, because the energy would first increase. If there are two different minima whose energies are nearly equal, a slight modification of the borders can move the system into the other solution. But it is very hard (time consuming) to perform this task numerically. The same phenomenon is also verified if the energy function is *flat*; that is, the energy gradient is small. Furthermore, the minimum can be difficult to reach if only one coordinate is varied at each step.

In general, the greatest concern with virtual reality simulators is not the exact solution of the problem but to obtain a reasonable solution in a short time interval, so that the simulation looks authentic and can be interactive. Therefore, it is indispensable to develop new numerical methods which can help to speed up the GW relaxation.

This article is organized as follows. First, existing GW models are briefly discussed and then the bending and the surface energy is calculated. The second-order surface energy term gives rise to the interaction between GW segments, which plays a fundamental role in the relaxations at multiple segments. Next the method is tested, optimized, and the influence of the higher order approximation on the physical model is analyzed. Then the proposed method is corroborated and the influence of surface parameters is analyzed. Finally, the main contributions are summarized.

## 2 Methods

Before presenting the method developed in this work, we define the GW model. In Section 2.2, an analytic expression for the GW segment orientation is determined, which minimizes the energy in the first-order approximation. Next, the minimization is carried out up to second order and the procedure is extended to multiple relaxation. Finally, an optimization of the parameters used in the multiple relaxation is performed.

### 2.1 Guidewire Models

The modeling of thin flexible objects, denoted as elastic rods, is an active field of research in computer graphics and animation (Bergou, Wardetzky, Robinson, Audoly, & Grinspun, 2008). In the field of medical simulations, elastic rods are frequently used to represent not only GWs but also threads, sutures, and catheters. If the GW is considered as a rod-like structure, there are different choices such as

Euler--Bernoulli beam theory (deformation due to bending),

Kirchhoff rod (the geometrically nonlinear generalization of the Euler--Bernoulli beam theory),

Timoshenko beam theory (deformation due to bending and to shear),

Cosserat rod (the geometrically nonlinear generalization of the Timoshenko beam theory).

To model the elastic rod deformations, nonlinear FEM (Duriez, Cotin, Lenoir, & Neumann, 2006; Li et al., 2001; Wang et al., 1998), mass-spring models (MSM) (Wang, Duratti, Samur, Spaelter, & Bleuler, 2007), and quasi-static approaches (Bosman & Alderliesten, 2005; Lawton, Raghavan, Ranjan, & Viswanathan, 2000) have been investigated.

The FEM is a numerical technique, which can be applied to model a deformable object like the GW inside the body. In this method, the equations incorporate the geometry and physical properties of the GW. Despite its numerical stability, the FEM requires a very high computational effort especially when the model is nonlinear (Meier, Lopez, Monserrat, Juan, & Alcañiz, 2005).

In MSM the instrument is considered as a network of masses connected to each other by springs/dampers. The main advantage of this method is its relative simplicity compared to FEM. However, it is not accurate and there is a risk of numerical instability, especially when the tissue is not soft in comparison with the instrument (Liu, Tendick, Cleary, & Kaufmann, 2003). Hence, the GW would demand high computational power to achieve stability, which is against the real-time requirements of simulators. The MSM is suitable only for pre-operative planning and designing purposes.

In rigid multibody link methods, the instrument is discretized into a set of rigid bodies connected to their neighbors by joints (Dawson, Cotin, Meglan, Shaffer, & Ferrell, 2000; Takashima et al., 2014). There are springs and dampers between the rigid bodies and the Newton--Euler equations are used to describe the translational and rotational dynamics. Moreover, it is easy to add other phenomena such as friction or material properties to each individual segment (Cai et al., 2004). A good approximation to the GW can be achieved using many small links to represent its high degree of flexibility. Since the rigid multibody link becomes slower as the length of the GW increases, the simulation is limited to a maximum number of segments, otherwise it will not run in real time. However, it is possible to have shorter segments in the distal end (because of more flexibility) and longer ones in the proximal end (due to more stiffness). This results in less computational time compared to MSM (Sharei, Alderliesten, Dobbelsteena, & Dankelman, 2018).

Lenoir, Meseure, Grisoni, and Chaillou (2002) simulated a surgical thread in real time using a dynamic spline. Cotin, Duriez, Lenoir, Neumann, and Dawson (2005) connected a set of beam elements that can model bending, twist, and other deformations. The model is based on a static finite element representation, which ensures an accurate simulation of devices such as catheters and GWs (Lenoir et al., 2006).

In contrast to the classical Frenet description, the Cosserat theory further considers the orientation of the rod cross-section. Pioneering work on Cosserat models was done by Pai (2002), which simulated the static configuration of a thread. This model takes into account all possible deformations of a one-dimensional object, but the numerical iterations to solve the problem are computationally expensive. Spillmann and Harders (2010) considered the effect of torsion in threads using a quaternion representation. They employed the method of Lagrange multipliers, which improves the stability of the method. However, assuming the torsion coefficient to be infinite results in a perfect torque control, which is close enough to reality (Sharei et al., 2018).

Tang, Wan, Gould, How, and John (2012) solved the dynamic deformations of GWs using a discrete differential geometry formulation. They considered a nonlinear elastic Cosserat rod (Antman, 2005) and used a parallel transport approach (Bergou et al., 2008) to compute elastic forces without stiff constraints, resulting in a stable and fast simulation. Ganji and Janabi-Sharifi (2009) used a forward kinematics approach to predict the catheter's tip position, but the model is non-physical (the catheter bends with constant curvature and there is no torsion).

Catheters and GWs are thin deformable objects that are inextensible in the axial direction and are prone to bend in the perpendicular direction. The GW model of the present work follows the idea of Konings, van de Kraats, Alderliesten, and Niessen (2003), which considers the GW as a discrete set of joints at positions $x0$, $\u2026$, $xN$, with $x0$ fixed. There are $N$ segments and the $i$-th segment $\lambda i=xi-xi-1$ is represented by a small rigid rod which is neither compressible nor bendable. Hence, the size $|\lambda i|=\lambda $ is the same for all segments (Alderliesten, Bosman, & Niessen, 2006). Nevertheless, it is possible to generalize the calculations to different sizes $\lambda $ as in the geometrical relaxation of Baier, Baier-Saip, Schilling, and Oliveira (2016).

When parts of the GW are inserted into the vessel, the current representation is adapted by adding segments and computing a new configuration with an optimization algorithm. The new segment is inserted just before the curved tip of the GW (Baier et al., 2016), because pushing the GW into the vessel mainly affects its distal end, while the rest of the scene looks almost static. The GW can also be rotated but this action only influences the tip which is intrinsically curved. The model allows for large rotations and the displacement is equal to one GW segment length $\lambda $, which is much larger than the step-size of Alderliesten et al. (2004).

### 2.2 Single Energy Minimization

The orientation of the GW is the result of interaction with the vascular wall and is mainly dominated by the forces experienced during propagation. These include the manipulation, contact with the vascular wall, and friction forces.

If friction is not zero, then a different path results in the vascular system (Alderliesten, Konings, & Niessen, 2007; Takashima et al., 2007). There are two forms of friction, kinetic and static, but since the velocities of the GW in the vessel are usually small, the velocity-dependent friction forces are neglected and the simulation is based on a quasi-static approach (Alderliesten et al., 2004; Tang et al., 2010; Tang, Wan, Gould, How, & John, 2012). For simplicity, the friction is not included here, although it can be incorporated in the model to calculate the force feedback in haptic devices. The friction would only change the total force acting on the GW, but the multiple segment relaxations developed in this work and its conclusions would remain the same.

Moreover, as in most of the studies on GW and catheter modeling, the effect of blood flow is neglected to reduce the complexity. However, in the presence of vascular malformations, modeling the blood flow might provide a better understanding of the pathological conditions (Cotin et al., 2000; Nowinski & Chui, 2001).

In order to compute the equilibrium of a segment, the sum of the forces arising from the bending and from artery--GW interactions must be equal to zero. Equivalently, this can be achieved through a minimization of the bending and surface energies, which is derived below.

#### 2.2.1 Bending Energy

*row*matrix and represents the transpose of the

*column*matrix $ui$. The scalar product of two vectors can be seen as the matrix product of a row matrix by a column matrix.

#### 2.2.2 Surface Energy

Biomechanical simulation with user interactions involves real-time computation of the deformation of soft tissues, collision detection, contact modeling, topological modification, and haptic feedback (Payan, 2012). The FEM provide high biomechanical realism, because the complex nonlinear behavior of soft tissue is directly accounted for through constitutive relations. However, in the case of catheterism, the artery deformations caused by the GW are usually tiny (Takashima et al., 2014), so that the relation between stress and strain is assumed to be linear.

Due to the discrete character of the GW representation, increasing $N$ some joints will come in contact with the surface and others will leave the surface; that is, the set $L$ will change as $\lambda $ varies. The surface energy is a continuous function of $\lambda $, but as one point enters in or leaves the surface, the derivative of the surface energy becomes discontinuous. Hence, the position $x\u2113$ of the minimum is shifted and there is a small jump in the minimum energy. For a high density of joints, the number entering minus the number leaving the surface is not *random* but follows a *regular pattern* as $\lambda $ varies. Furthermore, as $N$ increases the relative importance of a single joint decreases, so that the jump in the energy becomes smaller.

#### 2.2.3 Minimization

### 2.3 Multiple Energy Minimizations

The generalization to more than two simultaneous vector variations is straightforward. If $M$ represents the number of varied segments (VS) in a single step, then the matrix of dimension $8\xd78$ in Equation 33 is replaced by a matrix of dimension $4M\xd74M$.

#### 2.3.1 Separation between Relaxing Segments

Consider first the linear approximation. As shown in the discussion following Equation 22, $-b\u02dci.\u025bi$ will push the vector $(ui)new=ei,1+\u025bi$ in a direction parallel to $bi$. Normally $\kappa i-1\u223c\kappa i$ because the angle between $ui-1$ and $ui$ is not far from the angle between $ui$ and $ui+1$. Moreover, $-1\u2264\kappa i\u22641$ and so it is concluded that $7-\kappa i-1\u22487-\kappa i$. Hence, according to Equation 4 the vector $bi$ is nearly parallel to $ui-1+ui+1$. In other words, the bending energy will decrease if the unit vector $(ui)new$ is in the average direction of $ui-1$ and $ui+1$. Adding $s\u02dc(i,N).\u025bi$ to the energy changes the vector $bi$ to $qi=bi-s(i,N)$. In this way, the surface term pushes the unit vector in the direction of the applied surface force $-s(i,N)$ over the GW. There is a competition between the bending term which tends to keep the GW straight, and the surface term which tends to bend the GW so as to decrease the surface energy.

**Close segments:**If $j$ is close to $i$, then the probability that occurs a contact with the surface in this interval is negligible. Hence, it is likely that $s(i,N)=s(j,N)$, $R(i,N)=R(j,N)$, and the energy becomes

The double relaxation corrects half at segment $i$ and half at segment $j$. Since the quadratic term is positive, it is more advantageous to correct twice with $\u025bi=\u025bi(single)min/2$ and $\u025bj=\u025bi(single)min/2$ than to correct once with $\u025bi(single)min$.

If $j$ is not close to $i$, then the differences between $bi$, $bj$ and between $Ai$, $Aj$ will be large. As a result $\u025bi\xac\u2248\u025bj$ and the interaction term $\u025b\u02dci.R(j,N).\u025bj$ can be negative, decreasing the energy even more.

*close*and when $i$, $j$ are

*not close*, are illustrated in Figure 3a. Choosing the segments close from each other (e.g., $A$ and $B$), the bending forces will displace the GW in the same direction. However, the resulting displacement is opposed by the reaction force $-s1$ which will increase, because the GW is pushed against the surface. But if the segments are chosen at $A$ and $C$, the bending forces will displace the GW in opposite directions. The total vertical displacement to the right of $C$ is small, so neither the reaction force $-s1$ nor the reaction force $-s2$ increase substantially. Hence, $A$ and $C$ are

*almost free*to relax simultaneously, but this does not happen if $A$ and $C$ are relaxed separately because of reactions $-s1$ and $-s2$, respectively.

**Far segments:**Consider now the case in which the difference between $i$ and $j$ is large. If there are

*many*points in contact with the surface between joints $i$ and $j$, then

Figure 3b illustrates physically why multiple relaxation does not improve the outcomes when $j\u226bi$. Trying a single relaxation at $D$ will be hindered mainly by $-s4$, $-s5$, $-s7$, and $-s11$. Using a double relaxation at $D$ and $E$ will cancel the increase of $-s11$, but it has no effect on $-s4$, $-s5$, and $-s7$. The double relaxation is not efficient when the segments are on opposite sides of the contact point which causes troubles. Finally, a double relaxation at $D$ and $F$ cancels the increases of $-s4$, $-s5$, $-s7$, and $-s11$ affecting $D$. The increases of $-s6$ and of $-s10$ affecting $F$ are also canceled. Moreover, note that $-s8$ and $-s9$ do not come into play because they are perpendicular to the displacements.

#### 2.3.2 Algorithm

A very large number $M$ of VS should not be used, because the processing time to obtain and solve the systems of nonlinear equations increases. Depending on the number $M$ and the separations $\delta m$ ($m=1,\u2026,M-1$) between the VS, different *relaxation speeds* are obtained. The goal is to optimize $M$ and $\delta m$, so that the relaxation is as fast as possible.

Let $i$, $j$, $k$, $\u2026$ denote the segment indices where, in an update step, the multiple segment relaxations will be applied. Given the first index $i$, the other indices are obtained with $j=i+\delta 1$, $k=j+\delta 2$, $\u2026$ . In particular, note that for $Pij=0$ in Equation 31, the minimization problem of a double relaxation ($M=2$) becomes equivalent to two independent single relaxations ($M=1$) of segment $i$ and of segment $j$. Hence, if there is no contact with the surface between joints $i$ and $N$, no multiple relaxation is carried out. Hereunder, $Nlast$ is defined to be the last joint index where the GW is in contact with the surface.

It will be seen that in a double relaxation the fastest outcomes are obtained with $10\u2272\delta 1\u227220$. Here, a recipe is given to choose $\delta 1$ dynamically and the procedure is extended to higher order multiple relaxations.

Then the value of $j=i+5,i+6,\u2026,i+25$ is selected, which gives the smallest $\Delta Uij$. When $Nlast<i+25$ the comparison is carried out only for $j=i+5,i+6,\u2026,Nlast$. Further, if $Nlast<i+5$ no double relaxation is performed but just a single relaxation, because there is no substantial improvement in using the multiple relaxation when $j$ is close to $i$. Otherwise, $Nlast\u2265i+5$ and Algorithm 2 finds a reasonable segment index.

## 3 Results and Discussion

The major contribution of this work is the multiple relaxation which increases the numerical performance. The method is applied to a GW with $N$ segments following a sequential order, which is defined to be the relaxation cycle (Baier et al., 2016). The total number of updates in one cycle equals $12N(N+1)$. Furthermore, updating the $i$-th GW segment requires the evaluation of $N-i+1$ interactions with the surface.

The detection of each GW-surface collision and computation of the interaction force (Equation 12) takes on average 0.145 $\mu $s, but it increases slightly when $M>1$ because the matrix in Equation 30 must be calculated. In order to update the segments it is necessary to obtain and solve Equation 26, which takes about 10 $\mu $s for each segment. Again, if $M>1$ it is necessary to obtain and solve Equation 32.

The GW is a deformable object and as it slides along the vessel borders, a high number of contacts must be accurately modeled and computed. This is a challenge both in terms of accuracy and performance. In order to increase the processing speed, computing architectures have become increasingly parallel, with the ubiquitous use of multi-core CPUs, and the massive parallelism available in modern GPU (Courtecuisse & Allard, 2009). The latest generation of GPUs contains many processors and exploiting them fully requires creating many parallel tasks. Hence, using more powerful hardware and software, it is possible to decrease the computing time significantly so that a more realistic simulation can be achieved. Furthermore, increasing $\lambda $ and decreasing $N$ (i.e., decreasing the resolution of the GW to an acceptable level) results in a major reduction of the computation time.

In order to evaluate the efficiency of the method, the “lateral deviation” of the $i$-th joint $di$ is introduced. First, given an initially deformed GW, it is relaxed 200 cycles in order to obtain equilibrium. This reference curve (a 3D polyline) and the actual position of the $i$-th joint are projected in a plane and the minimum distance between them is found. As a measure of the deviation from equilibrium, the parameter $di$ is defined to be plus (minus) this distance if the joint is to the right (left) of the curve.

Updating the $h$-th segment changes the position of all the joints with $i\u2265h$, so that $di$ is calculated every time an update of a segment with index $h\u2208[1,i]$ is performed. In this way, the position of the joint relative to the reference curve can be tracked as a function of *time*. Note that for $h\u2208(i,N]$ the position does not change and $di$ is not recorded. Hence, the number of points (records) in a cycle equals $i(i+1)/2$ and the density of points in the relaxation increases with $i$.

Each GW segment has two degrees of freedom to describe the orientation of the unit vector $ui$ and one additional degree of freedom for the rotation around the symmetry axis. In order to make the analysis of the results simpler, the tubular artery is placed over a plane and only one projection is taken into account to track the position of the joint. However, it does not imply that the GW is planar because the artery is a 3D object.

Hereunder, the influence of the parameter $\delta m$ is analyzed in simple cases. Then the algorithm to determine $\delta m$ dynamically is evaluated and the stability of the method is investigated. Lastly, the convergence for increasing resolution is tested and the influence of the surface elastic constants on the GW shape is examined.

### 3.1 Parameter Optimization

Simulations with $M=2$, as in Figure 5a and in many other instances not shown here, have established that usually the fastest results are obtained for $10\u2272\delta 1\u227220$. On the other hand, for $2\u2264\delta 1\u22725$ and for $\delta 1\u227330$ the convergence is slow, but still faster than with $M=1$. Hence, there is no substantial improvement in using multiple relaxation if $\delta 1$ is very short or very large. Furthermore, choosing $M=3$, $\delta 1=\delta 2$ (Figure 5b) results in faster relaxations than with $M=2$, $\delta 1$ (Figure 5a) and in both cases the best outcomes are obtained for $10\u2272\delta 1\u227220$. However, note that for a smaller segment size $\lambda $ and a larger number of segments $N$, the optimal $\delta m$ increases because it is related to the separations of the contact points between the artery and the GW. The same effect is also observed if the artery has a bigger lumen diameter or a smoother axial curvature. In the next section, the possibility of different values for $\delta 1$ and $\delta 2$ chosen dynamically will be explored.

### 3.2 Algorithm Evaluation

$M$ . | 1 . | 2 . | 3 . |
---|---|---|---|

surface | 0.200 | 0.221 | 0.229 |

solution | 0.201 | 0.394 | 0.777 |

$M$ . | 1 . | 2 . | 3 . |
---|---|---|---|

surface | 0.200 | 0.221 | 0.229 |

solution | 0.201 | 0.394 | 0.777 |

*Note.* When the number of blocks $M$ is bigger than 1, it is necessary to consider additionally Equations 30 and 32. The calculations correspond to the GW in Figure 4, which has 200 segments. During one cycle the GW-surface interaction is computed $1.37\xd7106$ times and the solution is performed $20.1\xd7103$ times. The numerical experiments have been carried out with an Intel Core i7-4500U (1.8 GHz).

The parameter $d50$ of joint 50 is shown in Figure 6 for different numbers of VS. Comparing Figures 5 and 6, it can be seen that the double and triple relaxations with $\delta 1=15$ are indeed faster at the beginning, but the relaxations with the proposed algorithm are more stable.

A small segment index ($i\u226aN$) is particularly difficult to relax using a single relaxation, because *many* contact points oppose a rigid translation of the GW in one direction or another. However, using a double relaxation the two segments may relax in such a way that the sum of their displacements closely cancel and only a small rigid translation is carried out.

When three segments relax simultaneously, there is a greater chance that the sum of the displacements vanishes, so that the bending energy decreases without increasing the surface energy. In tests performed with $M=4$ and $M=5$ (not shown) the results have been very similar to $M=3$. But the processing time per cycle increases substantially, and in most cases the overall performance does not improve when $M>3$.

Lastly, Figure 8 shows parameter $d165$ and, since the relaxation is fast, only the 10 first cycles are shown. This segment is closer to the distal end of the GW and there are only few contact points which come into play. In general, it is easier to relax a segment closer to the end of the GW. The relaxations with $M=1$, 2, 3 do not take long and the fastest one is $M=3$. Moreover, in Figure 8 the noise is a bigger than in Figure 7, because the variation of $xi$ is the sum of the variations of $xi-1$ and of $\lambda ui$; that is, the noise builds up.

In a simulation, the physician does not wait 100 cycles to see the GW at equilibrium. Several actions can take place during a cycle: push, pull, or rotate the GW. The relevant part of the learning process is how the GW is guided inside the arteries and especially through the bifurcations (Grenon et al., 2011). It depends basically on the distal end of the GW and more relaxations are performed in this region, so that a realistic representation is obtained. In Figure 8, it can be seen that the relaxation is very fast at the beginning of the cycle. However, it is also important to achieve a fast relaxation in the rest of the GW, otherwise it will not look natural, for example, after a fast insertion (Baier et al., 2016).

Since the goal of the algorithm is to choose the separation between the VS, it is interesting to inspect how often the different values of $\delta m$ have been chosen in the relaxation. Figure 9 shows the relative frequency during the 100 cycles. In a double relaxation ($M=2$) the frequency of $\delta 1=5$ is somewhat higher than the frequency of $\delta 1=25$ and it shows a *monotonic* decreasing behavior between the extremes. On the other hand, for $M=3$ and $M=4$ some peaks appear, evidencing that the algorithm clearly chooses some key values. In particular, it can be deduced from the figure that the frequency of $\delta 2=5$ is large. After choosing the best value for $\delta 1=j-i$, the second $\delta 2=k-j$ will usually not be large because the probability decreases with the separation between segment $i$ and segment $k$.

Note that the combination of the physical and geometrical methods does not seem to be much faster than the physical method alone as described in Baier et al. (2016). The reason is that in the previous work, the distance is measured to the final equilibrium point and not to the equilibrium curve. After reaching the curve, the joint must *slide over the curve* to reach the final position and in this part the geometrical method is more efficient. However, what matters in the simulation is to see the joint over the equilibrium curve. When all joints are over this curve, they will also be in the equilibrium position.

It has been demonstrated that the multiple relaxations with a higher number of VS, relax in a shorter number of cycles. However, increasing $M$ also increases the processing time per cycle. Therefore, it is essential to know how the processing time depends on $M$. As an example, if $tM$ represents the average time of the first 10 cycles in Figure 4, then the following ratios are obtained: $t2/t1=1.52$, $t3/t1=2.67$, $t4/t1=4.83$. Thus, the processing time per cycle shows a significant increase when $M$ becomes larger. Considering the performance per cycle with different number of VS and the corresponding processing time per cycle, the best outcomes have been obtained with $M=2$ and $M=3$.

### 3.3 Static Solution

In Baier et al. (2016), it has been highlighted that in a static solution, the joints are not steady but they are bouncing on the surface of the artery: “*Only few points are effectively touching the surface, holding the wire to the equilibrium position. The points bounce in the artery and the total number of contacts vary from one step to the next.*” As pointed out in Section 2.2.2, this model is equivalent to a penalty method with $ks$ as coefficient. Using a high $ks$ may lead to an ill-conditioned system (Bui et al., 2019) and this seems to be the case in Baier et al. (2016). But a smaller $ks$ should not be used because the artery would be too soft.

Since the model has been slightly modified, it is interesting to examine what happens within this new formulation. In Figure 11 are the results of the calculations ($a$) with and ($b$) without the second order energy approximation, showing how much time each joint is in contact with the surface when the GW is in equilibrium. Furthermore, it displays the average modulus of the vector update $\u025bi=(ui)new-ui$. This figure is analogous to Figure 16 of Baier et al. (2016).

Comparing Figures 11a and 11b two remarkable differences appear. The first one is that in case ($a$), the joints are always in contact with the surface. This is not verified in case ($b$), although the joint indices are nearly the same in both cases. The second one is that the order of magnitude of $|\u025bi|$ is very different: in case ($a$) the magnitude is about $10-7$ times smaller than in case ($b$). Since $ui$ is virtually constant in case ($a$), the joints remain in the surface without bouncing. Thus, the calculation with quadratic terms is more stable.

In order to get a deeper insight into the origin of the difference, additional calculations have been performed varying the terms in Equation 17. The corresponding figures are very much like Figures 11a and 11b, and the results are only described.

($c$) Term $Ai$ was dropped but $R(i,N)$ was retained so that $Pii=R(i,N)$. The resulting figure is very similar to case ($a$): the joint indices in contact with the surface are the same and the curve of $|\u025bi|$ is comparable in shape and in magnitude.

($d$) Term $R(i,N)$ was dropped but $Ai$ was retained so that $Pii=Ai$. The figure is almost identical to case ($b$).

($e$) Moreover, the double and the triple relaxations with $Ai$ and $R(i,N)$ give the same outcome as in case ($a$).

Therefore, it can be concluded that the stability is due to matrix $R(i,N)$.

Before explaining this result, the relative size of $Ai$ and $R(i,N)$ will be compared. Both are symmetric positive definite matrices and the trace $TX=Tr(X)$ is used as a measure. Notice that if the spring constant $C$ in Equation 1 is the same for all GW segments, then $TAi$ does not depend on $i$. On the other hand, the value of $TRi$ is bigger for smaller $i$ and in the specific example considered here $TR1/TR200=18.0$. Further, it was found that the average of $TRi$ is 22.7 times larger than $TAi$, indicating that the second order surface energy contribution is more important than the second order bending energy contribution. Thus, the matrix $R(i,N)$ is dominant, especially for small $i$.

Figures 12b–d illustrate why the minimization in the second-order approximation is stable and the minimization in the first-order approximation is unstable. In the second-order approximation, the energy has a well defined minimum at $x2,min$ (Figure 12b), which is independent of $x0$. But in the first-order approximation, if $x0=xa$ is small (Figure 12c), then the energy has a minimum at $x1a,min$, which is bigger than $xa$. Thus, the joint is shifted to the right in the numerical minimization and the value of $x0$ increases. As a result, the new $U1st$ also increases and the minimum is shifted to the left; that is, an oscillating behavior is achieved. On the other hand, if $x0=xb$ is large (Figure 12d), then the minimum of the function for $x>0$ (blue curve) is at $x1b,min<0$ and the joint will leave the surface. In the next numerical step, the minimum of the function for $x<0$ (magenta curve) is at $xbend,min>0$ and the joint will move into the surface. Hence, the joint bounces on the surface and the outcome is an even worse oscillating behavior.

The inclusion of a second-order term in the surface energy makes the derivative of the energy continuous at the contact point, where the GW touches the artery. From a numerical point of view, it is simpler to find the minimum of a function with a continuous derivative. The discontinuity of the surface force in the first-order approximation causes problems and this instability propagates to the distal end of the GW. Since there are usually many points in contact with surface especially for small $i$, the end of the GW is notably affected.

It can be verified that within the first-order or within the second-order approximation, the bending energy is a smooth function. Hence, the inclusion of $Ai$ is not so relevant as $R(i,N)$. Lastly, since $Ai$ and $R(i,N)$ are positive definite matrices, they create a *potential well* with a minimum at the current position. Therefore, the modulus of update $\u025bi$ becomes smaller and the stability increases. No numerical artifact was introduced to bound $|\u025bi|$, but it arises naturally within the second-order energy approximation. This effect is more pronounced for small $i$ because, as it has been previously demonstrated, $R(i,N)$ is larger.

### 3.4 Verification of the Discrete Solution

The error in numerical simulation occurs in three ways (Bui, Tomar, Courtecuisse, Audette, et al., 2018). First, the inability of the model to represent the physical reality is expressed by the model error. Second, the inability of a numerical method to exactly solve the mathematical model leads to the discretization error. Third, the set of linear equations provided by the discretized model must be solved numerically. The associated error is known as the numerical error and includes the finite precision of computers and the round-off errors.

The GW deformation is a well studied problem in solid mechanics and the exact artery deformation is not relevant for the simulation, because it is usually much smaller than the GW deformation. Furthermore, the discretization error is much higher than the numerical error. Thus, the following analysis focuses on the discretization error and specifically demonstrates how the the size of $\lambda $ affects the results of the calculations. Note that $\lambda $ is the same for all GW segments; that is, in the present formulation it is not possible to increase the resolution of a specific wire interval while keeping the rest of the wire unchanged.

The *data scattering* in Figure 13 is real; that is, it is not due to round-off errors in the calculations. Indeed, varying the number of iterations the accuracy of the calculated $\Delta U/U$ has been estimated to be $\u223c3\xd710-5$. The same accuracy is obtained when the GW is perturbed and after that the equilibrium is found again. The change in $\Delta U/U$ observed in Figure 13 between consecutive values of $N'$ is on average $92\xd710-5$ with a maximum of $387\xd710-5$, which is much bigger than the accuracy of the calculations.

As pointed out at the end of Section 2.2.2, for bigger $N'$ the relative importance of a joint entering or leaving the surface becomes smaller, so that the jump from one point to the next should decrease as $N'$ increases. This is in accordance with the fact that the data in Figure 13 becomes *less disperse* as the resolution $\eta $ increases.

The accuracy of the results can be estimated increasing the resolution, that is, the number $N'$. The separation between the GW with 175 segments and the GW with 525 segments has an average value equal to 0.027 mm with a standard deviation of 0.041 mm. Hence, the calculated shape of the GW is almost identical in both cases. Since the standard deviation is bigger than the average deviation, the separation is smaller than the average for *most* points and considerably larger than the average at *few* points. In fact, a closer look (not shown here) reveals that the separation is bigger at larger GW curvatures. Please note that if the curvature is large, then the angle $\theta i$ between consecutive segments is also large and the error in the bending energy proportional to $\theta i6$ (Section 2.2.1) increases. On the other hand, at the contact points the GW is *clamped* and the separation between GWs calculated with different resolutions is small.

### 3.5 Parametric Uncertainty

The numerical simulation of soft tissue requires one to find a material model suitable to the patient under consideration and its parameters. In particular, the choice of parameters (material properties, boundary conditions, loading, etc.) used in a numerical model is one of the main difficulties (Hauseux, Hale, Cotin, & Bordas, 2018). The data is usually sparse, of unknown quality, and generated by different sources. Moreover, uncertainties are very difficult to determine because of inter-patient variability.

The least squares method is based on the squared difference between the measured data and the response of the model. In order to obtain the parameter values that give the best model response for the measured data, the squared difference is minimized with respect to the parameters. The result of such an approach is a deterministic estimate of the parameter values and the difference that remains at the identified values (the residual) is a measure for how well the model fits the experimental data. However, a true insight in the certainty of each identified parameter value is lacking (Rappel, Beex, Noels, & Bordas, 2019).

When considering the least-squares method, the confidence intervals are usually estimated under a normal distribution assumption (Elster & Wübbeler, 2015) and part of the experimental information is lost. On the contrary, Bayesian inference (BI) (Beck & Katafygiotis, 1998) constitutes a framework in which uncertainties in identified parameters naturally arise from the identification process itself under the form of a probability distribution function (PDF). Bayesian frameworks treat observations as realizations from a probability model and their output are probability distributions in terms of the model parameters (Rappel & Beex, 2019).

In the BI, initially assumed parameter PDF (prior) need to be specified, which are then updated using a likelihood function constructed from the observation data. The result is again a PDF (posterior), which represents the user's uncertainty about the model parameters and reflects the initial knowledge one has (Mohamedou et al., 2019).

If there are only few measurements (less than the number of parameters), the least-squares method has no unique solution. On the other hand, BI incorporates a regularization that makes cases with few measurements solvable (Rappel, Beex, & Bordas, 2018). Moreover, the parameters come with uncertainties and this is essential for the propagation of uncertainties in mechanical predictions. In case of a small number of measurements and obtaining an uncertainty, BI is unavoidable.

The Bayesian approach allows accounting for modeling uncertainty and the statistical noises of the experimental devices (Beck, 2010). It is possible to consider the output error, the input error, and the uncertainty of the model (Rappel, Beex, Noels, & Bordas, 2019). Incorporating both model uncertainty and the input error has a favorable influence on the identified parameters and the posterior predictions.

Different fibers often have different material properties, which are a realization from some PDF (Rappel & Beex, 2019). Not only the material randomness influences the macroscale responses of fibrous materials, but also the geometrical randomness. The BI can be used for the parameter identification of nonlinear constitutive models (Rappel, Beex, Hale, Noels, & Bordas, 2019). The approach is also employed to assess the quality of different models with respect to measured data (i.e., model selection), as for example, the hyperelastic constitutive models for soft tissue (Madireddy, Sista, & Vemaganti, 2015).

Furthermore, Koutsourelakis (2012) used BI to identify elastic material parameters which vary spatially. The intrusive formulation incorporates the various model equations in the likelihood (posterior) and is capable of inferring model discrepancies from noise in the data. In general, there are deviations between the physical reality (where measurements are made) and the idealized mathematical/computational description. Especially in the context of medical applications, it is important to account for the model discrepancy or inadequacy in order to infer the right material properties and make accurate diagnoses.

In forward uncertainty quantification, a constitutive material model is assumed a priori, but each material parameter follows a given probability distribution, for example, a normal distribution around a mean with a given standard deviation (Hauseux et al., 2018). The result of such a stochastic simulation is therefore also a probability distribution which can be characterized by statistical properties such as its mean and standard deviation. Given the fact that material parameters and models are so strongly variable for soft tissues, stochastic simulations can provide the user with confidence intervals. In a clinical environment, where safety-critical decisions must be made based on the output of simulations, it is of importance to propagate uncertainty through the model (Hauseux, Hale, & Bordas, 2017a).

The mechanical properties of tissues are very important in surgery simulations when tissue cutting is performed. However, in our case there are only small tissue deformations and they are not as significant as the GW deformations, whose physical model and parameters are well known. Nevertheless, through a Monte Carlo simulation, it will be estimated how much the GW deviates from the calculated deformation when the physical properties of the arteries are not homogeneous.

that is, a Fourier series with $G=9$ terms. The maximum index $\Omega =162$ is equal to the number of artery segments spanned by the GW in Figure 13a. The random parameters are in the intervals $ksg\u2208[0,2-gks]$ and $\varphi g\u2208[0,2\pi )$. In particular, with this choice for $ksg$ it follows that $0<Ks<2ks$; that is $Ks$ is always positive and bounded. The smallest and largest values are obtained with $ksg=2-gks$, $\varphi g=\pi $ and $ksg=2-gks$, $\varphi g=0$, respectively.

A few typical examples of $Ks$ are shown in Figure 14. Since $1\u2264\omega \u2264\Omega $ only half of the period is used and $Ks(\omega )$ does not exhibit a periodic behavior. Only the first 9 terms of the series are included, because higher-order terms would give rise to a non-smooth shape that is non-physical. Moreover, including the term $ks10$ would add a contribution smaller than $2-10ks$ to $Ks$ or less than 0.1% of $ks$. Also, it would be meaningless to make $G$ bigger than $\Omega /2$ as an arbitrary set of numbers ${Ks}$ containing $\Omega $ elements can be represented by a Fourier series with the first $\Omega /2$ terms. The PDF of $Ks$ generated in the Monte Carlo simulation is shown in Figure 15.

The numerical results have two sources of errors committed with respect to the undiscretized problem. The first error is due to the deterministic approximation and the second is due to the stochastic approximation (Monte Carlo estimator). The error in the deterministic approximation is far lower than that in the stochastic approximation, such that the error is dominated by the number of realizations used in the Monte Carlo estimator (Hauseux, Hale, & Bordas, 2017b).

For a given $Ks(\omega )$ curve, the equilibrium of the GW is calculated. Then the lateral deviation between the GW corresponding to the constant $ks$ (the reference GW) and the GW corresponding to $Ks(\omega )$ (the variable GW) is obtained. Figure 16 shows the time evolution of the Monte Carlo simulation of the mean deviation and of the rms value for some selected GW segments. It is seen that after 5000 steps the mean and the rms reach stable values, evidencing the convergence of the simulation.

## 4 Conclusions

It is crucial to have a fast relaxation so that the GW behavior looks realistic in a simulation, and a great effort has been done to increase the speed. Instead of using the traditional relaxation at a single segment, the method introduces multiple segment relaxations. Furthermore, the bending energy and the surface energy were approximated up to second order. In particular, the second-order term of the surface energy is responsible for a more stable convergence. The resulting formulas have been put in the form of matrices of dimension $4M\xd74M$, where $M$ represents the number of VS in an update step.

One important issue to be addressed in the method, is to decide the value of $M$ and of the separations $\delta m$ between the VS. In the simplest case, it has been found that the best outcomes are obtained for $10\u2272\delta 1\u227220$. For smaller values of $\delta 1$, two neighboring segments “tend” to relax in the same direction and the double relaxation does not improve substantially compared with the single relaxation. On the other hand, if $\delta 1=j-i$ is large, there may exist some joints in contact with the surface between $i$ and $j$ hindering the relaxation of the $i$-th segment. However, note that for a smaller segment size, a bigger lumen diameter, or a smoother axial curvature the optimal $\delta 1$ increases. An algorithm to choose dynamically $\delta 1,\delta 2,\u2026,\delta M-1$ in the multiple relaxation was proposed and analyzed.

The evaluations have shown that the number of cycles necessary to relax the GW decreases with $M$. Moreover, the number of cycles depends on the joint index being considered. It is usually tough to relax a segment with a small index because more GW points are in contact with the artery, hindering the relaxation. Since the processing time per cycle increases with $M$, it is not convenient to work with a large number of VS, and the best results have been obtained with $M=2$, 3. The GW relaxation is very fast, especially at the beginning of the cycle which is more relevant for dynamic simulations. Finally, it was shown that the model is not sensitive to the values of the elastic constants at the surface when they depend on the position, as in a real artery.

## Acknowledgments

This work was supported by INCT-MACC National Institute of Science and Technology on Medicine Assisted by Scientific Computing (CNPq, project 465586/2014-7).