CC BY
0Contents lists available at ScienceDirect
Journal of the Mechanical Behavior of Biomedical Materials
journal homepage: www.elsevier.com/locate/jmbbm
Simulation of arterial dissection by a penetrating external body using cohesive zone modelling
Christopher Noblea'*, Olaf van der Sluisb, Ruud M.J. Vonckenb, Oliver Burkec, Steve E. Franklinb Roger Lewisc, Zeike A. Taylor3
a CISTIB, Centre for Computational Imaging and Simulation Technologies in Biomedicine, INSIGNEO Institute for in silico Medicine, Department of Mechanical Engineering, The University of Sheffield, Sheffield, UK b Philips Research, Eindhoven, The Netherlands
c Department of Mechanical Engineering, The University of Sheffield, Sheffield, UK
CrossMark
ARTICLE INFO
Keywords:
Catheter induced dissection Porcine aorta Cohesive zone model Finite element method Critical energy release rate Critical opening displacement
ABSTRACT
In this paper, we study the dissection of arterial layers by means of a stiff, planar, penetrating external body (a 'wedge'), and formulate a novel model of the process using cohesive zone formalism. The work is motivated by a need for better understanding of, and numerical tools for simulating catheter-induced dissection, which is a potentially catastrophic complication whose mechanisms remain little understood. As well as the large deformations and rupture of the tissue, models of such a process must accurately capture the interaction between the tissue and the external body driving the dissection. The latter feature, in particular, distinguishes catheter-induced dissection from, for example, straightforward peeling, which is relatively well-studied. As a step towards such models, we study a scenario involving a geometrically simpler penetrating object (the wedge), which affords more reliable comparison with experimental observations, but which retains the key feature of dissection driven by an external body, as described. Particular emphasis is placed on assessing the reliability of cohesive zone approaches in this context. A series of wedge-driven dissection experiments on porcine aorta were undertaken, from which tissue elastic and fracture parameters were estimated. Finite element models of the experimental configuration, with tissue considered to be a hyperelastic medium, and evolution of tissue rupture modelled with a consistent large-displacement cohesive formulation, were then constructed. Model-predicted and experimentally measured reaction forces on the wedge throughout the dissection process were compared and found to agree well. The performance of the cohesive formulation in modelling externally driven dissection is finally assessed, and the prospects for numerical models of catheter-induced dissection using such approaches is considered.
1. Introduction
Catheterisation is a common procedure in treatments for a variety of cardiovascular diseases. As it is a minimally invasive procedure, it puts less strain on the patient and is far less costly than other treatment methods (i.e. if viable, stenting a diseased coronary artery is preferable to a coronary artery bypass graft and opening of the chest). However, in such procedures the catheter can damage the vessel wall, ranging from rubbing away at the glycocalyx and endothelial layer, to vessel trauma such as catheter-induced dissections (CIDs) - see Fig. 1. Coronary artery dissection that extends to the aortic root occurred in 9 out of 43143 cardiac catheterisations (approximately 0.02%) (Dunning et al.,
2000). However, this increased to 0.2% if patients were undergoing treatment for myocardial infarction. Another study (GOMez-Moreno
et al., 2006) found the same type of dissection in 15 of 12031 percutaneous coronary intervention (PCI) procedures (approximately 0.12%). Data for this type of dissection are most prevalent because of its high mortality rate, however for coronary artery dissection alone, or CID of other arteries, data are sparser. In addition, it is difficult to provide context for these figures in terms of total complications, as CID are often included as part of other categories such as vascular trauma (Dellimore et al., 2014). However, in Stathopoulos et al. (2009) the overall complication rate for PCI was estimated at 2.7% for the period 2003-2006, from which it may be loosely inferred that coronary CID extending to the aortic root occurs in 5% of PCI procedures. Furthermore, the mortality rate is relatively high: 67% of deaths from coronary catheterisation procedures reported in Jain et al. (2002) were from CID of the coronary arteries, while (Fiddler et al., 2015) reported
* Corresponding author. E-mail address: mta08cn@sheffield.ac.uk (C. Noble).
http://dx.doi.org/10.1016/jjmbbm.2017.03.004
Received 13 December 2016; Received in revised form 23 February 2017; Accepted 5 March 2017 Available online 06 March 2017
1751-6161/ © 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/).
Fig. 1. Cross-sectional schematic of an artery undergoing CID: the catheter enters the false lumen created by the dissected tissue flap, and acts further to propagate the dissection.
that 25% of aortic CIDs that require surgical intervention lead to death of the patient. In this study, as a step toward development of numerical tools for simulating the complex catheter-tissue interactions involved, we investigate a potential approach to modelling such dissection processes. Such tools may in turn be valuable in virtual design processes for new devices, aimed at mitigating catheter-induced trauma.
Unfortunately it is difficult to ascertain how CIDs are initiated or propagated. This is because, commonly, catheterisation procedures are observed via three dimensional rotational angiography, and the dissection can only be detected by its effect on the blood flow. Thus, visualising the blood vessel wall and the damage to it is very difficult. Furthermore, the spatial resolution of such imaging modalities is relatively low, meaning it is difficult to identify the dissection until its effect on the blood flow is pronounced. As a result there is little information in the literature on mechanisms by which CIDs in arteries are initiated or propagated. Therefore, means by which a CID is propagated can only be speculated (Fig. 1). The cross-sectional view of the arterial wall and hollow, cylindrical catheter can be seen as a dissection is propagated. The catheter forces itself between the flap and vessel wall and acts to force the two apart. Reliable and well validated models are a possible approach to investigate these processes.
Arterial dissection by other mechanisms has been well investigated by both ex vivo testing in a laboratory, and simulations with the finite element (FE) method. The first ex vivo arterial wall dissection was performed via liquid infusion into the artery media (Carson and Roach, 1990; Roach and Song, 1994; Tiessen and Roach, 1993). This was aimed at replicating blood flow-driven propagation of the dissection. Later, controlled peeling of the media was employed to gain further understanding of the force displacement behaviour, and in particular its anisotropy (Sommer et al., 2008; Tong et al., 2011). This methodology was then utilised on diseased arteries and aneurysms to assess disease-associated changes (Tong et al., 2014; Pasta et al., 2012). Most recently the role of collagen and elastin fibres in resisting dissection was investigated by partially digesting the proteins and analysing the dissection behaviour using controlled peel testing (Noble et al., 2016).
FE-based computational models of these controlled peeling experiments have been proposed. In each case, a so-called cohesive zone (CZ) formulation has been used to model evolution of the tissue rupture, which enabled the effect of collagen fibre bridging during dissection to be emulated. Gasser et al. used a cohesive law within the extended FE method to simulate the controlled peeling (dissection) experiments of Sommer et al. Gasser and Holzapfel (2006), Sommer et al. (2008). Later, an anisotropic cohesive law was applied by Ferrara and Pandolfi that accounted for the anisotropic behaviour seen in Sommer et al.'s peel tests (Ferrara and Pandolfi, 2010). However, no work has been performed to investigate arterial perforation or dissection by an external body either experimentally or computationally.
In the present work, dissection of arterial layers driven by an external body, i.e. a catheter, is of interest. This process is fundamentally different from a peeling scenario as detailed in Fig. 2. Comparable configurations exist during, for example, needle insertion and cutting, and the CZ approach has been used by various authors to model these
processes. Needle insertion into muscle tissue has been investigated with parameters extracted from experimental data (Misra et al., 2008), however, no validation against experimental results was presented in turn. This was partially addressed by Forsell and Gasser when modelling cardiac tissue penetration and failure as a result of perforation by a pacemaker wire (Forsell and Gasser, 2011). However, it was not clear how the governing critical energy release rate (Gc) was calculated from the experimental data, making it difficult to assess the fidelity of their FE model with respect to the latter. Finally, Oldfield et al. utilised an experimental-computational approach to estimate Gc during needle insertion into a gelatine soft tissue phantom (Oldfield et al., 2013). This approach yielded good agreement between experimental and computational force-displacement profiles at four different insertion rates and was more comprehensive than the previous studies mentioned. However, gelatine is a homogeneous isotropic material that does not exhibit the stochastic breaking of fibres/fibrils seen in perforation and damage of most soft tissues (Tomlinson and Taylor, 2015), and so it is still unclear how well the CZ formalism carries over to real soft tissues in this context. These observations motivate the present study.
As mentioned, we are motivated by the need for reliable numerical tools for simulating CID of arterial tissues. This is a complex problem involving, most importantly, interaction between tissues and catheters, and tissue rupture. The usual complexities associated with soft tissue models, large deformations and nonlinear constitutive responses, obviously are also involved. Compared with the peel test experimental configurations, and corresponding models, that have so far been studied, the salient feature is the tissue-catheter interaction that drives the dissection. The complex geometry of catheter tips further complicates such models, and in particular renders experimental validation difficult. Specifically, the apparent unphysical penetration of 'intact' (i.e. not actually ruptured) CZ elements by the rigid contact body that renders the predicted contact and fracture conditions ambiguous. In Oldfield et al. (2013) they noted that reducing the critical opening of the CZ elements can mitigate this but at the expense of model stability. Thus, it would be beneficial to understand the effect of the critical opening displacement on overall the model behaviour when utilising CZ elements to model problems such as soft tissue perforation and CID. We therefore study a mechanically simpler problem, as a first step, that nonetheless retains the key feature of interaction between rupturing tissues and a relatively stiff external body: dissection of aortic tissues by means of a flat, rigid wedge with a rounded edge. The simpler geometry of the contacting bodies in this configuration means the fidelity of the model with respect to experimental conditions may be ensured with more confidence. The scenario may also be approximated by a plane strain formulation, further simplifying analysis. We studied the performance of the CZ approach to rupture modelling in this context.
A series of wedge-induced dissection experiments was first conducted (Section 2) to obtain data with which to validate subsequent models. A FE model of the experimental configuration, employing a CZ description of the rupture process was then developed (Section 3). Next, a novel parameter estimation procedure was formulated, based
Fig. 2. 2D schematic illustrating the differing behaviour at the tissue dissection front for controlled peeling (left) and wedge driven dissection (right). Left, the direction the tissue is pulled is the same as the direction of tissue rupture due to the simple loading conditions. Right, the direction of tissue rupture is approximately perpendicular to the direction the tissue is pulled as the tissue deforms around the wedge. This will likely result in differing fracture behaviour, as a result of the different fracture modes in each case.
on simplified FE models of the various experiment phases (Section 4). Finally, the CZ model predictions were compared with experimental observations (Section 5). The remainder of the paper is organised under these headings, and concludes with a discussion of the viability and limitations of CZ formulations in such scenarios.
2. Wedge dissection tests
2.1. Experimental configuration
Ex vivo dissection of samples of porcine aorta by means of a rigid wedge was conducted using the apparatus illustrated in Fig. 3. Partially dissected tissue samples were clamped between the two towers, as shown, and an aluminium wedge was driven between them, forcing their dissection. Wedge dimensions were: length: 95 mm length, width: 40 mm, thickness: 2.5 mm. The thickness was chosen to approximate the diameter of an endovascular catheter. The front edge of the wedge that propagates the dissection was rounded.
2.2. Specimen preparation
Tissue samples were tested no more than 24 h after the slaughter of the animal and were stored in a refrigerator, in saline solution, in the interim. The tissue was prepared by removing excess connective tissue and cutting axially such that the vessel wall opened into a rectangular section. This was cut into 40 mm by 20 mm specimens and washed three times in saline solution. Each sample was partially dissected by creating a small incision with a scalpel and then manually pulling apart the newly created tongues, allowing enough free tissue on either side of the tear to be held by the testing apparatus (Fig. 3). The un-dissected sample length (distance between tear opening and end of the sample) was approximately 10 mm. This was chosen to provide sufficient tissue to dissect, allowing a plateau region in the force-displacement curve to develop, while ensuring wedge displacement, before full separation of the tissue halves, was not so large that the wedge reached the base of the holding apparatus. All samples were partially dissected immediately prior to the test and were again placed in saline solution to moisten before being placed in the holding apparatus. Samples were kept out of saline solution for as short as possible and if necessary rehydrated by pipetting saline solution onto the sample before testing. To ensure the tear stayed as close to the centre of the media as possible, the initial dissection was carefully performed such that the tissue tongues were of equal thickness. Additionally, the tissue was continually inspected during the test to ensure the two tongues remained of similar thickness. Finally, by allowing only 10 mm of tissue to be dissected the likelihood of the tear straying from the media was reduced. The tissue was initially dissected and mounted such that the
wedge would progress in the vessel axial direction. Reaction forces on the wedge were measured using a 20 N load cell as the wedge was displaced at 1 mm/s, and until the tissue was completely dissected. If there was an error the run was stopped and the data were discarded. Examples of errors include the tissue slipping from the clamp or the tissue tearing irregularly. Results for 11 specimens were recorded in this way.
2.3. Results
The resulting force displacement plot is presented in Fig. 4. The shape and magnitudes of the curves are highly erratic due to the stochastic dissection behaviour which has been noted previously (Sommer et al., 2008; Noble et al., 2016).
However, whilst there are differences in magnitude and profile there are shared features. All curves initially (up to displacement UI) show the classic increasing stiffness characteristic of fibrous tissues (Fung, 1981), which results from unfolding and subsequent stretching of collagen and elastin fibres. This is followed by a levelling off, as fracture begins (displacement UI to UT), and then a sporadic plateau region as fracture continues (displacement UT to UP). Finally there is a drop, either sudden or gradual, depending on the extent of fibre bridging, as the wedge fully separates the tissue halves at UE. These features are identified in Fig. 4.
The profiles of the plateau regions were significantly noisier than those obtained from controlled peel testing (Noble et al., 2016), and showed greater variation in magnitude between samples. Additionally, in the peel tests, there was a sudden drop at full separation of the tissue tongues allowing easy identification of the end of the plateau region. In the wedge tests this drop is more gradual, likely because of individual fibres breaking less easily in this loading regime, and of the complex interaction between tissue and wedge.
Raw experimental curves were averaged, as shown in Fig. 4, and the latter values were used in subsequent analyses.
3. Computational modelling of wedge-induced dissection
The wedge-induced dissection process was simulated using a CZ-based FE model, created within MSC MARC 2015. One half of the sample was modelled, reflecting symmetry about the X-axis. A schematic of the model, associated dimensions, and boundary conditions is shown in Fig. 5. Wall thickness was chosen based on a representative literature value of 2 mm (Gundiah et al., 2013; Li et al., 2004). It has been previously shown for the peel test data that there is no correlation between wall thickness and dissection force (Sommer et al., 2008), and it can therefore be assumed that wall thickness plays an insignificant role in the dissection response
Fig. 3. Top. Apparatus for holding tissue in place. Clamps hold the tissue tongues, and the tissue to be dissected lies in the gap between the two towers. Sand paper was used to prevent slippage. Middle. Schematic detailing the tissue dimensions before and after the initial partial dissection and before mounting on the apparatus. Bottom. 2D schematic of wedge dissection test before loading is applied (left) and during the test (right). The set-up extends out from the page. L is the initial tissue length (10 mm). The loading conditions in Noble et al. (2016) can be replicated from the bottom left schematic by removing the wedge and pulling the clamps horizontally.
itself. However, to confirm this assumption, we separately investigated the sensitivity to tissue thickness of Gc values estimated using our procedure, and found that it was insignificant. Plane strain conditions were enforced, with out-of-plane thickness of 20 mm, corresponding to real specimen dimensions. The model consisted of 12030 4-node quadrilateral elements and 1920 4-node CZ elements. Fine and coarse regions of the mesh were connected by nodal ties as shown in Fig. 5. Mesh convergence analyses were conducted for all simulations. The model was assumed frictionless because no behaviours associated with high friction (e.g. stick slip and sudden jolts due to release of tissue held by frictional forces) were observed during the dissection experiments (Section 2.1),
when viewed under a low magnification stereo microscope. Additionally, all models assumed symmetry along the dissection path. However, while the dissection propagated in the centre of the tissue, arterial dissection is a stochastic process that will result in asymmetric tearing. This was partially rectified by using average data such that it is hypothesised that the mean dissection path will be in the centre of the tissue. Finally, despite the methods used to prevent it, minor slippage may have occurred at the grips in the experimental procedure. Any tissue that has slipped will be stretched and thus a small slippage can result in a notable increase in UP. In this work fixed boundary conditions were applied to mimic the clamped tissue, however, there may be slight discre-
* = 2 -^^ a
p=1 ap
+ x2p + x;p -3) + -k(J1
Fig. 4. Experimentally measured reaction force/width versus wedge displacement for the 11 specimens tested. Raw experimental curves (grey) are overlaid with averaged curve (black). Also detailed on the average plot are the four identified stages in the wedge propagated dissection process: initial stretching (0 to UI), before tissue begins to split; transition (UI to UT), where damage is initiating, identified by a decrease in the curve gradient; plateau region (UT to UP), wherein tissue is splitting and dissection is propagating; separation (UP to UE), characterised by a relatively sudden drop in reaction force, as the tissue tongues separate completely. The suddenness of the latter phase depends on the density and strength of the remaining bridging fibres between the tissue tongues.
pancy between model and experimental separation resulting from this.
3.1. Elastic tissue modelling
Tissue elastic response was modelled using an Ogden strain energy potential W (Msc Software Corporation, 2014; Ogden, 1972). Trial fits suggested that a two-term formulation was adequate to achieve good agreement with experimental results:
where fip are shear moduli and ap are dimensionless constants such that the initial shear modulus ^ = 2"£p=i Mpap. ^3 are the
isochoric (volume preserving) principal stretches, J is the Jacobian, and K is the bulk modulus. K was assigned a value of 109 kPa, representing near incompressibility.
More complex strain energy functions, capable of describing a wider range of material responses (including, for example, anisotropic effects) are available: two such models, with differing phenomenologi-cal backgrounds, are described in Gasser and Holzapfel (2006), Bischoff et al. (2002), for example. However, these provide no advantage in the present scenario, in which the deformation and corresponding material response are relatively simple. It should be noted, though, that the CZ and bulk material behaviours are separate from one another, and the CZ formulation described below can be combined with more complex bulk material models where the application requires them.
3.2. Cohesive zone formulation
CZ models are based on the idea of a process zone within the material, ahead of the physical crack tip. Material fracture is regarded as a gradual process in which formation of new crack surfaces takes place along an extended virtual crack, and the separation of these surfaces is resisted by cohesive tractions (De Borst et al., 2012). The traction across the CZ increases with increasing separation until a maximum, after which it irreversibly decreases until it reaches zero, and is governed by a traction-separation law. The aortic media consists of highly organised repeating lamellae held together by bridging collagen fibres that show a gradual decrease of cohesive force after exceeding a load limit (Gasser and Holzapfel, 2006). This behaviour aligns well with the CZ concept and consequently dissection has often been simulated using this method (Gasser and Holzapfel, 2006; Ferrara and Pandolfi, 2010, 2008; Gasser and Holzapfel, 2007).Large interface separations, accompanied by gradual failure of bridging fibres, are encountered during the dissection process. Therefore, the large deformation CZ formulation presented by Van den Bosch et al. was employed (Van Den Bosch et al., 2007). In this formulation, unlike small strain formulations, the opening and traction vectors are resolved globally rather than with respect to a local basis, thus no distinction is
Fig. 5. Top. Wedge dissection FE model with CZ elements immediately before force displacement data were recorded. The tissue strip was initially horizontal. The nodes indicated with a dotted line were first displaced into the position indicated, to represent the conditions at the start of the test. CZ elements are of zero thickness and their location is indicated by the dot-dash line. The wedge contact body displacement direction is indicated and the contact body extends out to the right of the image. Bottom. Initial model configuration, illustrating the fine and coarse meshes, which were bonded by means of nodal ties.
Fig. 6. Schematics of the wedge dissection test in the initial stretching phase and at the end of the plateau region, identified in Fig. 4 (for Lp = Lp = L). Symmetry boundary conditions were imposed along the dot-dash line, the wedge contact body displacement direction is indicated and the contact body extends out to the right of the image. A: the model immediately prior to contact initiating. As for the complete dissection model (Fig. 4), the initially straight tissue strip is first displaced into the configuration shown, to represent conditions at the start of the test. B: the model at the end of the simulation, after displacing the wedge by UI. C: the model immediately prior to contact initiating, and after displacing the indicated nodes (dotted line) into the starting test configuration. D: the model at the end of the simulation, after displacing the wedge by UP.
made between normal and tangential directions. The traction vector t = te and separation vector Ô = Se, where e is a unit vector along the line between associated opposing points of the interface, are related by
t =GIS W
S I 5r I
where Sc is a critical opening displacement, and Gc is the critical energy release rate. Evaluating (2) at the critical opening point Sc yields a corresponding critical traction tmax for the material:
¿>exp(1) (3)
Gc can be acquired from the plateau region of the experimental force
displacement curve (UT to UP). Sc (or tmax) can be found from either the initiation region (0 to Uj) or from local deformations (Van Der Sluis et al., 2011). The ramifications of the parameter selection in this study are discussed in Section 5.
4. Model parameter estimation
4.1. Elastic tissue parameters
The elastic parameters (/p, ap) of the tissue are identified from the experimental data by an inverse procedure. Fig. 4 illustrates the main stages identified during wedge-driven dissection. The MATLAB optimisation toolbox was utilised to fit the constitutive response to data in
t......=
Fig. 7. Reaction forces and estimated energies for Lp and LP, identified by the respective numbering. The solid black curve shows the average measured wedge reaction force/width versus displacement; the grey dashed curve shows the model-predicted response during the initial stretching phase, using optimised Ogden constitutive parameters; the hashed area corresponds to the elastic energy Welastic stored in the stretched tissue just prior to separation; finally, the grey shaded area corresponds to Wext. UI, UT, UP and UE are again identified for clarity.
Table 1
Dissection energy Gc, critical opening displacement Sc and fitted 2-term Ogden strain energy function parameters calculated for Lp+ and Lp (columns indicated by 1 and 2, respectively). Dissection energies are compared to the mean ± standard deviation from Noble et al. (2016).
Wedge driven dissection
Noble et al. (2016)
Gc (mJ/mm2) 0.738 0.803 0.183 ±0.064
8c (mm) 0.417 0.417 -
№1 (kPa) 0.12 6.19 -
a1 (mm) 8.60 2.60 -
№2 (kPa) 13.69 21.27 -
a2 (mm) 5.41 4.04 -
relevant regions of the overall force-displacement profile; specifically, throughout the initial elastic deformation stage (up to displacement Ujr), and at the estimated point of separation (UP). Two simplified FE models, corresponding to the initial and separation portions of the tests were constructed to enable estimation of corresponding reaction forces (Fig. 6). Including the latter point (UP) in the fit data helped to improve agreement of subsequent model and experimental results across the entire profile. However, it must be noted that both UI and UP are identified from the (average) experimental curve by manual inspection, and that there is naturally therefore some level of uncertainty in their values. This is particularly so in the case of UP, since the averaged experimental curve displays a more gradual decay to zero force, rather than a sharp drop (probably, a consequence in turn of irregular breaking of the last fibrous components near the end of the dissection). This ambiguity in the selection of UP, and correspondingly in the true length of the dissected tissue at this point, is accounted for in the analysis below.
4.1.1. Models for initial and separation stretch phases
In the initial stretching phase, the two clamped tissue tongues are stretched as a result of the downward displacement of the wedge. The corresponding FE model, at the beginning and end of the process, is depicted in schematics A and B, respectively, in Fig. 6. The model
consisted of 1230 4-node quadrilateral elements. The contact body is displaced in the -X direction to a value of UI. The wedge-tissue interface is assumed frictionless.
Just prior to separation of the two tissue tongues, the wedge load is completely supported by the stretching of those strips of tissue. Therefore, their summed reaction forces should equal the corresponding reaction force on the wedge measured in the experiments. Their un-stretched length is LP + (initial tongue length), where LP is the (undeformed) length of tissue dissected by the wedge at UP. However, as described, LP is uncertain, because of ambiguity in the location of UP in the experimental force-displacement curve. Here, therefore, we assume it lies within some bounding values, Lp and Lp, as follows. The upper bound Lp+ assumes that at the selected UP the tissue has indeed completely separated and thus Lpp = L . For the lower bound Lpp, we assume firstly that complete separation does not in reality occur till the end of the curve (UE), and secondly that the tissue rupture progresses linearly with displacement between UT and UE. At displacement UP, then, we would have:
L(UP - UT)
For the average data in Fig. 4, (4) yields Lp = 6.61 mm. The second FE model, simulating the experimental configuration at this point, is shown in Fig. 6 for Lp = Lp = L. The models comprised 1389 and 1515 4-node quadrilateral elements, respectively, for the cases of Lp = L- = 6.61 mm and LP = Lp = L mm. The contact body is displaced in the -X direction to a value of UP. Again, the contact was assumed frictionless.
4.1.2. Curve fitting
The constitutive parameters a1, a2) were then fit to the average experimental (F-U) curve. This process was formulated as the optimisation problem
x = argmin|w1 || Fexp - Fc
■ w2(Fe
where x is the vector of tissue model parameters and ||^||2 denotes 2-norm. F"p is a vector of experimental reaction forces on the wedge in the initial stretching region (up to UI) and F"""p is a corresponding
Fig. 8. Top. FE model midway through the analysis. The light grey region shows the continuum elements and the dark grey the CZ elements. Bottom. Von Mises stress contour plot of model at the same time point. Scale is in MPa and, for clarity, cohesive zone elements have not been visualised.
1st region 2nd\ 3rd region 4th region
regidn
Fig. 9. Left.FE model from Fig. 8 zoomed in at crack tip. The four stages of the crack progression are identified. In the first region there is a more traditional "sharp" crack opening. In the second region the crack configuration becomes more complex, as the tissue begins to deform around the contact body. Subsequently, the tissue wraps around the contact body (region 3) and finally separates from the contact body (region 4). Right. Corresponding view of tissue deformation around contact body in the experiment, for comparison. While such an image cannot show the detailed region 1 and 2 behaviour directly, the wrapping of tissue around the leading edge, and its subsequent separation from the wedge (regions 3 and 4) can be seen.
vector of model-computed reaction forces (Fig. 6). Fsepexp is the mean Table 1. experimental force just prior to separation (displacement UP) and
Fsepcomp is the corresponding force within the "separation" FE model (Fig. 6). Finally, w1 and w2 are weighting parameters, manually tuned to achieve best fit to the initial elastic region whilst ensuring the force at separation closely matches the experimental value (w1 = 3 and w2 = 5). The fitting was performed over multiple, randomly assigned, initial guesses and the best fit was selected. This increases the likelihood that the estimated parameter values indeed represent a globally optimal set (Ogden et al., 2004). The fittings were performed for Lp and L- separately, yielding two sets of constitutive parameters.
The fitted curves, overlaid on the experimental data, are shown in Fig. 7 and the optimised constitutive parameters are presented in
4.2. Critical energy release rate
To compute the critical energy release rate Gc associated with dissection of the tissue, we adapted the expression from Sommer et al. (2008), Tong et al. (2011), Noble et al. (2016), developed for analysis of peel test data, to give
Ar (Wext - WeImtic)
where Ar is the energy expended per unit width and Aa the crack
S 0-3 g
II f ' 1 1 i
/ 1 1 1 1 ' 1 // '
— Exper mental data del data (G*) del data (G¿)
//' 1 1 - CZ mc
0 5 10 15 20 25 30 35 40
Displacement (mm)
Fig. 10. Plots of the average experimental force-displacement data with responses predicted by the CZ model for G+ and G— in Table 1.
length, at an arbitrary point during the fracture process. Evaluating at UP in the experimental procedure gives Ar = Wext - WeIastic and Aa = LP. Wext is the external work per unit width done on the system,
f P: F da
is the work expended by stretching the tissue up to the point of complete separation of the tissue layers, per unit width (D is the model domain, P the first Piola-Kirchhoff stress tensor, F the deformation gradient tensor and w the sample width). Wext is equated with the area beneath the force displacement curve up to the displacement UP. Wejastic is more complicated to compute, given the nonlinearity of the tissue mechanical response and the non-homogeneity of the deformation. In this study, it was estimated using FE analysis.
Wejastic is assumed as the total strain energy (per unit width) in the "separation" FE models (described in Section 4.1.1). Two values, W+astic and W—astic, corresponding to the bounding values for LP, were computed. The corresponding values for Gc are then:
(Wext — Welastic) - _ r . n
as represented in Fig. 7. The resulting Gc values are given in Table 1 alongside those presented in Noble et al. (2016), which were derived from peel testing of the same tissue type, and with similar specimen geometries.
It can be seen that G— is slightly larger than G+ (which reflects the slightly larger difference between L— and L+ than between W—astic and
Both G+ and G— are approximately four times larger than the dissection energy values reported in Noble et al. (2016). One reason for this difference is that frictional energy dissipation was assumed negligible in the wedge tests, whereas its contribution to Wext cannot be zero in practice. However, a previous study (Takashima et al., 2007) has identified a friction coefficient of 0.046 between steel and artery wall, thus we believe this contribution alone cannot fully account for the greater Gc in the wedge configuration. Moreover, as mentioned, the tissue was well hydrated throughout the experiments and no obvious friction-associated behaviours were observed. Additionally in Noble
Fig. 11. Illustration of the effect of 8c on the force/width vs displacement behaviour of the CZ models for the case of G+. Top: force/width vs displacement for 5c = 2.5/6 mm, 5c = 2.5/8 mm and 5c = 2.5/2 mm alongside the average experimental data. Bottom: FE CZ model with 5c = 2.5/2 mm. The CZ elements (dark grey) hug around the contact body due to the remaining traction.
et al. (2016), Wejastic was calculated by assuming a linear constitutive response (based on the presentation in Sommer et al. (2008)), while here the material behaviour was captured using a more realistic hyperelastic constitutive model. Thus, the Weiastic contribution is unlikely to be the same for each method. We suggest, however, that the dominant reason for the discrepancy is the difference in the loading patterns experienced by the rupturing tissues in each configuration. This follows from the consideration that the (highly heterogeneous) tissue's overall fracture behaviour will be governed by the microstructural arrangement and tissue's overall fracture behaviour will be governed by the microstructural arrangement and behaviour of its constituents (collagen/elastin fibres, ground substance, etc.). The latter may include both inherent constituent fracture energies and dissipative processes active within the so-called process zone of the fracture. If these phenomena were modelled explicitly (Vossen et al., 2016, in the context of metal-elastomer separation), different specimen loading regimes (as between peel tests and the current wedge configuration) would naturally lead to different overall fracture behaviour, since the corresponding constituent-level loads would similarly differ. In the adopted CZ approach, we effectively ignore these microscale processes, considering the tissue macroscopically and with fracture behaviour governed by a single fracture energy parameter Gc. Gc, in other words, lumps the effects of the microscale processes. Correspondingly, different loading configurations can be expected to manifest as different values for Gc in such a formulation, as observed here.
5. Dissection model results and discussion
The deformed CZ-based model configuration partway through the analysis is shown in Fig. 8. The behaviour reflects the experimental observations wherein the tissue is stretched away from the wedge once it has separated from its opposing half. Additionally, Fig. 9 shows a
close-up of the crack tip illustrating the non-uniform complex fracture behaviour described in Section 4.2.
Comparison of numerical and experimental results. The experimental results are compared with numerical predictions for our calculated value of Gp and Gp (Table 1) in Fig. 10. Good agreement can be observed. This lends confidence to the proposed approach for finding Gc and in the possibility of extending the approach to other types of damage by an external body, such as needle penetration or cutting with a sharper object such as a scalpel. Additionally, the result was verified by calculating Gp, according to (8), from the associated CZ model result. Good agreement between predicted and experimental Wext, Wpastic, and Gp values is again observed, lending further confidence to the proposed approach.
Effect of Sc value. Care must be taken when selecting Sc for problems in which a contact body acts on the CZ elements. For this study a value of Sc = 2.5/6 mm was chosen (i.e. 1/6 of the 2.5 mm wedge thickness); thus, opening the elements to the wedge thickness itself would produce S = 6Sc, which is sufficient to damage them to the extent that t is effectively zero. Further reducing Sc (while maintaining Gc) has little effect on the force-displacement response, as demonstrated by the result of using Sc = 2.5/8 mm, shown in Fig. 11. However, if Sc is enlarged such that the contact body thickness is <6Sc, the predicted plateau force noticeably decreases, despite use of the same Gc value (Fig. 11). The reason is that in this scenario, the wedge contact body is able to penetrate CZ elements without opening them to their full damage level (i.e. without completely separating the two tissue halves). The resulting model configuration is also shown in Fig. 11. It can be seen that the continuum elements hug around the contact body as a result of the relatively high remaining traction within the CZ elements. The energy associated with this penetration is correspondingly less than would be expected for a given value of Gc, and the estimated reaction force during fracture is lower. Clearly, this penetration of "unruptured" elements is non-physical. A possible solution is to adopt the approach of Pagani and Perego (2015). Therein, directional CZ elements are deployed with an additional string element between the corresponding crack surfaces, representing bridging fibres. A node midway along this string element, at the contact point with the contact body, ensures the contact body cannot penetrate the CZ element unopposed.
Time dependent behaviour. Arterial tissue visco-elasticity has been well documented (Patel et al., 1973) and it is likely that different loading speeds will effect the gradient of elastic regions of curves (between 0 and UT in Fig. 4). In addition, rate dependent damage behaviour has been identified by peel tests in Tong et al. (2014) and this will likely influence the magnitude of the plateau region in the experimental data. The numerical results should therefore be understood to represent the behaviour at 1 mm/s - chosen as an estimate of speeds encountered during catheterisation, though specific data on the latter are understandably difficult to find. However, we have no reason to expect that such rate-sensitivity would qualitatively change the results, only the involved magnitudes. Finally, we note that loading speeds and environmental conditions (tissue hydration, for example) were very similar in the current experiments and in our earlier peel tests (Noble et al., 2016), and are therefore confident at least in excluding viscous effects as the reason for the differences in Gc values.
6. Conclusions
The aim of this study was to assess the validity of CZ elements for modelling rupture of soft tissues driven by stiff external bodies. This was motivated by the need to model processes of CID of the aortic wall, a previously little studied complication of catheter cardiovascular intervention.
We developed an experimental set-up in which an aluminium wedge of 2.5 mm thickness was displaced between an initially dissected aortic wall to dissect the remaining tissue, while the force displacement
response was recorded. The results of this experiment showed qualitative similarities to results from controlled peel testing, though with greater variability between samples and a less distinct plateau region during dissection propagation.
Upper and lower values for the critical energy release rate (Gp and Gcp) were then calculated using a novel experimental-computational approach utilising the FE method with the optimisation toolbox within MATLAB. Gc values so derived were found to be at least 3-4 times greater than those obtained from peel test experiments on similar test specimens (Noble et al., 2016). We hypothesise that this is mostly a consequence of differences in material resistance to fracture (as quantified by Gc, in this case) in the different fracture modes present in each test configuration.
Finally, Gcp and Gcp were used in a CZ FE model of the experimental procedure. It was found that the experimental data lay within the model-computed force displacement responses, giving confidence in the use of CZ elements for modelling cutting and perforation of soft tissues. However, care must be taken when selecting the critical opening/traction (Sc/tmax) to ensure the energy of the CZ elements is expended in a physically consistent way, or alternatively methods such as those described by Pagani et al. must be used (Pagani and Perego, 2015).
Acknowledgements
This work was supported by the Engineering and Physical Sciences Research Council (Doctoral Training Grant) and the European Commission Framework Programme 7, Understanding Interactions of Human Tissue with Medical Devices (UNITISS, FP7-PEOPLE-2011-IAPP/286174).
References
Bischoff, J.E., Arruda, E.A., Grosh, K., 2002. A microstructurally based orthotropic hyperelastic constitutive law. J. Appl. Mech. 69 (5), 570. http://dx.doi.org/10.1115/ 1.1485754.
Carson, W., Roach, M.R., 1990. The strength of the aortic media and its role in the propagation of aortic dissection. J. Biomech. 3 (6), 579-588. http://dx.doi.org/ 10.1016/0021-9290(90)90050-D.
De Borst, R., Crisfield, M.A., Remmers, J.J., Verhoosel, C.V., 2012. Non-Linear Finite Element Analysis Of Solids And Structures. John Wiley & Sons, Ltd.. http:// dx.doi.org/10.1002/9781118375938.
Dellimore, K.H., Franklin, S.E., Helyer, A.R., 2014. A review of catheter related
complications during minimally invasive transcatheter cardiovascular intervention with implications for catheter design. Cardiovasc. Eng. Technol. 5 (3), 217-232. http://dx.doi.org/10.1007/S13239-014-0183-9.
Dunning, D.W., Kahn, J.K., Hawkins, E.T., O'Neill, W.W., 2000. Iatrogenic coronary artery dissections extending into and involving the aortic root. Catheter. Cardiovasc. Interv. 51 (4), 387-393.
Ferrara, A., Pandolfi, A., 2008. Numerical modelling of fracture in human arteries. Comput. Methods Biomech. Biomed. Eng. 11 (5), 553-567. http://dx.doi.org/ 10.1080/10255840701771743.
Ferrara, A., Pandolfi, A., 2010. A numerical study of arterial media dissection processes. Int. J. Fract. 166 (1-2), 21-33. http://dx.doi.org/10.1007/S10704-010-9480-Y.
Fiddler, M., Avadhani, S.A., Marmur, J.D., 2015. Guide catheter-induced aortic
dissection complicated by pericardial effusion with pulsus paradoxus: a case report of successful medical management. Case Rep. Med. 2015, 1-6. http://dx.doi.org/ 10.1155/2015/480242.
Forsell, C., Gasser, T.C., 2011. Numerical simulation of the failure of ventricular tissue due to deep penetration: the impact of constitutive properties. J. Biomech. 44 (1), 45-51. http://dx.doi.org/10.1016/J.Jbiomech.2010.08.022.
Fung, Y.C., 1981. Biomechanics: mechanical Properties of Living Tissues. SpringerVerlag, New York.
GOMez-Moreno, S., SabatE, M., JimENez-Quevedo, P., VAZquez, P., Alfonso, F., Angiolillo, D.J., HernANdez-AntoliN, R., Moreno, R., BaNUelos, C., Escaned, J., Macaya, C., 2006. Iatrogenic dissection of the ascending aorta following heart catheterisation: incidence, management and outcome. Eurointervention 2 (2), 197-202.
Gasser, T.C., Holzapfel, G.A., 2006. Modeling the propagation of arterial dissection. Eur. J. Mech. - A/Solids 25 (4), 617-633. http://dx.doi.org/10.1016/ J.Euromechsol.2006.05.004.
Gasser, T.C., Holzapfel, G.A., 2007. Modeling plaque fissuring and dissection during balloon angioplasty intervention. Ann. Biomed. Eng. 35 (5), 711-723. http:// dx.doi.org/10.1007/S10439-007-9258-1.
Gundiah, N., Babu, A.R., Pruitt, L.A., 2013. Effects of elastase and collagenase on the
nonlinearity and anisotropy of porcine aorta. Physiol. Meas. 34 (12), 1657-1673. http://dx.doi.org/10.1088/0967-3334/34/12/1657.
Jain, D., Kurowski, V., Katus, H.A., Richardt, G., 2002. Catheter-induced dissection of the left main coronary artery. Nemesis Invasive Cardiol. Case Report. Rev. Lit. 91 (10), 840-845. http://dx.doi.org/10.1007/S00392-002-0823-1, [Zeitschrift FÜR Kardiologie].
Li, A.E., Kamel, I., Rando, F., Anderson, M., Kumbasar, B., Lima, J.A.C., Bluemke, D.A., 2004. Using MRI to assess aortic wall thickness in the multiethnic study of atherosclerosis:distribution by race, sex, and age. Am. J. Roentgenol. 182 (3), 593-597. http://dx.doi.org/10.2214/ajr.182.3.1820593, [URL (http://www.ncbi. nlm.nih.gov/pubmed/14975953>].
Misra, S., Reed, K.B., Douglas, A.S., Ramesh, K.T., Okamura, A.M., 2008. Needle-Tissue Interaction Forces For Bevel-Tip Steerable Needles., Proceedings of the 2Nd Biennial Ieee/Ras-Embs International Conference on Biomedical Robotics and Biomechatronics, pp. 224-231 http://dx.doi.org/10.1109/Biorob.2008.4762872.
MSC Software Corporation. MARC 2014 Volume A: Theory and User Information 2014.
Noble, C., Smulders, N., Lewis, R., CarrE, M.J., Franklin, S.E., Macneil, S., Taylor, Z.A., 2016. Controlled peel testing of a model tissue for diseased aorta. J. Biomech. 49 (15), 3667-3675. http://dx.doi.org/10.1016ZJ.Jbiomech.2016.09.040.
Ogden, R.W., Saccomandi, G., Sgura, I., 2004. Fitting hyperelastic models to
experimental data. Comput. Mech. 34 (6), 484-502. http://dx.doi.org/10.1007/ S00466-004-0593-Y.
Ogden, R.W., 1972. Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids. Proc. R. Soc. Lond. A: Math. Phys. Eng. Sci. 328 (1575), 567-583. http://dx.doi.org/10.1098/Rspa.1972.0096.
Oldfield, M., Dini, D., Jaiswal, T., Rodriguez, F., Baena, Y., 2013. The significance of rate dependency in blade insertions into a gelatin soft tissue phantom. Tribol. Int. 63, 226-234. http://dx.doi.org/10.1016/J.Triboint.2012.08.021.
Oldfield, M., Dini, D., Giordano, G., Rodriguez, F., Baena, Y., 2013. Detailed finite element modelling of deep needle insertions into a soft tissue phantom using a cohesive approach. Comput. Methods Biomech. Biomed. Eng. 16 (5), 530-543. http://dx.doi.org/10.1080/10255842.2011.628448.
Pagani, M., Perego, U., 2015. Explicit dynamics simulation of blade cutting of thin elastoplastic shells using directional cohesive elements in solid-shell finite element models. Comput. Methods Appl. Mech. Eng. 285, 515-541. http://dx.doi.org/ 10.1016/J.Cma.2014.11.027.
Pasta, S., Phillippi, J.A., Gleason, T.G., Vorp, D.A., 2012. Effect of aneurysm on the mechanical dissection properties of the human ascending thoracic aorta. J. Thorac. Cardiovasc. Surg. 143 (2), 460-467. http://dx.doi.org/10.1016/ J.Jtcvs.2011.07.058.
Patel, D.J., Janicki, J.S., Vaishnav, R.N., Young, J.T., 1973. Dynamic anisotropic
viscoelastic properties of the aorta in living dogs. Circ. Res. 32 (1), 93-107. http:// dx.doi.org/10.1161/01.Res.32.1.93.
Roach, M.R., Song, S.H., 1994. Variations in strength of the porcine aorta as a function of location. Clin. Investig. Med. 17 (4), 308-318.
Sommer, G., Gasser, T.C., Regitnig, P., Auer, M., Holzapfel, G.A., 2008. Dissection
properties of the human aortic media: an experimental study. J. Biomech. Eng. 130 (2), 021007. http://dx.doi.org/10.1115/L2898733.
Stathopoulos, I., Jimenez, M., Panagopoulos, G., Kwak, E.J., Losquadro, M., Cohen, H., Iyer, S., Ruiz, C., Roubin, G., Garratt, K., 2009. The decline in PCI complication rate: 2003-2006 versus 1999-2002. Hell. J. Cardiol. 50 (5), 379-387, [URL (http://www.ncbi.nlm.nih.gov/pubmed/19767279>].
Takashima, K., Shimomura, R., Kitou, T., Terada, H., Yoshinaka, K., Ikeuchi, K., 2007. Contact and friction between catheter and blood vessel. Tribol. Int. 40 (2), 319-328. http://dx.doi.org/10.1016/J.Triboint.2005.10.010.
Tiessen, I.M., Roach, M.R., 1993. Factors in the initiation and propagation of aortic dissections in human autopsy aortas. J. Biomech. Eng. 115 (1), 123-125.
Tomlinson, R.A., Taylor, Z.A., 2015. Photoelastic materials and methods for tissue biomechanics applications. Opt. Eng. 54 (8), 81208. http://dx.doi.org/10.1117/ 1.Oe.54.8.081208.
Tong, J., Sommer, G., Regitnig, P., Holzapfel, G.A., 2011. Dissection properties and mechanical strength of tissue components in human carotid bifurcations. Ann. Biomed. Eng. 39 (6), 1703-1719. http://dx.doi.org/10.1007/S10439-011-0264-Y.
Tong, J., Cohnert, T., Regitnig, P., Kohlbacher, J., Birner-Gruenberger, R., Schriefl, A.J., Sommer, G., Holzapfel, G.A., 2014. Variations of dissection properties and mass fractions with thrombus age in human abdominal aortic aneurysms. J. Biomech. 47 (1), 14-23. http://dx.doi.org/10.1016AJ.Jbiomech.2013.10.027.
Van Den Bosch, M.J., Schreurs, P.J.G., Geers, M.G.D., 2007. On the development of a 3D cohesive zone element in the presence of large deformations. Comput. Mech. 42 (2), 171-180. http://dx.doi.org/10.1007/S00466-007-0184-8.
Van Der Sluis, O., Hsu, Y.Y., Timmermans, P.H.M., Gonzalez, M., Hoefnagels, J.P.M., 2011. Stretching-induced interconnect delamination in stretchable electronic circuits. J. Phys. D: Appl. Phys. 44 (3), 034008. http://dx.doi.org/10.1088/0022-3727/44/3/034008.
Vossen, B.G., Van Der Sluis, O., Schreurs, P.J.G., Geers, M.G.D., 2016. High toughness fibrillating metal-elastomer interfaces: on the role of discrete fibrils within the fracture process zone. Eng. Fract. Mech. 164, 93-105. http://dx.doi.org/10.1016/ J.Engfracmech.2016.05.019, [URL http://dx.doi.org/10.1016/j.engfracmech.2016. 05.019].
