Keywords

1 Introduction

The coronary artery disease results from a physiological ageing process causing a progressive narrowing (stenosis) of the artery lumen. In this study, we focus on the treatment of the disease by balloon inflation (i.e. angioplasty) and stent deployment. The procedure outcomes highly depend on the stent final geometry and on the accuracy of the deployment. Consequently, numerical simulation, if reliable, can be a valuable asset to minimise complications.

We find in literature an important number of studies using the Finite Element Method (FEM) to simulate the deployment of coronary stents. Several studies aim at the stent design optimisation: see for instance [8] (stent alone), [3] (stent and balloon) or [12] (stent, balloon, and artery). As design optimisation gives priority to precision over computational time (e.g. 48 h in [3]), the corresponding models can not be used directly in a clinical environment.

In this paper, we propose a faster simulation routine to model accurately the deployment of a balloon-expandable stent. Our objective is to reach an execution time compatible with clinical routine, so that our simulation method could be used in practice as a relevant prediction tool (typically in patient-specific applications). We decide to base our model on the discretisation of the stent geometry by 1D serially linked beam elements. Similar work was proposed by Čanić and Tambača in [11], where the stent structure is modelled by an assembly of 1D rod elements. Although limited to a linear elastic behaviour (auto-expandable stents), the simulation is much faster than for design optimisation.

We decide to use beam elements in a similar way, in order to model the more complex (nonlinear) deformation undergone by a balloon-expandable stent. The use of connected beam elements has already been proved efficient in medical simulation to model slender structures, such as endovascular coils for brain aneurysm in [4] or flexible needles in [1]. The use of beam elements in these studies allows to achieve low computational times (up to interactive simulations), but once again is limited to elastic deformations.

2 Method

2.1 Stent Discretisation

In order to benefit from the slender shape of the stent struts, we use serially linked 1D beam elements to discretise the entire stent structure.

Each beam element is represented by two nodes, each of them described by 6 Degrees of Freedom (DoF), 3 for position, and 3 for orientation. The approximation of the continuous medium between the nodes can be written in matrix form as:

$$\begin{aligned} \mathbf {u}(x, y, z) = \mathbf {N}(x, y, z)\mathbf {\tilde{u}}, \end{aligned}$$
(1)

where \(\mathbf {u}\) is the continuous displacement field inside the element, (x, y, z) are the material coordinates, and \(\mathbf {\tilde{u}}\) is the \(12\times 1\) vector of nodal displacement. \(\mathbf {N}\) is a \(3\times 12\) matrix containing Timoshenko interpolation shape functions, which expression is given in [2]. The Timoshenko beam model notably allows to take into account shear deformations and describe ‘thick’ beams more accurately. We encounter this type of beams in highly curved parts of the stent.

The mesh conception is made with Computer Aided Design. The dimensions of the stent are retrieved from the micro-CT acquisition of a crimped metallic stent. From these measures, we can reproduce a flat version of the stent mesh, composed of 2D edges. The mesh is then wrapped up into a cylindrical shape, each edge corresponding to a beam element.

2.2 Finite Element Workflow

Our Finite Element (FE) implementation entirely relies on the open simulation framework SOFAFootnote 1. A thorough description of most of the mechanisms involved in the SOFA simulation workflow can be found in [5].

Very briefly, we base our simulation on the (implicit) Backward Euler numerical scheme. This leads to solving at each time step the following mechanical system:

$$\begin{aligned} \underbrace{(\mathbf {M}- h\mathbf {B}-h^2 \mathbf {K})}_{\mathbf {A}} \varDelta \mathbf {v}_{t+h} = \underbrace{h\mathbf {f}_t(\mathbf {x}, \mathbf {v}) + h^2 \mathbf {K}\mathbf {v}_t}_{\mathbf {b}}, \end{aligned}$$
(2)

in which \(\mathbf {M}\) is the mass matrix of the discretised mechanical system (i.e. if the system is described by n DOFs, then \(\mathbf {M}\) is a \(n \times n\) matrix), \(\mathbf {x}\) is the global position vector, and \(\mathbf {v}\) the velocity. Subscript \(_t\) indicates that the variable is considered at time t, and h is the time step. \(\mathbf {f}\) stands for the forces, \(\mathbf {K}= \frac{\partial \mathbf {f}}{\partial \mathbf {x}}\) for the stiffness matrix and \(\mathbf {B}=\frac{\partial \mathbf {f}}{\partial \mathbf {v}}\) the damping matrix.

The system in (2) is solved in \(\varDelta \mathbf {v}_{t+h}\) at each time step. From this, the simulation can be progressed forward in time.

2.3 Mechanical Model

Elasticity. In elasticity, the relation between the stress tensor \(\varvec{\sigma }\) and the strain tensor \(\varvec{\epsilon }\) is given by Hooke’s law, and can be written in matrix form as:

$$\begin{aligned} \varvec{\sigma }= \mathbf {C}\varvec{\epsilon }. \end{aligned}$$
(3)

The explicit expression of \(\mathbf {C}\) as a matrix can be found in [9]. We stress out that \(\varvec{\epsilon }\) and \(\varvec{\sigma }\) are actually second order tensors, which can be written as \(6\times 1\) vectors without loss of information thanks to symmetry.

From (3) and the virtual work principal, the \(12\times 12\) element stiffness matrix \(\mathbf {K_{e}}\) can be expressed as:

$$\begin{aligned} \mathbf {K_{e}}= \int _{\varOmega } \mathbf {B_{e}}^{T} \mathbf {C}\mathbf {B_{e}}~ \text {d}\varOmega , \end{aligned}$$
(4)

where \(\varOmega \) is the 3D domain on which the element is defined, and \(\mathbf {B_{e}}\) is a \(6\times 12\) matrix obtained by spatial derivation of \(\mathbf {N}\). Expression of \(\mathbf {B_{e}}\) results from the small strain hypothesis, which is expressed in [9]. For details of the computation of \(\mathbf {K_{e}}\) from the principle of virtual work, we also refer the reader to [9].

Plasticity. In plasticity theory, the constitutive law (3) varies during the deformation, following two phases:

  • an elastic phase, during which the internal stress depends linearly on the strain, as in (3),

  • a plastic phase, during which part of the strain can be distinguished as resulting from plastic energy dissipation, and the stress-strain relation becomes nonlinear.

The transition from the elastic to the plastic phase is described by a yield criterion \(f: \mathbb {R}^6 \rightarrow \mathbb {R}\), defined on the stress space. In this study, we use the Von Mises yield function, which expression is given in [6].

The computation of an acceptable plastic stress state relies on the decomposition of the strain into a plastic and an elastic component: \(\varvec{\epsilon }= \varvec{\epsilon }^{el} + \varvec{\epsilon }^{pl}\). In the following, we adopt a model of perfect plasticity, as described in [6]. The term ’perfect’ indicates that the yield surface, described by \(f(\varvec{\sigma }) = 0\) in the stress space, also remains unchanged.

The evolution of the plastic strain is described by an associative flow rule, expressed as:

$$\begin{aligned} \text {d}\varvec{\epsilon }^{pl} = \text {d}\lambda \frac{\partial f}{\partial \varvec{\sigma }}, \end{aligned}$$
(5)

where \(\text {d}\lambda \) is called the plastic multiplier and defines the norm of \(\text {d}\varvec{\epsilon }^{pl}\). We use the radial return algorithm as described in [6], to compute a new plastic stress state \(\varvec{\sigma }_{t+h}\) at each time step. Details of the computation can be found in [6].

Once \(\varvec{\sigma }_{t+h}\) is known, we may finally compute the resulting internal forces, and linearised stiffness matrix, so that they are taken into account in the global mechanical system (2).

The new internal forces are simply computed by integration:

$$\begin{aligned} \mathbf {f}_{int}(\varvec{\sigma }_{t+h}) = \int _{\varOmega } \mathbf {B_{e}}^{T} \varvec{\sigma }_{t+h} ~ \text {d}\varOmega . \end{aligned}$$
(6)

The linearised version of the stiffness matrix (or tangent stiffness matrix) \(\mathbf {K}_t\) can be expressed in a similar analysis way as for (4), leading to:

$$\begin{aligned} \mathbf {K}_t = \int _{\varOmega } \mathbf {B_{e}}^{T} \left( \frac{\text {d}\varvec{\sigma }}{\text {d}\varvec{\epsilon }} \right) \mathbf {B_{e}}~ \text {d}\varOmega . \end{aligned}$$
(7)

Starting from Hooke’s law, we have:

$$\begin{aligned} \text {d}\varvec{\sigma }= \mathbf {C}\text {d}\varvec{\epsilon }^{el} = \mathbf {C}(\text {d}\varvec{\epsilon }- \text {d}\varvec{\epsilon }^{pl}) = \mathbf {C}(\text {d}\varvec{\epsilon }- \text {d}\lambda \frac{\partial g}{\partial \varvec{\sigma }}). \end{aligned}$$

Using the consistency condition (\(\text {d}f = 0\)), we can express \(\text {d}\lambda \) in terms of \(\frac{\partial f}{\partial \varvec{\sigma }}\) and \(\text {d}\varvec{\epsilon }\), and replace it in the expression of \(\text {d}\varvec{\sigma }\), to obtain:

$$\begin{aligned} \text {d}\varvec{\sigma }= \left( \mathbf {C}- \frac{\mathbf {C}\frac{\partial f}{\partial \varvec{\sigma }} \frac{\partial f}{\partial \varvec{\sigma }}^{T} \mathbf {C}}{\frac{\partial f}{\partial \varvec{\sigma }}^{T} \mathbf {C}\frac{\partial f}{\partial \varvec{\sigma }}} \right) \text {d}\varvec{\epsilon }= \mathbf {C}^{ep} \text {d}\varvec{\epsilon }. \end{aligned}$$
(8)

This allows to explicitly express the tangent stiffness as:

$$\begin{aligned} \mathbf {K}_t = \int _{\varOmega } \mathbf {B_{e}}^{T} \mathbf {C}^{ep} \mathbf {B_{e}}~ \text {d}\varOmega . \end{aligned}$$
(9)

In practice, we use Gaussian reduced integration to compute the stress and the element internal forces with (6), at each time step during the simulation. We update the tangent stiffness matrix with (9). Once \(\mathbf {K_{e}}\) is computed for each element, all the matrices are assembled into the global stiffness matrix \(\mathbf {K}\) in (2).

2.4 Balloon and Contact Constraints

We simulate the balloon indirectly by attaching a spring to each of the stent nodes. The other end of each spring is attached to a particle on a virtual surface, representing the balloon membrane. At time \(t=0\), the springs rest lengths are defined to maintain the system in equilibrium.

During the simulation, the surface particles are moved in the radial direction so that the diameter of the virtual membrane corresponds to the one given in the balloon compliance table. This table, provided by the stent manufacturer, gives the expected stent and balloon diameters for given pressures (typically every 1 atm from \(P=8\) to \(P=20\) atm). During the simulation, we gradually increment a virtual pressure to reproduce the actual increase in pressure carried out by the physician. The characteristics of the springs (stiffness and damping) are fixed so that the stent diameter in the simulation is coherent with the compliance table. As the table does not provide information for low pressures (<8 atm), we complete it with two entries, corresponding to 0 atm, and to the first pressure measured experimentally for which the stent becomes fully cylindrical (\(p_{cyl}\)). We then use a linear-by-parts model for pressure increasing in the simulation: the first part is a linear interpolation between 0 atm and \(p_{cyl}\), while the second part is fitted on the \(p_{cyl}\) value and the constructor compliance table.

We handle contacts between the topological primitives composing the stent (edges) and the artery (triangles) with a collision pipeline available in SOFA. At each time step, if a collision between two primitives is detected, a response is computed accordingly in the form of a Lagrangian constraint.

Briefly, the addition of constraints involves an additional step in the numerical solution:

  • At first, (2) is solved without constraints, giving a new position (free motion).

  • Then, a modified constrained system is solved from the free motion:

$$\begin{aligned} \mathbf {A}\varDelta \mathbf {v}_{t+h} = \mathbf {b}+ h\mathbf {H}^{T}\varvec{\lambda }, \end{aligned}$$
(10)

where \(\mathbf {H}\) is the Jacobian of the constraint expressions, and \(\varvec{\lambda }\) contains the Lagrange multipliers associated to each constraint. In SOFA, this system is typically solved iteratively using a Gauss-Seidel algorithm. A more detailed description of the process is available in [5].

2.5 Experimental Validation

To propose a first assessment of the simulation, we conducted an experiment involving the apposition of a stent in a perfectly straight coronary artery phantom. We chose silicone (Sylgard 184, Dow Corning, Midland, Michigan, U.S.A.) to reproduce the artery, as the material has already been successfully used to create phantom vessels. In order to exactly control the vessel geometry, we used additive manufacturing to create a plastic mould and a water-soluble stick made in PolyVinyl Alcohol (PVA). The mould gives the external shape of the phantom, while the lumen geometry is entirely define by the stick (once dissolved).

Following the phantom creation, we actually expand a stent inside the mock artery, under micro-CT supervision. We proceed by gradually incrementing the inflation pressure of the balloon, and realising a CT acquisition at each incremental step. The resulting 3D reconstructions contain the geometry of the stent in transient states, from a pressure of 0 atm (crimped stent) to the nominal pressure (fully expanded stent). The deformation can be compared to the output of the simulation.

As the artery phantom was obtained from 3D-printed elements, the 3D artery geometry to be included in the simulation is known. In this case, we use a surface mesh, on which we add an elastic behaviour [10], already implemented in SOFA.

3 Results

As a first validation, we simulated the expansion of a Synergy II coronary stent (Boston Scientific, Marlborough, Massachusetts, U.S.A) in a perfectly cylindrical artery model, which we reproduced experimentally. For Hooke’s Law, we use \(E =\) 203 GPa (Young’s modulus) and \(\nu =\) 0.285 (Poisson ratio), corresponding to the Pt-Cr alloy [7] used for the Synergy II. The artery diameter is chosen to match the stent nominal diameter (3 mm). Pressure in the experiment is increased from 0 atm to 11 atm, while carrying out micro-CT acquisition at 5, 6, 7, 9 and 11 atm. We increment the pressure in the simulation by following the linear-by-parts model mentioned above. We export the stent geometry at regular pressure increments (every 200 iterations, i.e. every 0.44 atm) between 0 and 11 atm.

Comparison is made between a pointset extracted from the micro-CT 3D reconstruction, and the stent nodes. We register the two pointsets in two steps: first by superimposing the principal axes of inertia and centres of gravity (rigid transform), and then by minimising the Euclidean distance between simulated points and their closest neighbours in the CT pointset (rigid rotation around the common axis).

We compare the poinset of each acquisition to all the exported stent geometries. We choose as best correspondence the one minimising the mean radial distance \(\bar{d_r}\) and Euclidean distance \(\bar{d_E}\) between the pairs of points. Table 1 gives the best matching simulated stent state for each micro-CT acquisition, with the number of iterations (it), the corresponding simulated pressure and diameter, and the minimised metrics values. 3D rendering of the registration output is given in Fig. 1 for a crimped and a deployed configurations.

Table 1. Comparison between experimental and simulated stent expansions.
Full size table
Fig. 1.
figure 1

3D rendering of the rigid registration between the last micro-CT acquisition (deflated balloon) and the simulation best correspondence. The CT reconstruction is displayed in grey, while the simulation mesh nodes are superimposed as red points (Color figure online).

Full size image

4 Discussion and Conclusion

The first remark which has to be made regarding our results is the significant difference between the experimental pressure and the simulated pressure in the best corresponding configuration. This discrepancy can be explained both by the fact that we export geometries only every 200 iterations (thus limiting the number of candidates for best correspondence), and by the complete lack of data in the compliance table, for \(P < 8\) atm. Consequently, the first part (\(P \le \) 5 atm) of the model that simulates the pressure increase is too approximative. This is evidenced by the almost identical diameters measured in the acquisitions for \(P =5\) atm and \(P =6\) atm, which can’t correspond to a linear model.

An issue with developing a more accurate pressure/diameter model for low pressure values is that the stent experience asymmetrical deformation for \(P \le \) 5 atm. As presented by De Beule et al. in [3], this phenomenon, called dog boning can only be reproduced if a realistically folded balloon is simulated.

If we don’t consider early transient states, Zahedmanesh et al. showed in [12] that an approximated balloon model is able to accurately retrieve the final configuration of the stent. In this study, it is indeed the case for the nominal pressure of 11 atm, for which the acquisition and the simulation are in accordance. In addition to that, considering only the diameter independently of the pressure, we can see that the simulation is also able to reproduce transient states as long as the stent is fully cylindrical. On the correspondences we obtain, the distance errors are thus in the same order of magnitude as the stent strut thickness (50 \(\upmu \)m for the Synergy II). We consider this error reasonable regarding stent apposition.

Regarding the computational time, the plastic mechanical model we implemented can be handled at 4 Frames Per Second (FPS), allowing a stent mesh composed of 1622 beam elements to reach the artery diameter in about 15 min, without any form of optimisation (single thread CPU execution on an Intel Xeon E3-1270 v5 (8CPUs) 3.6 GHz, 16 GB RAM, Windows10 64 bits). As soon as contact occurs, the number of FPS decreases, but the computational time stays significantly lower than for fully volumetric meshes.

In this paper, we present a complete methodology to realistically simulate the plastic deformation of a balloon-expandable stent with 1D serially linked beam elements, in a reasonable time. The long-term goal of this work is to be included in a workflow where we could retrieve a patient’s artery geometry and mechanical properties using endovascular Optical Coherence Tomography (OCT) imaging, run the simulation, and provide the physicians with a personalised prediction of stent apposition, before deployment. For this, we aim at a computational time lower than a few minutes, coherent with the first results we obtained.

Presently, clinical assessment of stent apposition almost exclusively relies on the evaluation by the physicians of the artery diameter, before and after the procedure. In this context, we believe that a stent apposition simulation could provide additional (mechanical) information to the physician. In order to progress towards this clinical use, we started to develop simulations and experiments with more realistic artery geometries. In parallel, the use of OCT to complete the clinical workflow is currently under investigation. Finally, we also intend to improve the present simulation by optimising our code (addressing the issue of computational time), and compare the results with equivalent 3D codes from the literature.