## Abstract

Blood vessels have unique properties that allow them to function together within a complex, self-regulating network. The contractile capacity of the wall combined with complex mechanical properties of the extracellular matrix enables vessels to adapt to changes in haemodynamic loading. Homogenized phenomenological and multi-constituent, structurally motivated continuum models have successfully captured these mechanical properties, but truly describing intricate microstructural details of the arterial wall may require a discrete framework. Such an approach would facilitate modelling interactions between or the separation of layers of the wall and would offer the advantage of seamless integration with discrete models of complex blood flow. We present a discrete particle model of a multi-constituent, nonlinearly elastic, anisotropic arterial wall, which we develop using the dissipative particle dynamics method. Mimicking basic features of the microstructure of the arterial wall, the model comprises an elastin matrix having isotropic nonlinear elastic properties plus anisotropic fibre reinforcement that represents the stiffer collagen fibres of the wall. These collagen fibres are distributed evenly and are oriented in four directions, symmetric to the vessel axis. Experimental results from biaxial mechanical tests of an artery are used for model validation, and a delamination test is simulated to demonstrate the new capabilities of the model.

## 1. Introduction

The arterial tree is a highly complex, self-regulating network. Active cellular processes control extracellular matrix remodelling and smooth muscle contraction and thereby allow the vessels to respond to changes in blood flow, blood pressure, circulating hormones and neural activity. Indeed, cellular control of the matrix is fundamental to ensuring both appropriate mechanical functionality and structural integrity, and there is an intricate communication between the cells and matrix. Of particular note, the microstructure endows the normal arterial wall with a unique stability against lateral bending (tortuosity) while maintaining constant axial force during changes in internal pressure, which is energetically favourable [1].

The wall of a normal elastic artery consists of three layers: the intima, media and adventitia. The intima is composed of a layer of endothelial cells that line the inside of the vessel and thus contact the flowing blood. The media is the middle layer; it consists of multiple layers of smooth muscle cells and elastin with admixed glycosaminoglycans and collagen fibres [2]. The adventitia is the outer layer; it consists primarily of collagen fibres. The elastin is highly compliant and resilient and carries most of the load at low pressures, while the stiffer collagen fibres carry load at high pressures and provide most of the strength of the wall. Collagen fibre dispersion differs between layers and specific artery types. In general, however, these fibres, along with the predominantly circumferential smooth muscle orientation, give arteries a macroscopic mechanical response that is nearly orthotropic with respect to their primary axes, circumferential and axial [3]. Whereas multi-layered models of the wall are particularly important in mechanobiological studies, single-layered homogenized models are often sufficient to capture the effects of the wall on the haemodynamics.

Continuum models of arteries have advanced significantly over the past three decades and are now able to describe well the overall anisotropic mechanical properties. One of the earliest and most inspirational models was the Fung model [4,5], which introduced an exponential strain energy function that was incorporated into many subsequent models. Another often used model [6,7] employs a multi-constituent framework, namely a neo-Hookean component representing an elastin-dominated isotropic matrix plus families of parallel fibres that represent fibrillar collagen. The fibres are often described by Fung-type strain energy functions. Two identical fibre families oriented symmetrically to the axis of the cylinder are often used to render the material macroscopically orthotropic. One variant accounts for dispersion in fibre orientation by using a parameter to characterize the distribution of fibre orientation angles. A different approach [8,9] introduces two additional fibre families, aligned axially and circumferentially to create a ‘four-fibre family’ model. Without considering dispersion explicitly, the mixture of the four fibre orientations provides a good phenomenological approximation of the macroscopic effects of collagen fibre dispersion and cross-links, which are not considered explicitly in current models [9]. The circumferential fibres also account for the passive elastic properties of the smooth muscle cells, which are typically aligned in the circumferential direction.

Notwithstanding the success of these and other recent models, it is worthwhile to consider additional ways to advance the characterization of arteries, leveraging what has been learned through prior continuum models. To this end, we propose a discrete particle model for an anisotropic artery based on the same fibre-reinforced elastic structure. While the continuum approach offers certain advantages, it presents challenges that can be avoided or at least lessened in discrete modelling. For instance, the problem of dissection, a pathological condition in which layers within the arterial wall separate, is difficult to simulate with continuum models but trivial to represent in a discrete model. Among discrete models, the peridynamic theory was proposed to overcome some of the intrinsic limitations of classical continuum theories when dealing with problems having discontinuous displacement fields, such as fracture [10,11]. Similarly, smoothed particle hydrodynamics (SPH) was introduced for solid mechanics as a high-strain Lagrangian method for modelling fractures and crack propagation [12].

While peridynamics and SPH formulations are motivated by methods of continuum mechanics (top-down approach), dissipative particle dynamics (DPD) [13,14] is a discrete particle mesoscale (bottom-up) approach in which particles represent clusters of molecules and finite-strain solid deformations are described using nonlinear spring potentials between neighbouring particles. DPD is an effective method for modelling the mechanics of soft matter (e.g. polymers) and complex fluids like blood, which suggests a potential advantage when developing fluid–solid interaction models for the vasculature. For example, existing DPD models of red blood cells and platelets [15–17] that accurately capture contributions of haematocrit to the viscosity of blood or contributions of platelet activation and aggregation to thrombosis at fine- and coarse-grained scales can be incorporated into flow simulations. Combining such capabilities with a model of dissection of the wall could thus enable a seamless means to study the complexities of the *in vivo* situation.

Given this ultimate goal, in this paper, we present an *orthotropic* DPD model of an artery composed of two integrated constituents: an isotropic elastin-dominated matrix within which is embedded multiple families of stiff collagen fibres. To account for consequences of collagen fibre dispersion and cross-linking we adopt the four-fibre method, which has been shown to describe biaxial data well and is more readily adaptable to discrete modelling than methods using a continuous distribution for dispersion (which ultimately models many discrete fibre directions). The elastin matrix is modelled to approximate a neo-Hookean (nonlinear) elastic behaviour, while equations for the fibre spring potentials are derived directly from the continuum equations of [7,9]. To the best of our knowledge, this is the first attempt to model fibre-reinforced or anisotropic solids using a discrete particle framework. This is also, to the best of our knowledge, the first attempt to model the macroscopic behaviour of the arterial wall using discrete particle methods, which would ultimately allow seamless coupling of discrete models of blood flow within a flexible vessel that can respond to the variations in flow or pressure. We begin with a two-fibre version of the model, which we use for numerical verification, and then we move to the four-fibre model, which we validate with experimental results for two separate biaxial tests that mimic *in vivo* mechanical loading conditions of a vessel. We then present results for delamination of a two-layer arterial wall which is constructed by connecting two separate fibre-reinforced layers.

## 2. Methods

### 2.1. Mathematical models

DPD is a coarse-grained discrete particle simulation method in which particles represent molecular clusters, as first developed by Hoogerbrugge & Koelman [14]. Solids are treated as networks of DPD particles connected by spring potentials. Two model structures are shown in figure 1. The arterial wall is modelled as a mesh of triangulated elastin matrix (grey) with embedded fibres (red/blue) arranged as a diamond mesh representing collagen. DPD particles are located at the vertices of the mesh. The two-fibre model uses two identical fibre families oriented symmetrically to the axis of the cylinder, thus rendering the overall material behaviour orthotropic in response to circumferential and axial loading. The four-fibre model includes additional fibres in the axial and circumferential directions. The DPD particles associated with the fibres will be referred to as fibre particles, while those in elastin matrix will be referred to as matrix particles.

In DPD models of solid material, results converge rapidly as the particle number increases (i.e. as the model is fine-grained). Unless the body undergoes severe shape changes and deformations, which require higher mesh resolutions to resolve small features including sharp corners, very fine-grained meshes are typically not used. For a detailed derivation of the scaling for a triangulated mesh refer to the study of Fedosov *et al*. [18].

Unlike multi-constituent continuum models, the fibres and elastin in the DPD model mechanically interact at the contact points between the fibre particles and the matrix triangle faces, as shown in figure 2. The fibre particle (red circle at *j*_{0} or *j* in the undeformed or deformed configuration) and its three matrix particle neighbours (located at the triangle vertices) are attracted by an elastic force defined in equation (2.3).

The total system energy is
2.1where *U*_{matrix} and *U*_{fibres} are energies in the matrix and fibre, given in equations (2.6) and (2.11). *U*_{int} is an interaction energy from the binding of the two constituents, which is inappropriately absent in most continuum models, but enters naturally in the DPD method [16]
2.2where *N*_{mf} is the number of connection points between the matrix and fibres; in this case, *N*_{mf} is equal to the number of fibre particles because the constituents are connected at points of intersection between the fibre particles and the triangle faces on the matrix layer. *d _{j}* is the distance between the vertex

*j*in the fibres and its corresponding projection point

*j*′ on the matrix layer; is the initial distance (in the unloaded material) between

*j*and

*j*′. Figure 2 shows (

*a*) the undeformed configuration and (

*b*) a deformed configuration.

*k*

_{mf}is the spring constant of the interaction potential. The corresponding elastic force on the vertex

*j*of the fibre layer is 2.3where

*t**is the vector from point*

_{j}*j*to point

*j*′.

Here *U*_{int} is a crucial term in our model for it represents the interaction energy between different constituents. Ideally, this term could vary for arteries from different regions as well as be a function of age, disease status and so forth; however, to compare directly to continuum models, we set *k*_{mf} to be a large positive number to prevent sliding between the fibres and matrix. In this way, *U*_{int} acts as a penalty term in the total energy that significantly restricts the inter-constituent motions and thus models a ‘constrained mixture.’ Specifically, we set the value of *k*_{mf} four orders of magnitude higher than the shear modulus of the elastin layer, *μ*_{0}, which was more than sufficient to prevent different constituents from separating. Further increases in *k*_{mf} did not change the simulated stretch tests (whereas much lower values of *k*_{mf} decreased the overall material stiffness). By considering weaker interaction forces between different constituents, we are able to simulate delaminations as discussed later. First, however, we develop a working constitutive model of a fibre-reinforced artery in DPD that can match results from existing continuum models for a healthy wall.

#### 2.1.1. Material thickness in a two-dimensional model

The arterial wall is technically treated as a two-dimensional material with stresses and strains uniform across the wall (in the radial direction). By assuming local incompressibility, however, we compute local changes in wall thickness, which affects the local structural stiffness. Thus, we avoid a membrane approximation.

The thickness at any location along the wall is calculated from the change in area of the triangular faces in the elastin. Thus, where *h*_{t} is the triangle thickness, *A*_{t} is the triangle area and the subscript ‘0’ denotes the unloaded value. The thickness of a matrix particle, *h*_{m}, is the average thickness of the triangles associated with it. Conversely, fibre particle thickness is computed as the weighted average of the thicknesses of its matrix particle neighbours, which are located at the vertices of triangular faces where fibre particles attach (figure 2). The thickness, *h*_{f} of fibre particle *f*, is computed as
2.4where *h*_{f} is the thickness of the fibre particle and *h _{i}* is the thickness of the matrix particle located at one of the three vertices of a triangle. The weights,

*w*depend on the distance between the matrix and fibre particles and are defined as 2.5where

_{i}*d*is the distance between the fibre particle and matrix particle

_{i}*i*, and

*D*is the sum of the three distances.

### 2.2. Matrix

The total energy in the matrix, *U*_{matrix} is
2.6where *U*_{s,m} is the total spring potential for the matrix, and *U*_{area} is the potential due to area constraints. Spring forces in the wall come from conservative elastic forces of the bonds. Following [15], we combine a wormlike chain (WLC) potential, *U*_{WLC}, with a repulsive power force potential, *U*_{POW}, to achieve a neo-Hookean behaviour. The former potential does not prevent compression, thus the *U*_{POW} potential is added as a repulsive force to resist compression of the springs (e.g. so that the material does not collapse along the width when pulled lengthwise). This approach has also been used effectively in modelling red blood cells [15]. For a mesh with *M*_{s} springs in the matrix, the potential of spring *j* depends on its length *l _{j}*, and

*U*

_{s,m}is the sum of all

*M*

_{s}springs: The spring energies are defined [15] as 2.7and 2.8where

*l*

_{max}is the maximum spring length (equilibrium length is

*l*

_{0}), and

*p*is the persistence length. The POW force coefficient is

*k*

_{p}and

*m*is the exponent. The persistence length and

*k*

_{p}are computed by balancing the spring forces at equilibrium () and from their relation to the shear modulus

*μ*

_{0}(see derivation in supporting material of [15]): 2.9where

*h*and

_{j}*p*are the current thickness and persistence length, respectively, of matrix particle

_{j}*j*. This is a slight variation from the expression in [15]; here, we include the thickness term

*h*, whereas the original expression was derived for a thin-walled material in which thickness did not vary, so

_{j}*μ*

_{0}was treated as the product of the shear modulus and the constant thickness.

Area constraints from [15,18] are defined as
2.10where *N*_{t} is the number of triangles in the mesh, and *k*_{a} is the area constraint coefficient. *A _{j}* is the current area of triangle

*j*, and

*A*

_{0}is the equilibrium value of the triangle area;

*h*

_{t}is the triangle thickness.

### 2.3. Fibre

The fibres are treated without area constraints because individual collagen fibres are ostensibly independent, hence the potential for the fibres comes entirely from the spring forces, Discretizing the continuum model stress function for fibres in [6] gives an equivalent expression for the DPD force of a single fibre spring (see derivation in appendix A):
2.11where *l* and *l*_{0} are the current and unloaded bond lengths, respectively; *d*_{f} is the orthogonal distance between two parallel fibres, and *h* is the thickness of the fibre; thus *d*_{f}*h* is the cross-sectional area of the fibre. The parameter *k*_{1} determines the initial stiffness and *k*_{2} determines the hardening behaviour. The Jacobian is set to *J* = 1 for incompressibility.

### 2.4. Experimental methods and simulation procedures

We first simulate pressurization (radial tractions only) of a thin-walled tube using a two-fibre DPD model and compare results against the theoretical continuum results to verify our implementation. Using the same two-fibre model, we also perform simulations of a planar sheet undergoing biaxial stretching similar to the experiments of Keyes *et al*. [19] and Kural *et al*. [20]. We then consider the more complete four-fibre model in DPD and validate it against experimental biaxial data collected from a common carotid artery (CCA) of the mouse. The advantage of using biaxial stretching tests for calibration is that arteries experience both luminal pressure and significant axial stretch in their *in vivo* state. Thus, biaxial tests reveal more about the *in vivo* mechanical behaviour than do uniaxial tests. For example, a unique mechanical property of arteries is that they maintain a constant level of axial force during pressurization when maintained at their *in vivo* length [21–23]. When extended beyond the *in vivo* axial stretch, increased pressure causes an increase in axial force; when held below the *in vivo* axial stretch, increased pressure causes a decrease in axial force [8]. Below, we describe both experimental and simulation procedures for two biaxial tests: pressure–diameter tests at and around the *in vivo* length (*P*−*d*), and force–length (*f*−*l*) tests at various internal pressure levels.

#### 2.4.1. Pressure–diameter (*P*−*d*) and force–length (*f*−*l*) tests

Full details of our approved (Institutional Animal Care and Use Committee) animal protocols, experimental methods and data analysis have been reported previously [8]. Briefly, murine common carotid arteries were excised, cleaned and placed within a custom computer-controlled biaxial device. Following preconditioning, cyclic *P*−*d* tests were performed over a range of pressures from 0 to 140 mmHg at three different fixed axial stretches, *in vivo* and ±5% of *in vivo*. Similarly, cyclic axial *f*−*l* tests were performed over a range from 0 to 10 mN at four different fixed pressures: 10, 60, 100 and 140 mmHg. Together, these seven protocols provide biaxial data sufficient for robust nonlinear parameter estimation to characterize the highly nonlinear and anisotropic behaviour. Assuming plane stress, theoretical expressions for luminal pressure and axial force can be determined as a function of unknown material parameters for each experimental data point. A nonlinear least-square regression, based on the Levenberg–Marquardt algorithm, is used to estimate the material parameters by minimization of a properly defined objective function. For more details on parameter estimation, refer to [8].

In *P*−*d* test simulations using DPD, the vessel is pressurized over a range of 0–150 mmHg in three separate trials in which it is held at each of three levels of axial stretch, *λ _{z}* = length/unloaded length = 1.64, 1,73, 1.81, where 1.73 is the estimated

*in vivo*axial stretch. Consistent with most continuum modelling, in our DPD simulations, we ignore the cannulated ends of the artery and simulate only the middle segment (1 mm length) where the radius is nearly uniform and where the experimental data are measured. The axial force is computed by summing the material forces on each DPD particle at the ends of the segment. In this way, the entire length of the vessel is free to move radially. The vessel is first extended axially until it reaches one of the three values of

*λ*, at which point it is restricted from further axial deformation for the rest of the simulation so that it can maintain a constant axial stretch. After the axial stretch, pressure is applied in small steps and the resulting tube thickness, diameter and axial force are computed.

_{z}We use a slightly different but equivalent procedure for the *f*−*l* tests. We begin by pressurizing the DPD vessels to one of four levels (10, 60, 100 or 140 mmHg) and allowing the length, diameter and thickness to change freely in response to the pressure. This step is equivalent to the first three steps in the experimental procedure: the artery was extended to a large axial stretch, then it was pressurized and finally it was returned to a lower axial stretch at which axial force was nearly 0. In both procedures, the final state is that the vessel is pressurized and experiences no applied axial force. Because creep, relaxation and hysteresis are negligible during such tests on arteries [8,9,24], the two procedures are equivalent. We confirmed this equivalence in DPD simulations, but chose the first procedure to reduce computational time. From this point, the simulation procedure followed the experiment: an axial force of 0–18 mN was applied to the ends of the tube (distributed uniformly over the circumference).

In contrast with the experiment, the simulation used an equivalent but different method of computing the axial force due to pressure as this approach was more convenient to model: in the experiment, the luminal pressure in the cannulated vessel acted on a closed-ended vessel, which meant that the axial force experienced by the vessel was equal to the sum of *f*_{t}, the axial force applied to the vessel to stretch it axially, plus or the pressure times the inner cross-sectional area of the tube. In DPD, pressure was applied to an open-ended tube, and we corrected for the true axial force of a closed tube after the simulation. Thus, the axial force due to pressure on a closed tube, is subtracted from axial force computed in simulation, where *P* is the pressure used in the simulation (in units of pascals), and *r*_{in} is the current inner radius. Essentially, by simulating pressurization of an open-ended tube, we are doing the equivalent of pressurizing a closed tube and then applying axial force of

#### 2.4.2. Delamination tests

In order to explore applications of the proposed DPD model in dealing with discontinuities arising from failure in soft tissues, we simulate intramural delamination to replicate a peeling experiment in which a rectangular, two-layer strip of arterial tissue is separated by applying opposite forces to the top edges of the layers (as though peeling apart two strips of adhesive tape). We enforced a rigid boundary condition along the short edges of the sheets where the pulling force is applied, mimicking conditions in an experimental testing device in which the mounted specimen would be constrained on each end [25].

Here we use a high binding force value for the fibre–matrix interaction, *k*_{mf} to ensure no slip between these components. Further, since the focus is the inter-layer interactions in an arterial wall, we note that there are mainly two parameters involved: the inter-layer binding strength *k*_{mm} which is similar to *k*_{mf} in equations (2.2) and (2.3), and a new quantity defined as the critical bond length before breaking, *l*_{crit}. Normally, these quantities will not be constant throughout the arterial wall; thus to be more realistic, we assume that *k*_{mm} and *I*_{crit} follow a normal distribution with mean values and standard deviations.

## 3. Results

### 3.1. Pressurization of a thin-walled tube

We pressurized a discrete particle tube embedded with two diagonal fibres symmetric across the axis for verification against continuum results in fig. 7 from [6]. As in the example from [6], we use a cylinder with a mean radius *R* = 4.745 mm and wall thickness *H* = 0.43 mm in the stress-free configuration. We compute the axial and circumferential stretches, *λ _{z}* and

*λ*, respectively, due to internal pressurization using three tubes with different fibre angles,

_{θ}*α*, where

*α*is the reference angle between each diagonal fibre and the axis. Table 1 gives the mesh and the parameters used in the simulation.

Results for the three tubes are shown in figure 3 along with the interpolated results from [6]. The top plot shows the axial stretch due to pressure, while the bottom plot shows circumferential stretch. As the tube is pressurized, causing the circumference to increase, the axial length of the unconstrained tube decreases due to the incompressibility constraint. The tube with the 50° fibre angle experiences the smallest circumferential stretch, because these fibres are closest to being aligned with the circumferential direction, thus they support more of the circumferential load than the 40° or 30° fibres. The 50° fibre angle also experiences the smallest axial shortening (top plot) because the circumferential deformation is less. The DPD model matched well the continuum model for the circumferential stretch. While there is slight discrepancy for the axial stretch, this could be explained by slight differences in surface area constraints and bending stiffness, which are modelled differently in the discrete than the continuum models.

### 3.2. Biaxial mechanical tests

Because arterial tissue is subjected to biaxial loads under normal physiological conditions, it is important to test the model's performance in biaxial mechanical tests. Figure 4 shows the results of a planar biaxial stretch of two 10 mm by 10 mm samples with 0.5 mm thickness, identical except for their fibre angles, 30° and 40°. Table 2 gives the mesh values and parameters. The sheets were deformed with an equal force applied to the horizontal and vertical edges, with free boundary conditions in both *x*- and *y*-directions, thus preserving a rectangular shape. The top plot shows the stretch along the length; the bottom plot shows the stretch along the width. Solid black curves are spectral element (continuum) results based on two-fibre family model of Gasser *et al*. [6], and open circles are DPD results, which match well except at the largest deformations. Compared with the 40° fibre angle, the 30° fibres are closer to being aligned with the vertical (lengthwise) direction, so they support a smaller portion of the horizontal load than the 40° fibres. For this reason, the sheet with the 30° fibre angle experiences significantly higher stretch along the width (bottom plot), thus resulting in increased shortening along the length (upper plot). The extreme deformation at smaller fibre angles is likely the source of discrepancy between DPD and continuum results.

While the two-fibre model provides important verification benchmarks, the four-fibre model [9] has proved better suited to capture the complex, biaxial mechanical response of murine arteries. Here we move to the four-fibre model following [9], as shown in the lower half of figure 1, and present biaxial *P*−*d* and *f*−*l* data for a CCA as described above in Methods. Model parameters for the traditional four-fibre family continuum model and the DPD model are obtained through a best-fit estimation following the procedures outlined in [8]. Fits of the continuum model, in figure 5*d–f*, were very good, with estimated values of The unloaded (stress-free) dimensions of the artery are 5.57 mm length, 447 µm outer diameter and 65 µm wall thickness. Parameters and mesh for the DPD model are given in table 3.

DPD results for *P*−*d* tests are shown in figure 5*a*,*b*, with *R*^{2}-values (>0.93) lower than those for the continuum model, but still very good. Figure 5*a* shows the diameter versus pressure when the tube is held at three levels of axial stretch, and figure 5*b* gives the axial force versus pressure. In figure 5*b*, the results at the *in vivo* axial stretch are relatively flat, so that as the pressure increases, the axial force remains constant. When the vessel is held at an axial stretch above the *in vivo* value, increased pressure causes an increase in the axial force. When the vessel is held slightly below the *in vivo* axial stretch, the axial force decreases with increased pressure. Hence, the DPD model captured this complex, fundamental biaxial behaviour well, similar to most continuum models reported in the literature.

The results for the *f*−*l* tests are shown in figure 5*c* along with the corresponding *R*^{2}-value. These results match the experiment for axial force versus length well, especially at 100 and 140 mmHg. The intersection point of all four curves (roughly 1.73) is the *in vivo* axial stretch value. Overall, the description of data with the DPD model is similar to that of the best continuum models, as exemplified in figure 5*d*–*f* for the same *P*−*d* and *f*−*l* tests.

### 3.3. Modelling delamination of arterial layers

As described in the Methods section, we use a high binding force value for the fibre–matrix interaction, *k*_{mf}, whereas we introduce weaker inter-layer binding strength *k*_{mm} along with a quantity defined as the critical bond length before breaking, *l*_{crit}. Further, we assume that *k*_{mm} and *l*_{crit} follow a normal distribution with mean values and standard deviations given in table 4. This table also lists parameters used for the peeling experiment of figures 6 and 7. We consider two different orientations for the fibres, 60° and 30° relative to the pulling direction (*y*-axis in figure 6).

Figure 6 shows a sequence of snapshots for the peeling process of a double-layered square patch. We plot the peeling force/width versus delamination length in figure 7 for the above-mentioned fibre orientations for five simulations along with their mean curves. The results are qualitatively similar to the experimental results of [25] with the peeling force for 30° orientation higher than that for the 60° orientation. This difference could be due to more aligned fibres along the pulling direction in 30° configuration, which results in higher stiffness of the loaded specimen and an increase in the peeling force.

## 4. Discussion

In this work, we adopted a middle ground between a two-dimensional approximation and a physically more realistic three-dimensional approach. We determined the local thickness of the wall based on volume incompressibility while ignoring transmural variations in stress; this is equivalent to using plane-stress elements in FEM. This approach was critical to our model's ability to produce accurate results for complex biaxial stretching, both in tubular and planar specimens. Further, we believe that future modelling of fluid–structure interactions in DPD will be more efficient using a two-dimensional wall model. All the simulations were performed on a single computer node with 16 cores, and the runs took no longer than several minutes due to the small number of particles needed to represent the arterial walls. We do not expect computation time to become a major challenge even if we use a three-dimensional model of arterial layers. Using DPD to model three-dimensional mechanics of an arterial wall consisting of several layers of particles spanning the wall thickness remains for future study. This approach may be preferred in some cases, such as in studying the mechanobiology and, more importantly, the delamination/dissection process. A truly three-dimensional DPD model of the arterial wall will also require significantly more experimental data, including bond connectivity along the third dimension and the distribution and alignment of fibres across the thickness. There is also less information on ‘bond energies’ in compression than in tension, noting that the continuum stress is compressive in the radial direction. Hence, further extension of the DPD model must await additional experimental data.

For our model, we assumed that the fibres and elastin matrix were perfectly adhered so that no inter-constituent sliding or separation was permitted regardless of the deformation or applied force. This provided a model analogous to its continuum predecessors, which is useful for verification and robust parameter estimation. As seen via comparisons with simulated (planar) and actual (tubular) biaxial data, the DPD model captured the salient results well. The small quantitative difference in the biaxial stretch results may be explained by small discrepancies between the matrix equations in the continuum model and the DPD worm-like chain model. While the two are essentially equivalent at small deformations, they diverge slightly at large deformations. The matrix also introduces some uncertainty in the DPD model as it has three additional parameters that the continuum model does not: *k*_{a}, the local area constraint; *l*_{max}, the maximum WLC spring extension; *m*, the WLC power term. Because *m* had very low sensitivity, we set its value at 2, the same as in [15]. The other parameters were calibrated by optimizing *k*_{a} for a range of values of *l*_{max} and then selecting the best combination. The overall material stiffness could be increased by decreasing *l*_{max} or increasing *k*_{a}, but each choice affected the shape of the response curves differently. For the two-fibre model, we found that the best matches were obtained when *l*_{max} was maximized at 5 (above which, results were unchanged), balanced with a large *k*_{a} optimized to provide the appropriate stiffness. For the four-fibre model, we found the opposite: the best results were obtained by setting *k*_{a} to 0 and decreasing *l*_{max} to 2. The difference may be explained by the way the microstructure of the fibres affects the macroscopic properties of the material. The addition of two new fibre families, in changing the nature of the mechanical response, may alter the sensitivities of the matrix parameters. In the validation experiments with the four-fibre model, the parameter fitting may also have contributed to the error because it was calibrated to the continuum model.

Layer interactions are unavoidable in some pathological conditions. For example, in carotid artery dissection, a leading cause of stroke in young adults [26], intramural layers are separated and are often contiguous with an intimal tear that allows blood within the wall. While the continuum approach poses difficulties in modelling dissections, DPD could provide an effective framework particularly when coupling the wall mechanics and haemodynamics. The present model is suitable to model dissection by adding additional multi-constituent layers and binding them to each other (but using a lower binding force in equation (2.3) to allow a delamination/dissection). By contrast, *ad hoc* models (such as cohesive potential energy) are needed in continuum models and displacement discontinuities have to be resolved by adding more nodal degrees of freedom in the FEM mesh [27].

We performed simple illustrative simulations for the peeling of a two-layered arterial wall to demonstrate one possible use of the DPD model in studying damage and failure in soft tissues, as shown in figure 6. The layer interactions can be modelled (in a discrete sense) by considering harmonic bonds connecting the two layers at the location of matrix particles. Although the peeling force versus delamination length in figure 7 is qualitatively similar to the peeling forces reported in [25], performing a close comparison with that data is not possible due to lack of information on the matrix/fibre properties and fibres orientation. Nonetheless, we examined and presented the sensitivity of the simulations to the values of *k*_{mm} and *l*_{crit} (figure 8). First we increase the mean value of *k*_{mm} twofold to 0.01(N m^{–1}) while keeping its variation at 10%. We observe an increase in the peeling force in figure 8*a* as expected since the interlayer bond strength has been increased. In addition, we test the effect of *l*_{crit} in figure 8*b* by introducing a higher variation through a higher standard deviation of 100%. We find a slight increase in both the amplitude of the peeling force and its mean value. The above analysis suggests that the local microstructure can be effectively modelled in this DPD framework.

## 5. Conclusion

We presented a DPD model of a multi-constituent (fibre-reinforced) artery that mimics its primary microstructural features. To the best of our knowledge, this is the first discrete particle model of such a biosolid. In simple pressurization tests, the two-fibre DPD model demonstrated excellent verification against theoretical results. The four-fibre model also performed well in biaxial stretch tests when validated against experimental results; it was able to reproduce the unique mechanical behaviour of arteries in which, when held at *in vivo* axial length, the axial force nearly remains constant at all pressures, while the axial force rises with pressure above the *in vivo* length and decreases with pressure below the *in vivo* length. Further, we were able to model the delamination process of a multi-layered specimen, which suggests DPD is a promising approach for modelling failure and damage in soft tissues. Most importantly, however, this approach promises seamless particle-based fluid–solid interaction simulations in the future, particularly in cases of localized tissue damage that may result in aneurysms.

## Authors' contributions

A.W. and A.Y. carried out the discrete modelling and simulations, performed analysis, participated in the design of the study, drafted the manuscript and the subsequent revisions; A.Y. carried out the continuum simulations; Z.P. participated in the design of the study and the theoretical model formulation; C.B. carried out the experiments and collaborated on the parameter estimation; J.D.H. helped design the study and helped draft and finalize the manuscript; G.E.K. helped design and coordinate the study and helped draft the manuscript. All authors gave final approval for publication.

## Competing interests

Authors have no competing interests.

## Funding

This work was supported by an NIH grant (1U01HL116323). Computational support was provided by the Center for Computation and Visualization (CCV) at Brown University, ALCF through a DOE/INCITE award, and TACC/STAMPEDE through an XSEDE grant (TG-DMS140007)

## Appendix A. Derivation of fibre spring forces in dissipative particle dynamics

To derive a DPD bond-type equivalent to the stress function in the continuum model of Gasser *et al.* [6], consider stretching a two-dimensional sheet consisting of *n*_{f} parallel fibres in the *x*-direction from initial length *l*_{0} to *l* with an effective continuum cross-sectional area of *A*_{f}. The fibre orientation vector is and the deformed fibre orientation vector is

For the no dispersion case (*κ* = 0), the deformed structure tensor [6] is given as
A 1and the Green–Lagrange strain-like quantity is given as
A 2The stress function in continuum (see table 2 in [6]) is given as
A 3where *k*_{1} and *k*_{2} are two material parameters that determine the initial stiffness and hardening behaviour, respectively.

Thus, the projected Kirchhoff stress tensor for the fibres (eqn 4.8 in [6]) is
A 4where the Kirchhoff stress tensor is related to through Here **p** is the fourth-order projection tensor defined as After some algebra, we find

The Cauchy stress in the *x*-direction is given by
A 5where *J* is the Jacobian, which presents the volume change. If volume is incompressible, *J* = 1.

Since the summation of fibre forces (in DPD or reality) equals the stress times the cross-sectional area (in effective continuum media), i.e.
A 6the final exponential form of the DPD bond force is
A 7where *l* and *l*_{0} are the current and unloaded bond lengths, respectively; *J* = 1 for volume incompressibility.

Also, it may be noted that *A*_{f}/*n*_{f} is the cross-sectional area (thickness, *h* times the width of the sheet) divided by the number of fibres; thus where *d*_{f} is the distance between two parallel fibres along the normal direction, which can be calculated from the initial fibre configuration along with volume incompressibility, as we describe here. First, we assume local volume incompressibility such that the local volume is preserved. Local volume is equal to the fibre bond length *l* times the distance to the nearest parallel fibre, *d*_{f}, times the thickness *h*. Thus,
A 8where the subscript ‘0’ denotes the initial (unloaded) value. From this, we get an expression for *d*_{f}*h*:
A 9

The values *l*_{0} and *h*_{0} are known, and *d*_{f0} can be computed from the initial configuration of the diamond grid
A 10where *α* is the fibre angle and is the distance along the orthogonal direction between two parallel edges of an equilateral diamond with sides having length *l*_{0} and angles 2*α* and 180° − 2*α*. The final expression for d_{f}*h* is
A 11

- Received November 4, 2015.
- Accepted December 18, 2015.

- © 2016 The Author(s)