## 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×4M$, 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$, $…$, $xN$, with $x0$ fixed. There are $N$ segments and the $i$-th segment $λi=xi-xi-1$ is represented by a small rigid rod which is neither compressible nor bendable. Hence, the size $|λi|=λ$ is the same for all segments (Alderliesten, Bosman, & Niessen, 2006). Nevertheless, it is possible to generalize the calculations to different sizes $λ$ 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 $λ$, 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

Let $ui$ be a unit vector parallel to the GW segment $λi$; that is, $λi=λui$ where $λ$ is a constant representing the length of the segments. The rod corresponds to an Euler--Bernoulli beam which is homogeneous, isotropic, and follows a linear elastic stress-strain law; that is, $σ=Eɛ$. The bending energy of the arc at the $i$-th joint (see Figure 1) is given by Baier et al. (2016):
$Ubend,i(θi)=12EIRi2Riθi=EIλθisinθi2≈C1213-14cosθi+cos2θi,$
where $EI$ is the flexural rigidity and $C=EI/λ$ represents a spring constant. The error in the approximation is of the order of $31280Cθi6$. Hence, to a high degree of accuracy the energy can be approximated by
$Ubend,iui,ui+1=C1213-14u˜i.ui+1+(u˜i.ui+1)2.$
(1)
The vector $u˜i$ is a 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.
Figure 1.

Bending energy of the GW segment (red arc). Three consecutive joints are located at the coordinates $xi-1$, $xi$, and $xi+1$ (blue dots). Observe that $λui=xi-xi-1$ and $λui+1=xi+1-xi$, where $λ$ represents the distance between joints and the vectors $ui$, $ui+1$ are unitary. The angle $θi$ of the arc can be calculated with the formula $cosθi=u˜i.ui+1$ and the radius with the formula $Ri=λ/[2sin(θi/2)]$.

Figure 1.

Bending energy of the GW segment (red arc). Three consecutive joints are located at the coordinates $xi-1$, $xi$, and $xi+1$ (blue dots). Observe that $λui=xi-xi-1$ and $λui+1=xi+1-xi$, where $λ$ represents the distance between joints and the vectors $ui$, $ui+1$ are unitary. The angle $θi$ of the arc can be calculated with the formula $cosθi=u˜i.ui+1$ and the radius with the formula $Ri=λ/[2sin(θi/2)]$.

In an update step, the orientation of the $i$-th segment is changed $ui→ui+ɛi$ and the other unit vectors $uj$ are kept constant, in such a way that the total energy decreases. The new bending energy at the $i$-th joint is
$Ubend,iui+ɛi,ui+1=C1213-14(u˜i+ɛ˜i).ui+1+(u˜i+ɛ˜i).ui+12=Ubend,iui,ui+1+C12[-14ɛ˜i.ui+1+2κiɛ˜i.ui+1+(ɛ˜i.ui+1)2],$
(2a)
with $κi=u˜i.ui+1$. Further, at the $(i-1)$-th joint the energy is also modified
$Ubend,i-1ui-1,ui+ɛi=Ubend,i-1ui-1,ui+C12[-14ɛ˜i.ui-1+2κi-1ɛ˜i.ui-1+(ɛ˜i.ui-1)2].$
(2b)
The bending energy associated with the $i$-th segment is the sum of Equations 2a and 2b
$Ubend(ɛi)=Ubend(0)+C6(-7+κi)u˜i+1.ɛi+C12ɛ˜i.ui+1u˜i+1.ɛi+C6(-7+κi-1)u˜i-1.ɛi+C12ɛ˜i.ui-1u˜i-1.ɛi=Ubend(0)-b˜i.ɛi+12ɛ˜i.Ai.ɛi,$
(3)
the vector $bi$ and the tensor $Ai$ being defined by
$bi=C6(7-κi-1)ui-1+(7-κi)ui+1Ai=C6ui-1u˜i-1+ui+1u˜i+1.$
(4)
Notice that $uu˜$ represents a vector direct product (a dyadic) and it can be viewed as a square $3×3$ matrix.
Increasing the number of GW segments $N$ and decreasing the segment length $λ$ in such a way that $Nλ=$ const (the total length of the GW), should improve the representation of the GW. It is therefore advisable to analyze what happens with the energy when the resolution increases. Using Equation 1 the total bending energy is
$Ubend=∑i=0N-1C1213-14cosθi+cos2θi,$
(5)
For small angles $θi$ the energy can be approximated by
$Ubend=∑i=0N-1EI2λθi2-θi424=NEI2λθ2-θ424,$
(6)
where the angle brackets represent averages.
Increasing the number of GW segments by a factor $η=N'/N$, decreases the segment size to $λ'=λ/η$ and the angle approximately to $θ'=θ/η$. Hence, the new energy is
$Ubend'=N'EI2λ'θ'2-θ'424=NEI2λθ2-θ424η2.$
(7)
The relative energy variation can be calculated with the approximate formula
$ΔUbendUbend=Ubend'-UbendUbend=1-1η2θ424θ2,$
(8)
showing that the energy has a small increase as the resolution $η$ increases.

#### 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.

Changing $λi$ by $λɛi$ affects the coordinates $xi,…,xN$ by the same amount; that is, $xℓ→xℓ+λɛi$ for $ℓ∈[i,N]$. Hence, due to the interactions between the GW and the artery, the surface energy also varies. Let $xℓ$ be in contact with the artery which is discretized by a set of triangles with linear springs at the vertices (Baier-Saip, Baier, Schilling, & Oliveira, 2017; Baier, Srinivasan, Baier-Saip, Voelker, & Schilling, 2015). Each surface element is characterized by a plane normal $nℓ$ (pointing outwards from the artery), a point $rℓ$ on this plane, and an effective spring stiffness $ks$. The interaction energy is given by
$Usurf,ℓ(xℓ+λɛi)=12ksn˜ℓ.(xℓ+λɛi-rℓ)2=12ksn˜ℓ.(xℓ-rℓ)2+ksλn˜ℓ.(xℓ-rℓ)n˜ℓ.ɛi+12ksλ2ɛ˜i.nℓn˜ℓ.ɛi=Usurf,ℓ(xℓ)+s˜ℓ.ɛi+12ɛ˜i.Rℓ.ɛi,$
(9)
where
$sℓ=ksλn˜ℓ.(xℓ-rℓ)nℓRℓ=ksλ2nℓn˜ℓ.$
(10)
In particular, the contact force between the GW and the artery, that is, the gradient of the surface energy
$Fℓ=sℓλ=ksn˜ℓ.(xℓ-rℓ)nℓ$
(11)
can be viewed as a penalty method (Courtecuisse et al., 2014). The spring stiffness $ks$ is termed as a penalty coefficient and $n˜ℓ.(xℓ-rℓ)$ is the interpenetration distance.
If there are several points in contact with the surface in the interval $i≤ℓ≤N$, then $sℓ$ and $Rℓ$ must be replaced by the sum
$s(i,N)=∑ℓ=iNsℓR(i,N)=∑ℓ=iNRℓ.$
(12)
Observe that if $n˜ℓ.(xℓ-rℓ)<0$ the joint does not touch the surface and $sℓ$, $Rℓ$ are set equal to zero. Finally, the surface energy becomes
$Usurf(ɛi)=Usurf(0)+s˜(i,N).ɛi+12ɛ˜i.R(i,N).ɛi.$
(13)
Increasing the number of GW segments $N$ and decreasing the segment length $λ$, varies the surface energy at equilibrium
$Usurf=∑ℓ∈L12ksn˜ℓ.(xℓ-rℓ)2,$
(14)
where the sum extends only over the set $L$ of joints in contact with the surface. Assuming that in a first approximation the shape of the GW does not change as the resolution varies, the position $xℓ$ is a continuous function of $λ$ which can be taken as a curve parameter. Decreasing $λ$ moves the joints toward the beginning of the GW and new ones appear at the distal end (see Figure 2). Since the density of joints increases, there will be more joints in contact with the surface. In order to minimize the sum of the bending energy and the surface energy, the shape of the GW curve must change slightly, varying the bending energy.
Figure 2.

GW in contact with the surface. The left joint (blue) is fixed and the shape of the GW does not change. As the separation between the joints $λ$ decreases, the joints move to the left as indicated by the arrows. ($a$) The red joint and the green joints are inside the surface and outside the surface, respectively. ($b$) As the joints move, the red joint gets out of the surface and the green joints get in the surface. Furthermore, in order to keep the total GW length constant, new joints are created at the right (open circles).

Figure 2.

GW in contact with the surface. The left joint (blue) is fixed and the shape of the GW does not change. As the separation between the joints $λ$ decreases, the joints move to the left as indicated by the arrows. ($a$) The red joint and the green joints are inside the surface and outside the surface, respectively. ($b$) As the joints move, the red joint gets out of the surface and the green joints get in the surface. Furthermore, in order to keep the total GW length constant, new joints are created at the right (open circles).

Note that Equation 12 can be rewritten in form
$s(i,N)=ksλNcin˜ℓ.(xℓ-rℓ)nℓR(i,N)=ksλ2Ncinℓn˜ℓ,$
(15)
where $Nci$ represent the number of joints in contact with the surface for $i≤ℓ≤N$ and the angle brackets are averages taken over these contact joints. Since for a large $N$ (= total number of joints) $λ∝1/N$, $Nci∝N$, and the averages are nearly constants, it follows that in a first approximation $s(i,N)$ is constant and $R(i,N)∝λ$. Hence, the minimum of the surface energy has a linear variation in $λ$.
Because $ΔUsurf=Usurf(λ')-Usurf(λ)$ is equal to zero for $λ'=λ$, the lowest order of the energy variation can be written as
$ΔUsurf=1-λ'λϒ=1-1ηϒ,$
(16)
where
$ϒ=-λ'dUsurfdλ'|λ'=λ$
(16′)
is a constant.

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 $λ$ varies. The surface energy is a continuous function of $λ$, but as one point enters in or leaves the surface, the derivative of the surface energy becomes discontinuous. Hence, the position $xℓ$ 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 $λ$ 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

The total energy is the sum of Equations 3 and 13
$U(ɛi)=U(0)-q˜i.ɛi+12ɛ˜i.Pii.ɛi,$
(17)
with $qi=bi-s(i,N)$ and $Pii=Ai+R(i,N)$. The problem to be solved is then to find the vector $ɛi$ which minimizes $U$. In particular, for $Pii=0$ Equation 17 reduces to Equation 10 of Baier et al. (2016).
Let $ei,1$, $ei,2$, $ei,3$ be three orthonormal vectors with $ei,1=ui$. In this base, the vector $ɛi$ can be written in the form
$ɛi=(vi,1-1)ei,1+vi,2ei,2+vi,3ei,3=vi-ei,1,$
(18)
so that the new vector after the energy minimization is given by
$(ui)new=ui+ɛi=vi.$
(19)
Since the modulus of $(ui)new$ is equal to 1, the vector components $vi,1$, $vi,2$, $vi,3$ are not independent
$vi,12+vi,22+vi,32=1.$
(20)
First, considering the case $Pii=0$, the minimization of Equation 17 taking into account Equation 20 can be performed with a Lagrange multiplier $σi$
$f(vi,1,vi,2,vi,3,σi)=U(0)-qi,1(vi,1-1)+qi,2vi,2+qi,3vi,3+12σivi,12+vi,22+vi,32-1,$
(21)
with $qi,ℓ=q˜i.eiℓ$. Using standard procedures, the minimization gives
$vi,1ei,1+vi,2ei,2+vi,3ei,3=qi,1ei,1+qi,2ei,2+qi,3ei,3σiσi=qi,12+qi,22+qi,32.$
(22)
Thus, the new vector is equal to $q^i=qi/|qi|$. The result is analogous to a dipole moment $m$ in a magnetic field $B$, where the magnetic energy $-B˜.m$ is minimized when $m$ is parallel to $B$.
If the correction is small, the variation $(ui)new-ui$ is equal to the component of $q^i$ perpendicular to $ui$ in accordance with Baier et al. (2016). Nonetheless, in a frontal collision between the GW and artery $|s(i,N)|$ can be bigger than $|bi|$, and the vector $q^i$ eventually points in a direction nearly antiparallel to $ui$. Thus, the change of the unit vector will be very large and in order to fix this problem, the following formula is used instead of Equation 19
$(ui)new=(1-α)ui+αvi|(1-α)ui+αvi|,$
(23)
where $0<α<1$. Ideally $α→1$ and the criterion to choose $α$ is Equation 18 of Baier et al. (2016)
$|(ui)new-ui|≤ΔumaxN+1-i,$
(24)
with $Δumax∼0.1$.
If $Pii$ is different from zero, then Equation 21 is replaced by
$f(vi,1,vi,2,vi,3,σi)=U(0)-qi,1(vi,1-1)+qi,2vi,2+qi,3vi,3+12Pii,11(vi,1-1)2+Pii,22vi,22+Pii,33vi,32+Pii,12(vi,1-1)vi,2+Pii,13(vi,1-1)vi,3+Pii,23vi,2vi,3+12σivi,12+vi,22+vi,32-1,$
(25)
where $Pii,ℓm=e˜ℓ.Pii.em$ are the elements of a symmetric matrix. Hence, the equations to be solved are 20 and
$σi+Pii,11Pii,12Pii,13Pii,12σi+Pii,22Pii,23Pii,13Pii,23σi+Pii,33vi,1vi,2vi,3=qi,1+Pii,11qi,2+Pii,12qi,3+Pii,13.$
(26)
The solution $vi,1$, $vi,2$, $vi,3$, $σi$ can be found, for example, using the Newton--Raphson method for system of nonlinear equations (Press, Teukolsky, Vetterling, & Flannery, 1992). In one iteration step, the corrections $δvi,1$, $δvi,2$, $δvi,3$, $δσi$ are calculated by solving the equation
$σi+Pii,11Pii,12Pii,13vi,1Pii,12σi+Pii,22Pii,23vi,2Pii,13Pii,23σi+Pii,33vi,3vi,1vi,2vi,30δvi,1δvi,2δvi,3δσi=-Fi,1Fi,2Fi,3Fi,σ,$
(27)
where
$Fi,1=(σi+Pii,11)vi,1+Pii,12vi,2+Pii,13vi,3-(qi,1+Pii,11)Fi,2=Pii,12vi,1+(σi+Pii,22)vi,2+Pii,23vi,3-(qi,2+Pii,12)Fi,3=Pii,13vi,1+Pii,23vi,2+(σi+Pii,33)vi,3-(qi,3+Pii,13)Fi,σ=vi,12+vi,22+vi,32-12.$
(28)
After each step, it is convenient to renormalize the vector $vi,1ei,1+vi,2ei,2+vi,3ei,3$ and so $Fi,σ$ becomes equal to zero. As an initial guess for $vi,1$, $vi,2$, $vi,3$, and $σi$, the values in Equation 22 can be used; that is, the first-order solution is used as an starting point to find the second-order solution. The convergence is very fast because, at least for large $i$, the matrix elements of $Pii$ are small compared with the components of the vector $qi$.

### 2.3 Multiple Energy Minimizations

Consider the situation where $ui$ and $uj$ are varied simultaneously, with $j>i+1$. Then the equations are the same as in the previous section with the exception of the surface energy. In this case, Equation 13 is replaced by
$Usurf(ɛi,ɛj)=Usurf(0,0)+s˜(i,j-1).ɛi+12ɛ˜i.R(i,j-1).ɛi+s˜(j,N).(ɛi+ɛj)+12(ɛ˜i+ɛ˜j).R(j,N).(ɛi+ɛj)=Usurf(0,0)+s˜(i,N).ɛi+s˜(j,N).ɛj+12ɛ˜i.R(i,N).ɛi+12ɛ˜j.R(j,N).ɛj+ɛ˜i.R(j,N).ɛj.$
(29)
Defining the symmetric square matrices
$Pii=Ai+R(i,N)Pjj=Aj+R(j,N)Pij=R(j,N),$
(30)
the minimization proceeds as before. There are two restrictions $|vi|2=1$, $|vj|2=1$, and two Lagrange multipliers $σi$, $σj$ must be introduced. Equation 25 is replaced by
$f(vi,vj,σi,σj)=U(0,0)-q˜i.(vi-ei,1)-q˜j.(vj-ej,1)+12(v˜i-e˜i,1).Pii.(vi-ei,1)+12(v˜j-e˜j,1).Pjj.(vj-ej,1)+(v˜i-e˜i,1).Pij.(vj-ej,1)+12σi(v˜i.vi-1)+12σj(v˜j.vj-1),$
(31)
and Equation 26 is replaced by
$σiI+PiiPijPijσjI+Pjjvivj=qiqj+PiiPijPijPjjei,1ej,1.$
(32)
Furthermore, the Newton--Raphson method becomes
$σiI+PiiPijvi0PijσjI+Pjj0vjv˜i0˜000˜v˜j00δviδvjδσiδσj=-FiFjFi,σFj,σ,$
(33)
where
$FiFj=σiI+PiiPijPijσjI+Pjjvivj-qiqj-PiiPijPijPjjei,1ej,1.$
(34)
Note that $I$ represents the $3×3$ identity matrix and $0$ ($0˜$) is the zero column (row) vector of dimension three.

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×8$ in Equation 33 is replaced by a matrix of dimension $4M×4M$.

#### 2.3.1 Separation between Relaxing Segments

Omitting the constant terms, the relevant part of the energy in a double relaxation can be written as (see Equations 3 and 29)
$Uij=-b˜i.ɛi+12ɛ˜i.Ai.ɛi+s˜(i,N).ɛi+12ɛ˜i.R(i,N).ɛi-b˜j.ɛj+12ɛ˜j.Aj.ɛj+s˜(j,N).ɛj+12ɛ˜j.R(j,N).ɛj+ɛ˜i.R(j,N).ɛj,$
(35)
where $ɛi$ and $ɛj$ satisfy the constrains $|ei,1+ɛi|=|ej,1+ɛj|=1$. The first and third lines of Equation 35 represent the bending energies, the second and fourth lines represent the surface energies, and the last line represents the interaction between segments $i$ and $j$.

Consider first the linear approximation. As shown in the discussion following Equation 22, $-b˜i.ɛi$ will push the vector $(ui)new=ei,1+ɛi$ in a direction parallel to $bi$. Normally $κi-1∼κi$ because the angle between $ui-1$ and $ui$ is not far from the angle between $ui$ and $ui+1$. Moreover, $-1≤κi≤1$ and so it is concluded that $7-κi-1≈7-κ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˜(i,N).ɛi$ 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
$Uij=-b˜i.ɛi+12ɛ˜i.Ai.ɛi-b˜j.ɛj+12ɛ˜j.Aj.ɛj+s˜(j,N).(ɛi+ɛj)+12(ɛ˜i+ɛ˜j).R(j,N).(ɛi+ɛj).$
(36)
Since matrix $R(j,N)$ is positive definite, the interaction term between segment $i$ and segment $j$ will clearly help to decrease the energy if $ɛi+ɛj$ is small, or $ɛi≈-ɛj$. In other words, the corrections to the vectors $ui$ and $uj$ must be opposite and so the position of the joints $xk$ does not vary for $k≥j$.
When the difference between $i$ and $j$ is small and the GW is not far from equilibrium, the bending at joint $i$ and at joint $j$ will be similar, so that $bi≈bj$ and $Ai≈Aj$. Due to the symmetry dependence of the energy on $ɛi$ and on $ɛj$, the minimum will be at $ɛi≈ɛj$. Substituting $bi=bj$, $Ai=Aj$, and $ɛi=ɛj$ in Equation 36 and comparing with the single energy variation at $i$, it can be verified that in each term $ɛi(single)→2ɛi$, with the exception of $12ɛ˜i(single).Ai.ɛi(single)$ which is replaced by
$212ɛ˜i.Ai.ɛi=122ɛ˜i.Ai.2ɛi-142ɛ˜i.Ai.2ɛi.$
Thus, the energy $Ui(single)$ is identical to the energy of the double relaxation $Uij$ without the term $-142ɛ˜i.Ai.2ɛi$. If the minimum of $Ui(single)$ is at $ɛi(single)min$, then the energy difference between the minimum of $Ui(single)$ and the minimum of $Uij$ is at least $14ɛ˜i(single)min.Ai.ɛi(single)min$. It will be shown that the quadratic surface energy $12ɛ˜i.R(i,N).ɛi$ is much bigger than the quadratic bending energy $12ɛ˜i.Ai.ɛi$, so that the overall energy minimization is only marginal.

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 $ɛi=ɛi(single)min/2$ and $ɛj=ɛi(single)min/2$ than to correct once with $ɛi(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 $ɛi¬≈ɛj$ and the interaction term $ɛ˜i.R(j,N).ɛj$ can be negative, decreasing the energy even more.

The cases when $i$, $j$ are 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.
Figure 3.

GW subjected to bending forces (curved arrows) and surface forces $-si$ (red straight arrows). For simplicity, the artery is not shown and the arrows indicate the directions where the forces will displace the GW. ($a$) The two segments at points $A$ and $B$ are close and ($b$) the two segments at $D$ and $E$ are far from each other. A better choice for a double relaxation would be $A$, $C$ in case ($a$) and $D$, $F$ in case ($b$).

Figure 3.

GW subjected to bending forces (curved arrows) and surface forces $-si$ (red straight arrows). For simplicity, the artery is not shown and the arrows indicate the directions where the forces will displace the GW. ($a$) The two segments at points $A$ and $B$ are close and ($b$) the two segments at $D$ and $E$ are far from each other. A better choice for a double relaxation would be $A$, $C$ in case ($a$) and $D$, $F$ in case ($b$).

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
$ɛ˜i.R(i,N).ɛi=ɛ˜i.R(j,N).ɛi+∑ℓ∈Lijksλ2n˜ℓ.ɛi2≫ɛ˜i.R(j,N).ɛj,$
where the sum extends only over the set $Lij$ of joints in the interval $[i,j)$ which are in contact with the surface. Hence, the contribution of the interaction is only marginal compared to the single relaxation problem of segment $i$.

Figure 3b illustrates physically why multiple relaxation does not improve the outcomes when $j≫i$. 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 $δm$ ($m=1,…,M-1$) between the VS, different relaxation speeds are obtained. The goal is to optimize $M$ and $δm$, so that the relaxation is as fast as possible.

Let $i$, $j$, $k$, $…$ 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+δ1$, $k=j+δ2$, $…$ . 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≲δ1≲20$. Here, a recipe is given to choose $δ1$ dynamically and the procedure is extended to higher order multiple relaxations.

The best way to determine the optimal value of $δ1$ would be to calculate the energy variations for $δ1=5,6,…,25$ and to select the value which gives the lowest value. However, this procedure is highly inefficient. Instead of calculating $Uij$ in Equation 35 for $j=i+δ1$, only an estimation will be found. Considering the energy up to first order, the change in the unit vectors are $ɛi,1=q^i-ui$ and $ɛj,1=q^j-uj$. Assuming these values for $ɛi$ and $ɛj$, the variation of the energy when segment $j$ is included, is up to second order
$ΔUij=-q˜j.ɛj,1+12ɛ˜j,1.Pjj.ɛj,1+ɛ˜i,1.R(j,N).ɛj,1.$
(37)
The terms containing only $ɛi,1$ have been dropped, since they do not vary when different values of $j$ are taken into account. Observe that $qj=bj-sj$, $Pjj=Aj+R(j,N)$ are required in Equation 37, and that the surface terms $sj$, $R(j,N)$ were previously obtained during the calculation of $si$, $R(i,N)$. In order to determine $ΔUij$, the bending terms $bj$, $Aj$ still have to be computed, which is straightforward. Algorithm 1 summarizes the estimation of the energy variation.

Then the value of $j=i+5,i+6,…,i+25$ is selected, which gives the smallest $ΔUij$. When $Nlast the comparison is carried out only for $j=i+5,i+6,…,Nlast$. Further, if $Nlast 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≥i+5$ and Algorithm 2 finds a reasonable segment index.

The extension to a triple relaxation is straightforward. Suppose that $j$ has already been determined as described above and let $k=j+5,j+6,…,j+25$. Then calculate
$ΔUijk=-q˜k.ɛk,1+12ɛ˜k,1.Pkk.ɛk,1+(ɛ˜i,1+ɛ˜j,1).R(k,N).ɛk,1$
(38)
and select the value of $k$ which gives the smallest $ΔUijk$. In order to apply Algorithm 2, replace $j→i$ at the input and use the updated $ɛi,1$ which has been obtained in the previous function call. As before, if the return value is false, do not perform a triple relaxation but only a double relaxation.

## 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 $μ$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 $μ$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 $λ$ 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≥h$, so that $di$ is calculated every time an update of a segment with index $h∈[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∈(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 $δm$ is analyzed in simple cases. Then the algorithm to determine $δ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.

Figure 4.

Initial GW (green) and final GW at equilibrium (blue) inside an artery of diameter 2 mm. There are $N=200$ segments of length $λ=1$ mm and the positions of the joints with indices 50, 100, 165 are indicated by arrows.

Figure 4.

Initial GW (green) and final GW at equilibrium (blue) inside an artery of diameter 2 mm. There are $N=200$ segments of length $λ=1$ mm and the positions of the joints with indices 50, 100, 165 are indicated by arrows.

### 3.1 Parameter Optimization

A huge amount of data has been analyzed, but only a few selected tests are presented. As an example, consider Figure 4 which shows the initial GW and the GW at equilibrium. A representative joint index $i=50$ was selected and the time evolution of the parameter $d50$ is shown in Figure 5. The calculations were performed for a single relaxation ($M=1$), for a double relaxation ($M=2$), and for a triple relaxation ($M=3$).
Figure 5.

Lateral deviation of joint index $i=50$ in Figure 4 as function of the number of cycles up to 100. The results are shown for different numbers $M$ of VS and different separations $δ1$, $δ2$ between the VS.

Figure 5.

Lateral deviation of joint index $i=50$ in Figure 4 as function of the number of cycles up to 100. The results are shown for different numbers $M$ of VS and different separations $δ1$, $δ2$ between the VS.

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≲δ1≲20$. On the other hand, for $2≤δ1≲5$ and for $δ1≳30$ the convergence is slow, but still faster than with $M=1$. Hence, there is no substantial improvement in using multiple relaxation if $δ1$ is very short or very large. Furthermore, choosing $M=3$, $δ1=δ2$ (Figure 5b) results in faster relaxations than with $M=2$, $δ1$ (Figure 5a) and in both cases the best outcomes are obtained for $10≲δ1≲20$. However, note that for a smaller segment size $λ$ and a larger number of segments $N$, the optimal $δ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 $δ1$ and $δ2$ chosen dynamically will be explored.

### 3.2 Algorithm Evaluation

As an example of the application of the algorithm, consider again the relaxation of the GW shown in Figure 4. The GW has 200 segments (which is indeed a very large number) and the information of one cycle is summarized in Table 1. From the measurements, it follows that as the number of blocks increases, the time to solve the equations becomes dominant over the time to compute the interactions with the surface. Hence, a very large $M$ should not be used.
Table 1.
Time (in Seconds) to Compute One Cycle of the GW-Surface Interaction (Equation 12) and to Obtain and Solve Equation 26
$M$123
surface 0.200 0.221 0.229
solution 0.201 0.394 0.777
$M$123
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×106$ times and the solution is performed $20.1×103$ times. The numerical experiments have been carried out with an Intel Core i7-4500U (1.8 GHz).

Figure 6.

Lateral deviation of joint index $i=50$ in Figure 4 as function of the number of cycles (up to 100). Each relaxation curve has a total of $132×103$ points. The results are shown for different number of VS.

Figure 6.

Lateral deviation of joint index $i=50$ in Figure 4 as function of the number of cycles (up to 100). Each relaxation curve has a total of $132×103$ points. The results are shown for different number of VS.

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 $δ1=15$ are indeed faster at the beginning, but the relaxations with the proposed algorithm are more stable.

A small segment index ($i≪N$) 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.

Figure 7.

Lateral deviation of joint index $i=100$ in Figure 4 as function of the number of cycles. Each relaxation curve has $515×103$ points and a thicker curve indicates that the joint oscillates with bigger amplitude (i.e., a noise). The results are shown for different number of VS.

Figure 7.

Lateral deviation of joint index $i=100$ in Figure 4 as function of the number of cycles. Each relaxation curve has $515×103$ points and a thicker curve indicates that the joint oscillates with bigger amplitude (i.e., a noise). The results are shown for different number of VS.

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$.

The parameter $d100$ of the same GW relaxation is shown in Figure 7. This corresponds to an intermediate segment index and there are less contact points opposing the rigid translation of the GW. Hence, all relaxations are faster than in Figure 6 and the results with $M=1$, 2, and 3 give similar relaxation times. Nonetheless, the single relaxation is oscillatory because the slow relaxation of segments with indices $h≪i$, gradually move the joints $h. Segment $i$ must rotate in one direction or another to continuously accommodate the rest of the GW for $h∈[i,N]$. Furthermore, the thickness of the line indicates that the joint is oscillating around the equilibrium position, with an oscillation amplitude typically within the range between $10-4$ mm and $10-2$ mm. This is a very small displacement and cannot be observed at all in simulations, since the point remains in the same pixel during the rendering. At about 50 cycles the noise in the lateral deviation corresponding to $M=2$ and to $M=3$ is very small, because most segments with indices $h are already close to the equilibrium.
Figure 8.

Lateral deviation of joint index $i=165$ in Figure 4 as function of the number of cycles up to 10. Each relaxation curve has $139×103$ points. The green ($M=2$) and the blue ($M=3$) curves show a very sharp decrease in the beginning.

Figure 8.

Lateral deviation of joint index $i=165$ in Figure 4 as function of the number of cycles up to 10. Each relaxation curve has $139×103$ points. The green ($M=2$) and the blue ($M=3$) curves show a very sharp decrease in the beginning.

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 $λ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).

Figure 9.

Relative frequency (in percentage) of the number of times the algorithm chooses $δ1$, $δ2$, or $δ3$ in the interval [5, 25]. The results correspond to the GW relaxations with $M=2$, $M=3$, and $M=4$. The last one does not appear in Figures 68, but it is quite similar to $M=3$.

Figure 9.

Relative frequency (in percentage) of the number of times the algorithm chooses $δ1$, $δ2$, or $δ3$ in the interval [5, 25]. The results correspond to the GW relaxations with $M=2$, $M=3$, and $M=4$. The last one does not appear in Figures 68, but it is quite similar to $M=3$.

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 $δ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 $δ1=5$ is somewhat higher than the frequency of $δ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 $δ2=5$ is large. After choosing the best value for $δ1=j-i$, the second $δ2=k-j$ will usually not be large because the probability decreases with the separation between segment $i$ and segment $k$.

As a final test consider the GW in Figure 15 of Baier et al. (2016). Figure 10 shows the relaxation of joint index $i=92$ which is located in the middle of the gray box. Besides the single, double, and triple relaxation, the geometrical relaxation is also included. The triple is the fastest one and the double also shows a sharp decrease within the first cycle. The relaxation method developed by Konings et al. (2003) corresponds essentially to the red curve in Figure 10, although the recipe for the single energy minimization in the present work (Equation 22) is by far more efficient. Moreover, comparisons with FEM and MSM are not given as they are known to result in longer computation times (Sharei et al., 2018).
Figure 10.

Lateral deviation of joint index $i=92$ in Figure 15 of Baier, Baier-Saip, Schilling, and Oliveira (2016) up to 80 cycles. The results are shown for different numbers of VS. For comparison, the geometrical method (black curve) has been included in a single relaxation.

Figure 10.

Lateral deviation of joint index $i=92$ in Figure 15 of Baier, Baier-Saip, Schilling, and Oliveira (2016) up to 80 cycles. The results are shown for different numbers of VS. For comparison, the geometrical method (black curve) has been included in a single relaxation.

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 $ɛi=(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 $|ɛi|$ 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.

Figure 11.

The red bars represent the percentage of the time in which each joint (blue curve in Figure 4) is in contact with the artery's surface. The black curve represents the average modulus of the vector update $ɛi$. ($a$) The calculations include the second order terms $Ai$ and $R(i,N)$. ($b$) The calculations do not include $Ai$ and $R(i,N)$; that is, $Pii$ vanishes in Equation 17.

Figure 11.

The red bars represent the percentage of the time in which each joint (blue curve in Figure 4) is in contact with the artery's surface. The black curve represents the average modulus of the vector update $ɛi$. ($a$) The calculations include the second order terms $Ai$ and $R(i,N)$. ($b$) The calculations do not include $Ai$ and $R(i,N)$; that is, $Pii$ vanishes in Equation 17.

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 $|ɛi|$ 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$.

Consider a simplified one dimensional version of the surface energy in Equation 9
$U2nd=12ksx2=12ksx02+2x0δx+(δx)2forx>00forx<0,$
(39a)
with $x0=n˜ℓ.(xℓ-rℓ)$ and $δx=λn˜ℓ.ɛi$. If the term $(δx)2$ is dropped, the first-order approximation results in
$U1st=12ksx02+2x0δx+12ksx02=ksx0xforx>00forx<0,$
(39b)
where the “arbitrary constant” $12ksx02$ has been included in the first line so as to make the energy $U1st$ continuous at $x=0$. Equations 39a and 39b are plotted in Figure 12a. In particular, the derivative of $U2nd$ is continuous at $x=0$, but the derivative of $U1st$ (which depends on $x0$) is discontinuous at $x=0$. Hence, in the second-order approximation, a class $C1$ curve is obtained; however, in the first order approximation, a class $C0$ curve is obtained. Moreover, the derivative of $U2nd$ at $x=x0$ is equal to the derivative of $U1st$.
Figure 12.

Energy when the joint is inside the surface $x>0$ (red, green, blue curves) and outside the surface $x<0$ (magenta curve). ($a$) Surface energy $Usurf$ in the second order approximation $U2nd$, in the first order approximation $U1st,a$ for a small $x0=xa$ and $U1st,b$ for a large $x0=xb$. The derivatives of $U2nd$ and $U1st,a$ coincide at $xa$, and the derivatives of $U2nd$ and $U1st,b$ coincide at $xb$. Performing a numerical minimization of the total energy $Ubend+Usurf$ in the region $x>0$ results in ($b$) a minimum at $x2,min>0$, ($c$) a minimum at $x1a,min>0$, and ($d$) a minimum at $x1b,min<0$. The numerical minimization of the energy $Ubend$ in the region $x<0$ results in a minimum at $xbend,min>0$.

Figure 12.

Energy when the joint is inside the surface $x>0$ (red, green, blue curves) and outside the surface $x<0$ (magenta curve). ($a$) Surface energy $Usurf$ in the second order approximation $U2nd$, in the first order approximation $U1st,a$ for a small $x0=xa$ and $U1st,b$ for a large $x0=xb$. The derivatives of $U2nd$ and $U1st,a$ coincide at $xa$, and the derivatives of $U2nd$ and $U1st,b$ coincide at $xb$. Performing a numerical minimization of the total energy $Ubend+Usurf$ in the region $x>0$ results in ($b$) a minimum at $x2,min>0$, ($c$) a minimum at $x1a,min>0$, and ($d$) a minimum at $x1b,min<0$. The numerical minimization of the energy $Ubend$ in the region $x<0$ results in a minimum at $xbend,min>0$.

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 $ɛi$ becomes smaller and the stability increases. No numerical artifact was introduced to bound $|ɛi|$, 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 $λ$ affects the results of the calculations. Note that $λ$ 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.

Figure 13 shows that the energy varies less than 2.5% as the GW resolution $η=N'/N$ increases by a factor 3. Each dot represents the calculation of $ΔU/U$ for $N'$ segments, with $175≤N'≤525$. According to Equations 8 and 16, the total energy variation as the GW depends on the resolution as follows
$ΔUU=c11-1η+c21-1η2,$
(40)
where the constant $c1$ is associated with the surface energy and the constant $c2$ is associated with the bending energy. Fitting the points in Figure 13 gives $c1=0.0263$ and $c2=0.00554$, showing that the model is more sensitive to variations of the GW-surface interaction than to variations of the GW deformation. In particular, note that as the resolution $η$ increases, $ΔU/U$ approaches to the constant value $c1+c2$ indicating that the numerical solution converges, as expected when a problem is solved right.

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 $ΔU/U$ has been estimated to be $∼3×10-5$. The same accuracy is obtained when the GW is perturbed and after that the equilibrium is found again. The change in $ΔU/U$ observed in Figure 13 between consecutive values of $N'$ is on average $92×10-5$ with a maximum of $387×10-5$, which is much bigger than the accuracy of the calculations.

Figure 13.

($a$) GW with $N=175$ segments. This number is increased to $N'$ and $λ=1$ mm is decreased to $λ'$ in such a way that $Nλ=N'λ'$; that is, the total GW length remains constant. ($b$) Relative energy variation $ΔU/U$ as function of the GW resolution $η=λ/λ'=N'/N$ for $N'∈[175,525]$ (red dots). The blue curve represents the fitting corresponding to Equation 40.

Figure 13.

($a$) GW with $N=175$ segments. This number is increased to $N'$ and $λ=1$ mm is decreased to $λ'$ in such a way that $Nλ=N'λ'$; that is, the total GW length remains constant. ($b$) Relative energy variation $ΔU/U$ as function of the GW resolution $η=λ/λ'=N'/N$ for $N'∈[175,525]$ (red dots). The blue curve represents the fitting corresponding to Equation 40.

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 $η$ 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 $θi$ between consecutive segments is also large and the error in the bending energy proportional to $θ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.

Figure 14.

A few typical examples showing how $Ks$ (see Equation 41) varies along the surface of the artery. In the Monte Carlo simulation, a set of 5000 random curves have been generated to perform the calculations, and for each one the equilibrium of the GW was found.

Figure 14.

A few typical examples showing how $Ks$ (see Equation 41) varies along the surface of the artery. In the Monte Carlo simulation, a set of 5000 random curves have been generated to perform the calculations, and for each one the equilibrium of the GW was found.

Instead of a constant $ks$, the effective spring stiffness of the artery segments with index $ω=1,2,…,Ω$ is replaced by
$Ks(ω)=ks+∑g=1GksgcosgπωΩ+ϕg,$
(41)

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

A few typical examples of $Ks$ are shown in Figure 14. Since $1≤ω≤Ω$ only half of the period is used and $Ks(ω)$ 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 $Ω/2$ as an arbitrary set of numbers ${Ks}$ containing $Ω$ elements can be represented by a Fourier series with the first $Ω/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).

Figure 15.

The PDF of $Ks$ (see Equation 41). The total number of points is $5000×162$ and $96.6%$ of the points are in the interval $(0.5ks,1.5ks)$. Moreover, the mean is equal to $1.0012ks$ and the standard deviation is equal to $0.236ks$.

Figure 15.

The PDF of $Ks$ (see Equation 41). The total number of points is $5000×162$ and $96.6%$ of the points are in the interval $(0.5ks,1.5ks)$. Moreover, the mean is equal to $1.0012ks$ and the standard deviation is equal to $0.236ks$.

For a given $Ks(ω)$ 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(ω)$ (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.

The final result for each GW segment is detailed in Figure 17. It can be seen that the modulus of the mean is mostly smaller than the rms, showing that the reference GW is inside the range of positions calculated for the variable GWs. Furthermore, Figure 18 shows the resulting PDF of the lateral deviation in a Monte Carlo simulation. When the elastic constant $Ks$ is in the range $±50%$ of $ks$ (Figure 15), the lateral deviation varies by about $±0.005$ mm or less, which is negligible for the simulator. Hence, there is no significant difference between the calculations performed with constant $ks$ and with variable $Ks(ω)$.
Figure 16.

Monte Carlo simulation with variable $Ks$. The green curve is the mean deviation between a segment of the variable GW and the reference GW as function of the number of steps. The red curve is the rms value of the deviation. The top figure corresponds to the GW segment 30, the middle figure corresponds to the GW segment 100, and the bottom figure corresponds to the GW segment 175 (the tip).

Figure 16.

Monte Carlo simulation with variable $Ks$. The green curve is the mean deviation between a segment of the variable GW and the reference GW as function of the number of steps. The red curve is the rms value of the deviation. The top figure corresponds to the GW segment 30, the middle figure corresponds to the GW segment 100, and the bottom figure corresponds to the GW segment 175 (the tip).

Figure 17.

Monte Carlo simulation with variable $Ks$ after 5000 steps. The green curve represents the mean deviation between a segment of the variable GW and the reference GW as function of the segment number. The red curve represents the corresponding rms value.

Figure 17.

Monte Carlo simulation with variable $Ks$ after 5000 steps. The green curve represents the mean deviation between a segment of the variable GW and the reference GW as function of the segment number. The red curve represents the corresponding rms value.

Figure 18.

The PDF of the lateral deviation $d$ when the elastic constant $Ks$ is varied in a Monte Carlo simulation.

Figure 18.

The PDF of the lateral deviation $d$ when the elastic constant $Ks$ is varied in a Monte Carlo 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×4M$, 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 $δm$ between the VS. In the simplest case, it has been found that the best outcomes are obtained for $10≲δ1≲20$. For smaller values of $δ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 $δ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 $δ1$ increases. An algorithm to choose dynamically $δ1,δ2,…,δ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).

## REFERENCES

Alderliesten
,
T.
,
Bosman
,
P. A. N.
, &
Niessen
,
W. J.
(
2006
).
Towards a real-time minimally-invasive vascular intervention simulation system
.
IEEE Transactions on Medical Imaging
,
26
,
128
132
.
Alderliesten
,
T.
,
Konings
,
M. K.
, &
Niessen
,
W. J.
(
2004
).
Simulation of minimally invasive vascular interventions for training purposes
.
Computer Aided Surgery
,
9
,
3
15
.
Alderliesten
,
T.
,
Konings
,
M. K.
, &
Niessen
,
W. J.
(
2007
).
Modeling friction, intrinsic curvature, and rotation of guide wires for simulation of minimally invasive vascular interventions
.
IEEE Transactions on Biomedical Engineering
,
54
,
29
38
.
Al-Moghairi
,
A. M.
, &
Al-Amri
,
H. S.
(
2013
).
Management of retained intervention guide-wire: A literature review
.
Current Cardiology Reviews
,
9
,
260
266
.
Anderson
,
J. H.
,
Raghavan
,
R.
,
Wang
,
Y. P.
,
Mullick
,
R.
, &
Kong
,
C. C.
(
1997
).
DaVinci—A vascular catheterization simulator
.
Journal of Vascular and Interventional Radiology
,
8
,
261
.
Antman
,
S. S.
(
2005
).
Nonlinear problems of elasticity
(2nd ed.).
New York
:
Springer
.
Baier-Saip
,
J. A.
,
Baier
,
P. A.
,
Schilling
,
K.
, &
Oliveira
,
J. C.
(
2017
).
Approximate artery elasticity using linear springs
.
Journal of Medical and Biological Engineering
,
37
,
899
911
.
Baier
,
P. A.
,
Baier-Saip
,
J. A.
,
Schilling
,
K.
, &
Oliveira
,
J. C.
(
2016
).
Simulator for minimally invasive vascular interventions: Hardware and software
.
Presence: Teleoperators and Virtual Environments
,
25
,
108
128
.
Baier
,
P. A.
,
Srinivasan
,
L.
,
Baier-Saip
,
J. A.
,
Voelker
,
W.
, &
Schilling
,
K.
(
2015
).
Surfaces for modeling arteries in virtual reality simulators
.
IFAC-PapersOnLine
,
48
,
031
036
.
Balaji
,
N. R.
, &
Shah
,
P. B.
(
2011
).
.
Circulation
,
124
,
e407
e408
.
Basdogan
,
C. D. S.
,
Kim
,
J.
,
Muniyandi
,
M.
, &
Srinivasan
,
M. A.
(
2004
).
Haptics in minimally invasive surgical simulation and training
.
IEEE Computer Graphics and Applications
,
24
,
56
64
.
Beck
,
J. L.
(
2010
).
Bayesian system identification based on probability logic
.
Structural Control Health Monitoring
,
17
,
825
847
.
Beck
,
J. L.
, &
Katafygiotis
,
L. S.
(
1998
).
Updating models and their uncertainties. I: Bayesian statistical framework
.
Journal of Engineering Mechanics
,
124
,
455
461
.
Bergou
,
M.
,
Wardetzky
,
M.
,
Robinson
,
S.
,
Audoly
,
B.
, &
Grinspun
,
E.
(
2008
).
Discrete elastic rods
.
ACM Transactions on Graphics
,
27
,
63:1
63:12
.
Bosman
,
P. A. N.
, &
Alderliesten
,
T.
(
2005
).
Evolutionary algorithms for medical simulations—A case study in minimally-invasive vascular interventions
. In
Genetic and Evolutionary Computation Conference Workshop
(pp.
125
132
).
Bui
,
H. P.
,
Tomar
,
S.
, &
Bordas
,
S. P. A.
(
2019
).
Corotational cut finite element method for real-time surgical simulation: Application to needle insertion simulation
.
Computer Methods in Applied Mechanics and Engineering
,
345
,
183
211
.
Bui
,
H. P.
,
Tomar
,
S.
,
Courtecuisse
,
H.
,
Audette
,
M.
,
Cotin
,
S.
, &
Bordas
,
S. P. A.
(
2018
).
Controlling the error on target motion through real-time mesh adaptation: Applications to deep brain stimulation
.
International Journal for Numerical Methods in Biomedical Engineering
,
34
,
e2958
.
Bui
,
H. P.
,
Tomar
,
S.
,
Courtecuisse
,
H.
,
Cotin
,
S.
, &
Bordas
,
S. P. A.
(
2018
).
Real-time error control for surgical simulation
.
IEEE Transaction on Biomedical Engineering
,
65
,
596
607
.
Cai
,
Y. Y.
,
Chui
,
C.-K.
,
Ye
,
X.
,
Anderson
,
J. H.
,
Liew
,
K.-M.
, &
Sakuma
,
I.
(
2004
).
Simulation-based virtual prototyping of customized catheterization devices
.
Journal of Computing and Information Science in Engineering
,
4
,
132
139
.
Cardoso
,
F. M.
, &
Furuie
,
S. S.
(
2016
).
Guidewire path determination for intravascular applications
.
Computer Methods in Biomechanics and Biomedical Engineering
,
19
,
628
638
.
Carter
,
F. J.
,
Schijven
,
M. P.
,
Aggarwal
,
R.
,
Grantcharov
,
T.
,
Francis
,
N. K.
,
Hanna
,
G. B.
, &
Jakimowicz
,
J. J.
(
2005
).
Consensus guidelines for validation of virtual reality surgical simulators
.
Surgical Endoscopy and Other Interventional Techniques
,
19
,
1523
1532
.
Cheng
,
X.-R.
,
Song
,
Q.-K.
,
Xie
,
X.-L.
,
Cheng
,
L.
,
Wang
,
L.
,
Bian
,
G.-B.
, …
Prasong
,
P.
(
2017
).
A fast and stable guidewire model for minimally invasive vascular surgery based on lagrange multipliers
. In
Seventh International Conference on Information Science and Technology
(pp.
109
114
).
Chui
,
C. K.
,
Nguyen
,
H. T.
,
Wang
,
Y. P.
,
Mullick
,
R.
,
Raghavan
,
R.
, &
Anderson
,
J. A.
(
1996
).
Potential field of vascular anatomy for realtime computation of catheter navigation
. In
R.
Banvard
(Ed.),
First Visible Human Conference
,
Bethesda, MD, USA
.
Coles
,
T. R.
,
Meglan
,
D.
, &
John
,
N.
(
2011
).
The role of haptics in medical training simulators: A survey of the state of the art
.
IEEE Transactions on Haptics
,
4
,
51
66
.
Cotin
,
S.
,
Dawson
,
S. L.
,
Meglan
,
D.
,
Shaffer
,
D. W.
,
Ferrell
,
M. A.
,
Bardsley
,
R. S.
, …
Wendlandt
,
J.
(
2000
).
ICTS, an interventional cardiology training system
. In
Medicine Meets Virtual Reality 2000
(
Vol. 70
, pp.
59
65
).
Studies in Health Technology and Informatics
.
Cotin
,
S.
,
Duriez
,
C.
,
Lenoir
,
J.
,
Neumann
,
P.
, &
Dawson
,
S.
(
2005
).
.
Computer Aided Surgery
,
8
,
534
542
.
Courtecuisse
,
H.
, &
Allard
,
J.
(
2009
).
Parallel dense Gauss-Seidel algorithm on many-core processors
. In
HPCC '09 Proceedings of the 2009 11th IEEE International Conference on High Performance Computing and Communications
(pp.
139
147
).
Courtecuisse
,
H.
,
Allard
,
J.
,
Kerfriden
,
P.
,
Bordas
,
S. P. A.
,
Cotin
,
S.
, &
Duriez
,
C.
(
2014
).
Real-time simulation of contact and cutting of heterogeneous soft-tissues
.
Medical Image Analysis
,
18
,
394
410
.
Dawson
,
S. L.
,
Cotin
,
S.
,
Meglan
,
D.
,
Shaffer
,
D. W.
, &
Ferrell
,
M. A.
(
2000
).
Designing a computer-based simulator for interventional cardiology training
.
Catheterization and Cardiovascular Interventions
,
51
,
522
527
.
Duriez
,
C.
,
Cotin
,
S.
,
Lenoir
,
J.
, &
Neumann
,
P.
(
2006
).
.
Computer Aided Surgery
,
11
,
300
308
.
Elster
,
C.
, &
Wübbeler
,
G.
(
2015
).
Bayesian regression versus application of least squares—An example
.
Metrologia
,
53
,
S10
S16
.
Engum
,
S. A.
,
Jeffries
,
P.
, &
Fisher
,
L.
(
2003
).
Intravenous catheter training system: Computer-based education versus tradition learning methods
.
The American Journal of Surgery
,
186
,
67
74
.
Ganji
,
Y.
, &
Janabi-Sharifi
,
F.
(
2009
).
.
IEEE Transactions on Biomedical Engineering
,
56
,
621
632
.
Grenon
,
S. M.
,
Reilly
,
L. M.
, &
Ramaiah
,
V. G.
(
2011
).
Technical endovascular highlights for crossing the difficult aortic bifurcation
.
Journal of Vascular Surgery
,
54
,
893
896
.
Hauseux
,
P.
,
Hale
,
J. S.
, &
Bordas
,
S. P. A.
(
2017a
).
Accelerating Monte Carlo estimation with derivatives of high-level finite element models
.
Computer Methods in Applied Mechanics and Engineering
,
318
,
917
936
.
Hauseux
,
P.
,
Hale
,
J. S.
, &
Bordas
,
S. P. A.
(
2017b
).
Calculating the Malliavin derivative of some stochastic mechanics problems
.
PLOS ONE
,
12
,
e0189994
.
Hauseux
,
P.
,
Hale
,
J. S.
,
Cotin
,
S.
, &
Bordas
,
S. P. A.
(
2018
).
Quantifying the uncertainty in a hyperelastic soft tissue model with stochastic parameters
.
Applied Mathematical Modelling
,
62
,
86
102
.
Konings
,
M. K.
,
van de Kraats
,
E. B.
,
Alderliesten
,
T.
, &
Niessen
,
W. J.
(
2003
).
Analytical guide wire motion algorithm for simulation of endovascular interventions
.
Medical & Biological Engineering & Computing
,
41
,
689
700
.
Koutsourelakis
,
P. S.
(
2012
).
A novel Bayesian strategy for the identification of spatially varying material properties and model validation: An application to static elastography
.
International Journal for Numerical Methods in Engineering
,
91
,
249
268
.
Lawton
,
W.
,
Raghavan
,
R.
,
Ranjan
,
S. R.
, &
Viswanathan
,
R. R.
(
2000
).
Tubes in tubes: Catheter navigation in blood vessels and its applications
.
International Journal of Solids and Structures
,
37
,
3031
3054
.
Lenoir
,
J.
,
Cotin
,
S.
,
Duriez
,
C.
, &
Neumann
,
P.
(
2006
).
Physics-based models for catheter, guidewire and stent simulation
.
Studies in Health Technology and Informatics
,
119
,
305
310
.
Lenoir
,
J.
,
Meseure
,
P.
,
Grisoni
,
L.
, &
Chaillou
,
C.
(
2002
).
. In
ESAIM: PROCEEDINGS
(
Vol. 12
, pp.
102
107
).
Li
,
S.
,
Guo
,
J.
,
Wang
,
Q.
,
Meng
,
Q.
,
Chui
,
Y.
,
Qin
,
J.
, &
Heng
,
P.
(
2012
).
A catheterization-training simulator based on a fast multigrid solver
.
IEEE Computer Graphics and Applications
,
32
,
56
70
.
Li
,
Z.
,
Chui
,
C. K.
,
Anderson
,
J. H.
,
Chen
,
X.
Ma
,
X.
,
Hua
,
W.
, …
Nowinski
,
W. L.
(
2001
).
Computer environment for interventional neuroradiology procedures
.
Simulation & Gaming
,
32
,
404
419
.
Lim
,
H. L.
,
Shetty
,
B. R.
,
Chui
,
C. K.
,
Wang
,
Y. P.
, &
Cai
,
Y. Y.
(
1998
).
Real-time interactive surgical simulator for catheter navigation
. In
Proceedings of SPIE
(
Vol. 3262
, pp.
4
14
).
Liu
,
A.
,
Tendick
,
F.
,
Cleary
,
K.
, &
Kaufmann
,
C.
(
2003
).
A survey of surgical simulation: Applications, technology, and education
.
Presence: Teleoperators and Virtual Environments
,
12
,
599
614
.
Lunderquist
,
A.
,
Ivancev
,
K.
,
Wallace
,
S.
,
Enge
,
I.
,
Laerum
,
F.
, &
Kolbenstvedt
,
A.
(
1995
).
The acquisition of skills in interventional radiology by supervised training on animal models: A three year multicenter experience
.
,
18
,
209
211
.
,
S.
,
Sista
,
B.
, &
Vemaganti
,
K.
(
2015
).
A Bayesian approach to selecting hyperelastic constitutive models of soft tissue
.
Computer Methods in Applied Mechanics and Engineering
,
291
,
102
122
.
Meier
,
U.
,
Lopez
,
O.
,
Monserrat
,
C.
,
Juan
,
M. C.
, &
Alcañiz
,
M.
(
2005
).
Real-time deformable models for surgery simulation: A survey
.
Computer Methods and Programs in Biomedicine
,
77
,
183
197
.
Mohamedou
,
M.
,
Zulueta
,
K.
,
Chung
,
C. N.
,
Rappel
,
H.
,
Beex
,
L.
,
,
L.
, …
Noels
,
L.
(
2019
).
Bayesian identification of Mean-Field Homogenization model parameters and uncertain matrix behavior in non-aligned short fiber composites
.
Composite Structures
,
220
,
64
80
.
Mori
,
T.
,
Hatano
,
N.
,
Maruyama
,
S.
, &
Atomi
,
Y.
(
1998
).
Significance of hands-on training in laparoscopic surgery
.
Surgical Endoscopy
,
12
,
256
260
.
Nowinski
,
W. L.
, &
Chui
,
C. K.
(
2001
).
. In
International Workshop on Medical Imaging and Augmented Reality
(pp.
87
94
).
Pai
,
D. K.
(
2002
).
STRANDS: Interactive simulation of thin solids using Cosserat models
.
Computer Graphics Forum
,
21
,
347
352
.
Payan
,
Y.
(Ed.) (
2012
).
Soft tissue biomechanical modeling for computer assisted surgery
.
Berlin
:
Springer
.
Pelletier
,
M. P.
,
Kaneko
,
T.
,
Peterson
,
M. D.
, &
Thourani
,
V. H.
(
2017
).
From satures to wires: The evolving necessities of cardiac surgery training
.
The Journal of Thoracic and Cardiovascular Surgery
,
154
,
990
993
.
Press
,
W. H.
,
Teukolsky
,
S. A.
,
Vetterling
,
W. T.
, &
Flannery
,
B. P.
(
1992
).
Numerical recipes in C
(2nd ed.).
Cambridge
:
Cambridge University Press
.
Rappel
,
H.
, &
Beex
,
L. A. A.
(
2019
).
Estimating fibres' material parameter distributions from limited data with the help of Bayesian inference
.
European Journal of Mechanics A/Solids
,
75
,
169
196
.
Rappel
,
H.
,
Beex
,
L. A. A.
, &
Bordas
,
S. P. A.
(
2018
).
Bayesian inference to identify parameters in viscoelasticity
.
Mechanics of Time-Dependent Materials
,
22
,
221
258
.
Rappel
,
H.
,
Beex
,
L. A. A.
,
Hale
,
J. S.
,
Noels
,
L.
, &
Bordas
,
S. P. A.
(
2019
).
A tutorial on Bayesian inference to identify material parameters in solid mechanics
.
Archives of Computational Methods in Engineering
.
https://doi.org/10.1007/s11831-018-09311-x
Rappel
,
H.
,
Beex
,
L. A. A.
,
Noels
,
L.
, &
Bordas
,
S. P. A.
(
2019
).
Identifying elastoplastic parameters with Bayes' theorem considering output error, input error and model uncertainty
.
Probabilistic Engineering Mechanics
,
55
,
28
41
.
Schneider
,
P. A.
(Ed.) (
2003
).
Endovascular skills: Guidewire and catheter skills for endovascular surgery
(2nd ed.).
New York
:
Marcel Dekker
.
Sharei
,
H.
,
Alderliesten
,
T.
,
Dobbelsteena
,
J. J.
, &
Dankelman
,
J.
(
2018
).
Navigation of guidewires and catheters in the body during intervention procedures: A review of computer-based models
.
Journal of Medical Imaging
,
5
,
010902
.
Spillmann
,
J.
, &
Harders
,
M.
(
2010
).
Inextensible elastic rods with torsional friction based on Lagrange multipliers
.
Computer Animation and Virtual Worlds
,
21
,
561
572
.
Takashima
,
K.
,
Shimomura
,
R.
,
Kitou
,
T.
,
,
H.
,
Yoshinaka
,
K.
, &
Ikeuchi
,
K.
(
2007
).
Contact and friction between catheter and blood vessel
.
Tribology International
,
40
,
319
328
.
Takashima
,
K.
,
Tsuzuki
,
S.
,
Ooike
,
A.
,
Yoshinaka
,
K.
,
Yu
,
K.
,
Ohta
,
M.
, &
Mori
,
K.
(
2014
).
Numerical analysis and experimental observation of guidewire motion in a blood vessel model
.
Medical Engineering & Physics
,
36
,
1672
1683
.
Tang
,
W.
,
,
P.
,
Gould
,
D.
,
Wan
,
T.
,
Zhai
,
J.
, &
How
,
T.
(
2010
).
A realistic elastic rod model for real-time simulation of minimally invasive vascular interventions
.
The Visual Computer
,
26
,
1157
1165
.
Tang
,
W.
,
Wan
,
T. R.
,
Gould
,
D. A.
,
How
,
T.
, &
John
,
N. W.
(
2012
).
A stable and real-time nonlinear elastic approach to simulating guidewire and catheter insertions based on Cosserat rod
.
IEEE Transactions on Biomedical Engineering
,
59
,
2211
2218
.
Tsang
,
J. S.
,
Naughton
,
P. A.
,
Leong
,
S.
,
Hill
,
A. D. K.
,
Kelly
,
C. J.
, &
Leahy
,
A. L.
(
2008
).
Virtual reality simulation in endovascular surgical training
.
The Surgeon
,
6
,
214
220
.
Wang
,
F.
,
Duratti
,
L.
,
Samur
,
E.
,
Spaelter
,
U.
, &
Bleuler
,
H.
(
2007
).
A computer-based real-time simulation of interventional radiology
. In
29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society
(pp.
1742
1745
).
Wang
,
Y.
,
Chui
,
C.
,
Lim
,
H.
,
Cai
,
Y.
, &
Mak
,
K.
(
1998
).
Real-time interactive simulator for percutaneous coronary revascularization procedures
.
Computer Aided Surgery
,
3
,
211
227
.
Willaert
,
W. I. M.
,
Aggarwal
,
R.
,
Herzeele
,
I.
,
Cheshire
,
N. J.
, &
Vermassen
,
F. E.
(
2012
).
Recent advancements in medical simulation: Patient-specific virtual reality simulation
.
World Journal of Surgery
,
36
,
1703
1712
.
Wu
,
J. X.
,
Chen
,
G. C.
,
Chang
,
C. W.
, &
Lin
,
C. H.
(
2016
).
Development of virtual-reality simulator system for minimally invasive surgery (MIS) using fractional-order vascular access
. In
2016 SAI Computing Conference
.