## Abstract

A three-dimensional cell-based mechanical model of coronary artery tunica media is proposed. The model is composed of spherical cells forming a hexagonal close-packed lattice. Tissue anisotropy is taken into account by varying interaction forces with the direction of intercellular connection. Several cell-centre interaction potentials for repulsion and attraction are considered, including the Hertz contact model and its neo-Hookean extension, the Johnson–Kendall–Roberts model of adhesive contact, and a wormlike chain model. The model is validated against data from *in vitro* uni-axial tension tests performed on dissected strips of tunica media. The wormlike chain potential in combination with the neo-Hookean Hertz contact model produces stress–stretch curves which represent the experimental data very well.

## 1. Introduction

In-stent restenosis is an unwanted growth of a smooth muscle cell (SMC) layer in a coronary artery after treatment of a stenosis with balloon angioplasty and stenting [1,2]. We aim to numerically predict the risk of this dangerous side effect which may lead to further surgeries. Tissue growth is provided by proliferation and migration of SMC to the vessel's lumen [2–4]. As far as migration occurs as a result of intercellular bond degradation and expelling of cells to the lumen, correct modelling of mechanical interactions in the tunica media is crucial for the adequate simulation of the phenomenon. As a starting point for the research on the mechanics of SMC under in-stent restenosis condition, we present a cell-based model of the mechanics of post-mortem tissue and verify the model against published *in vitro* stretch tests.

We model individual SMC as agents which interact with each other via repulsive and attractive forces, while at the same time each cell goes through a cell cycle, taking the cell's mechanical and biochemical environment into account [1]. This model of the tunica media is embedded in an overall multiscale model for in-stent restenosis [2] which has been applied to study many aspects of this pathology [3–7]. So far we only applied our model in two-dimensional simulation of the in-stent restenosis response. We have realized a three-dimensional version (e.g. [8]), and first results of the three-dimensional simulation are in good agreement with porcine data [8]. However, so far the cell-based model of the arterial wall was not validated against known mechanical properties of coronary arteries. That is the main topic of this paper.

Discrete tissue models are generally used in biomechanical and biochemical studies (e.g. [10–16]), including discrete particle dynamics (DPD) models with ‘blobs’ of material as basic model units [16] and cell-based models [10–15] well suitable for tumour growth simulation. We introduce a cell-based model of tunica media of a coronary artery, which we verify against experimental stretch tests [17], where strip samples of tunica media dissected from left anterior descending human coronary arteries were uni-axially stretched in longitudinal and circumferential directions. Here we focus on mimicking the mechanics of quiescent SMC. The cell-cycle transitions and interactions between synthetic SMC, all important for the process of neointima formation in in-stent restenosis, are beyond the scope of this paper. The experimental data published in [18] prove that N-cadherin junctions between synthetic SMC are weakened to provide vessel healing via SMC migration and neointima formation. Such mechanical interactions between synthetic SMC and cell migration modelling are the subject of a separate study.

Taking the single cell at the mesoscale as a starting point for multiscale models of physiology [19] and then, in a ‘middle out fashion’ [20], coupling relevant processes at larger tissue/organ scales or smaller subcellular scales, is now a well-accepted paradigm. In this light we have proposed the Virtual Artery [21] as a multiscale model at the physics–chemistry–biology (PCB) interface [22] where the PCB mesoscale model represents relevant single cell processes in the arterial wall (endothelial cells, SMC) and in the flowing blood (red blood cells, platelets, white blood cells).

The Virtual Artery builds on models of individual cells and their interactions. For instance, we modelled the mechanical properties of single red blood cells and platelets with a boundary element method, and coupled them with flow of blood plasma relying on an immersed boundary method [23,24]. We applied this to study the transport of platelets in aneurisms [25] as well as the shear induced diffusion of red blood cells [26]. Currently we are adding biology and chemistry to capture processes in relation to thrombosis.

The present paper is organized as follows: §2 describes several mathematical models of cell-to-cell interactions considered in this study, next, anisotropic effects modelling is described and, finally, a resulting system of differential equations of cell motion is presented. Section 3 sums up the numerical tests performed in the study; §4 presents simulation results followed by a discussion; and §5 gives the conclusions.

## 2. Mathematical models

### 2.1. Cell-to-cell interactions

A scheme of a healthy coronary artery cross section is shown in figure 1*a*, with spindle-shaped SMC packed in concentric spirals in the tunica media layer.

In our model of the tunica media, we assume SMC to be elastic spheres forming a hexagonal close-packed lattice in the initial, unloaded configuration (figure 1*b*). In the unloaded configuration, a cell (shown green in figure 1*b*) has maximum twelve interaction neighbours: six are in the same pack layer (shown blue) and three neighbours in the upper and lower pack layers (shown with transparent grey). The arterial tunica media model is formed by bending a planar layer of close-packed cells into a tube (figure 1*c*). The cells are uniform, with radii equal to 36 µm, five layers of cells forming the tissue (this is typical for human coronary arteries, e.g. [27]).

Interactions in a pair of cells are assumed to be cell-centred. The balance between repulsive and attractive forces determines the equilibrium configuration of cells. Anisotropy of the tissue is taken into account by modifying stiffness of springs according to their orientation (see §2.2 for details).

Assuming that intercellular repulsion occurs as a result of mechanical compression of elastic cells (modelled as spheres) [28], we apply a neo-Hookean extension of Hertz contact model [29] to describe the repulsive force *F*_{hertz} acting between a pair of cells
2.1where *d* is the indentation computed as mutual penetration of two elastic spheres (in figure 3*a*,*b*), *d* = *R*_{1} + *R*_{2} − *x*, *R*_{1}, *R*_{2} are cells' radii, *x* is the distance between cells' centres, —radius of the contact circle, —effective cell radius, —elastic constant, and *d* ≥ 0.

The interaction potential can be expressed as the integral of (2.1) over indentation . To avoid numerical integration or complex analytical calculation of this integral, a partial sum of Taylor series approximation of (2.1) was used for computing sample deformation energy *Π* in the numerical implementation employing the nonlinear conjugate gradients method (§2.5). Thus, assuming the magnitude of mutual indentations to be smaller than the spheres' effective radius (and hence ), we use a polynomial Taylor approximation for the denominator in (2.1). The approximation is limited by a cubic term:
2.2Using the approximation (2.2) in (2.1) and keeping the terms of the order not higher than 3 when multiplying the square brackets below in (2.3), we obtain an easily integrated polynomial expression for neo-Hookean extension of the Hertz force:
2.3

For indentations smaller than , the relative error of this approximation is less than 2.4%, which is quite good (figure 2 for the comparison of *F*_{hertz}(*d*) and *F*_{hertz_appr}(*d*) in the range of indentations ). In simulations presented below, the mutual indentations never exceeded the value (§4.1).

The neo-Hookean Hertz interaction potential is now approximated as follows:
Besides Hertz repulsion, we also considered a power law repulsion potential proposed in [16]:
2.4where *k*_{2} is a constant parameter.

Next, we assume that intercellular attraction forces in the tunica media are provided by stretch stiffness of protein cytoskeleton and extracellular matrix (ECM) fibres. The intercellular attraction forces are non-zero in a limited zone of negative indentations, where the contact between elastic spheres is open. The intercellular attraction is described with a wormlike chain potential proposed in [16] for a DPD model of an isotropic ‘matrix’ in the arterial wall composition. We apply it on a cell scale to the anisotropic model of tunica media:
2.5where *k _{BT}* =

*k*

_{B}

*Tl*

_{max}/4

*p*,

*k*

_{B}is the Boltzmann constant,

*T*is absolute temperature,

*p*is a persistence length,

*l*

_{max}is the maximal chain length,

*z*=

*x*/

*l*

_{max}is the relative distance of a pair of cells.

By summing the wormlike chain attraction force (2.5) with the neo-Hookean Hertz repulsion force (2.1) and, alternatively, with the power-law repulsion force (2.4), we obtain the tissue models referenced below as ‘wmlc + hertz’ and ‘wmlc + pow’ models.

Besides the ‘wmlc + hertz’ and ‘wmlc + pow’ models, we also studied the adhesive contact interaction described by the Johnson–Kendall–Roberts (JKR) model [29]. Equation (2.6) describes JKR interaction force *F*_{JKR} and potential *U*_{JKR}(*d*) combining Hertz repulsion of isotropic linear elastic spheres and adhesive attraction in the contact circle:
2.6where —effective Young's modulus, , —effective cell radius, *γ* is surface interaction constant; *a* is the contact zone radius related to indentation *d* implicitly:
2.7Expressions (2.6) cover both positive and negative indentations. A critical gap at which the cells separate is defined as ; a critical tension force is [30]. After separation, interaction is resumed when the spheres get into contact (at *d* ≥ 0)—the corresponding hysteresis loops are shown with black arrows in figure 3 presenting cell-to-cell interaction potentials and forces.

To avoid numerical solution of equation (2.7), a polynomial approximation of *F*_{JKR}(*d*) proposed in [30] has been used:
2.8The approximation (2.8) works well for indentations with an order of magnitude comparable to that of *d*_{cr} [30].

The potential wells and interaction forces are shown in figure 3, for the three interaction models studied (‘wmlc + hertz’, ‘wmlc + pow’, JKR). The equilibrium configurations marked with filled circles in figure 3*b* refer to the states of positive indentation of the cells. So, in the absence of external loads, spheroids representing the cells in a tissue are pre-compressed. The red-filled area shows zones of the open contact between the cells, where the ‘wmlc + hertz’ and ‘wmlc + pow’ models mimic the response of ECM and cell cytoskeleton to stretching. The parameter values used for cell-to-cell interaction visualization are presented in table 1.

### 2.2. Model anisotropy

The experimental stretch–stress curves in [17] show significant anisotropy of tunica media, with circumferential stiffness about three times higher than axial stiffness, due to the circumferentially elongated shape of smooth muscle cells. We assume stretch stiffness in the radial direction to be equal to axial stretch stiffness. Anisotropy is taken into account by modifying elastic constants according to the direction of an intercellular bond. Thus, the elastic constant *C* in the Hertz repulsion force (2.1) is computed as follows: *C* = *C _{l}* + (

*C*−

_{c}*C*) · |cos

_{l}*α*|, where

*α*is the angle between the intercellular bond and circumferential direction in the initial (unloaded) configuration,

*C*,

_{l}*C*are elastic constants for the bonds oriented along axial and circumferential directions, respectively. The stiffness parameters

_{c}*k*

_{2}in (2.4),

*k*and

_{BT}*l*

_{max}in (2.5), in (2.6) are scaled in a similar way, depending on the orientation of an intercellular bond towards the circumferential direction.

### 2.3. Residual stresses modelling

Arterial walls are subjected to bending stress *in vivo*. This is confirmed by experiments on measuring opening angles (e.g. [31]). The isolated layers of tunica adventitia and tunica media show the presence of residual bending stresses, as well. A possible explanation of residual bending stresses in tunica media is its heterogeneity across the wall width; higher concentration of ECM fibres at the outer radius makes the outer layer contract more intensively when the ring is cut. In the presented model, the residual bending was simulated by intensifying attraction forces in the outer layer of SMC (see §4.3 for details).

### 2.4. Neighbour detection algorithm

In the non-deformed configuration, each cell in the hexagonal pack has 12 contacting neighbours (figure 1*b*, six in the same layer—the blue cells—and three in both the upper and lower layer—the grey cells). During loading, the neighbour list is dynamically updated at each loading step and repulsive interactions are added if new penetrations occur. In the *in silico* stretch tests, the attractive bonds were never deleted, because the simulated range of stretches (from 1 to 1.4) does not imply any tissue ruptures. According to [17], the ultimate stretch values in tunica media are 1.81 and 1.74 for circumferential and axial tension, correspondingly. The margins have been obtained in [17] for normal samples of tunica media composed of contractile SMC. For synthetic cells however, the margins can be lower. As discussed in §1, we now focus only on the quiescent SMC mechanics.

Note that although the repulsive bonds list is updated dynamically at loading, the attractive bonds list is kept static and is determined by the initial, unloaded configuration of the tissue. Doing so, we mimic the elastomeric properties of the tissue. If the neighbour detection algorithm would entirely be based on positive mutual penetrations in the actual configuration (with no interaction in the open contacts), plastic deformation patterns would be mimicked, better representing amorphous media (as, for example, in [32],where plastic deformation patterns were obtained as a result of neighbour detection specifics).

### 2.5. Solving the equations of motion

A tissue is modelled as a system of discrete particles interacting with their neighbours. The motion of the system is described by Newton's law:
where *m _{i}*,

*a*are mass and acceleration of a particle, and

_{i}*F*is the force vector acting on a particle number

_{i}*i*.

We assume an over-damped motion and fully neglect inertia. Hence, cell accelerations are zero and the force vector is formed from viscous and elastic components. The resulting system of ordinary differential equations describing tissue mechanics is expressed as follows:
2.9where *X = X*(*t*) is the global vector of cell centre positions, *b* is the damping factor and *F*(*X*) is the global nonlinear elastic forces vector.

Cell centre positions are specified in the reference (initial) frame by three Cartesian coordinates per cell and stored in the global vector *X*. The global elastic forces vector is formed in a loop over the cells: based on the computed distances between a cell and its neighbours, the forces in each pair of cells are calculated by formulae (2.1)–(2.4), and then projected onto reference frame axes; the projections are accumulated in the global force vector.

In the general case of finding a non-stationary solution of (2.9), the system is integrated over time using the explicit 4-stage Runge–Kutta scheme [33]. We can also immediately solve for the stationary solution *F*(*X*) = 0 by approaching the equilibrium configuration with a nonlinear conjugate gradients method [34]:
where *s _{i}* are search directions,

*Π*(

*X*) is elastic energy of the system, −∇

*Π*(

*X*) =

_{n}*F*(

*X*),

_{n}*β*is the Polak–Ribière parameter [34].

_{n}## 3. Problem set-up and simulation procedure

Geometric parameters of a stripe sample of tunica media extracted from human left anterior descending coronary artery were taken from [17]. The parameters are presented in table 2. The model value of cell radii was 36 µm, so that there were five layers of cells across the thickness of tunica media.

The *in silico* tests of the tunica media model included the following stages:

(1) Uni-axial stretching of a strip of tunica media and qualitative comparison of simulated stretch–stress data with experimental data from [17]. The comparison has shown that the ‘wmlc + hertz’ model is the most suitable for mimicking the mechanical properties of tunica media (see §4.1). Further numerical tests were performed for the ‘wmlc + hertz’ model only.

(2) Calibrating the ‘wmlc + hertz’ model parameters by matching simulated and experimental curves with the least-squares method. Six parameters in total have been calibrated:

*l*_{max},*k*and_{BT}*C*in circumferential and axial directions.(3) Analysis of sensitivity of stretch–stress curves to variation of model parameters.

(4) Tests on mimicking opening angle formations by introducing heterogeneity into the model: increasing the stiffness of intercellular bonds in the outer layer of SMC.

(5) Simulation of stent deployment, wall stress analysis.

Stress components at stretching were calculated as *σ _{a}*

_{,c}=

*F*

_{a}_{,c}/

*A*, where

*σ*

_{a}_{,c}are the axial and circumferential stress components,

*F*

_{a}_{,c}the axial and circumferential loads, and

*A*the sample cross-sectional area in the deformed configuration. Stretch was calculated as

*λ*=

*L*/

*l*, where

*l*and

*L*are the unloaded and deformed sample length, respectively. The initial configuration of the sample was stress-free.

## 4. Results

### 4.1. Uni-axial stretch tests, model parameters calibration

Deformation of the *in silico* strip sample is shown in figure 4*b*–*e*, for axial and circumferential stretching tests (the initial configuration of the sample is shown in figure 4*a*). The colour palette represents the distribution of longitudinal and circumferential projections of intercellular forces SX and SZ (N). The *in silico* and *in vitro* stress–stretch curves are compared in figure 5*a*,*b*. The experimental curves obtained in [17] show a rapid nonlinear growth of tunica media strength above physiological ranges of stretch; this highly nonlinear response is provided by the ECM and cytoskeleton fibres. As opposed to the experimental data, the ‘JKR’ model produces stretch–stress curves with a saturation zone occurring due to reduction of cellular contact area at tension. At a stretch equal to 1.25, the *in silico* JKR sample ruptured. The ‘wmlc + pow’ and ‘wmlc + hertz’ curves qualitatively agree with the experimental ones. Finally, the ‘wmlc + hertz’ model was chosen as the most appropriate for further calibration because the Hertz contact model has some experimental background in cell mechanics [28]. The ‘wmlc + pow’ model therefore was not calibrated for the precise matching with the experimental data.

Calibration of the ‘wmlc + hertz’ model parameters was performed with the least-squares method. The resulting optimal parameters are presented in table 3, for longitudinal/radial and circumferential directions of transversal anisotropy. The calibrated ‘wmlc + hertz’ stress–stretch curves are shown in figure 5*a*,*b*.

Circumferential and axial cell-to-cell interaction potentials in the calibrated ‘wmlc + hertz’ model are presented in figure 6. Figure 6 shows that the equilibrium indentation for the axial potential and for the circumferential potential; therefore, the Hertz force approximation proposed in §2.1 is relevant and the error of approximation is less than 3%.

The JKR model parameters used for simulation have the following values: *γ* **=** 0.75 × 10^{−3} mN µm^{−1}, ; the wmlc + pow model parameters are as follows: *l*_{max_c}/(2*R*) = *l*_{max_x}/(2*R*) = 1.1, *k _{bT}*

_{_x}=

*k*

_{bT}_{_c}= 7 × 10

^{−2}mN µm

^{−1},

*k*

_{2_x}=

*k*

_{2_c}= 70 mN μm

^{2}.

The big difference between calibrated values of Hertz repulsion coefficients *C _{c}* and

*C*(

_{x}*C*= 23 kPa,

_{c}*C*=

_{x}*C*= 0.24 kPa,

_{r}*C*/

_{c}*C*≈ 100, table 3) in fact simulates the elongated shape of contacting cells in a real tissue. If we refer to the Hertz theory of elliptical contact and consider two contact configurations in a pair of identical ellipsoids with semi-axes of 50 × 5 × 5 micrometres (these values represent typical dimensions of SMC), the contact force will depend on the local Gaussian curvature of the contacting surfaces [30] (figure 7).

_{x}For configuration 1 in figure 7, the Gaussian curvature radius of each ellipsoid is computed as , and this value is then used for computing effective radius in formula (2.1) of the Hertz theory. For configuration 2 in figure 7, the curvature radius equals *b*.

According to formula (2.1), the main term in the expression for the Hertz force is proportional to the following expression: . At a fixed level of stretch *ɛ*, *d* ∼ *bɛ* for case 1 and *d* ∼ *aɛ* for case 2. Therefore, for case 1 and for case 2.

Next, taking into account different amount of contacts per unit cross-sectional area in axial and circumferential directions, we get the following estimates for the total force response produced by a unit area:

From this follows: , which gives us a ratio of 177 for a given ellipsoid dimensions.

As already mentioned above, the elastic constants in the Hertz force scale by a factor of almost 100: *C _{c}* = 23 kPa,

*C*

_{x}*=*

*C*= 0.24 kPa. This factor of 100 agrees in the order of magnitude with the value obtained above from the Hertz theory for homogeneous elastic ellipsoids, proving that the elliptical shape of the cells has been indirectly accounted for in our current model.

_{r}According to experimental data, tunica media is an incompressible material [17]. In discrete models, the bulk incompressibility of arterial wall model is often tuned by adjusting radial stiffness of the tissue across the wall, so that the incompressibility condition *λ _{c}* ·

*λ*

_{r}*λ*= 1 is satisfied (e.g. [16]). However, this procedure does not allow controlling pairwise ratios of strain components, especially the relation between axial and circumferential strains at axial or circumferential tension. We assumed radial spring stiffness to be equal to axial stiffness. Therefore, bulk incompressibility was not simulated. At

_{x}*in silico*stretching, the value of

*λ*·

_{c}*λ*·

_{r}*λ*product varied in the ranges of [1; 1.14] and [1; 1.19] for axial and circumferential stretching, respectively.

_{x}### 4.2. Sensitivity of simulated stretch–stress curves

Sensitivity of stretch–stress curves of the ‘wmlc + hertz’ model against parameter variation is presented in figure 8. Variation of *l*_{max_c} and *l*_{max_x} was performed with a 1% step and those of the other parameters with a 10% step. Variation of maximal chain length for circumferentially oriented springs *l*_{max_c} significantly affects both circumferential and axial stiffness (figure 8*a*), with circumferential stiffness showing a larger response. Variation of maximal chain length for axially oriented springs *l*_{max_x} affects axial stiffness moderately and does not significantly alter circumferential stiffness (figure 8*b*). The difference between sensitivity of the model to variation of circumferential and axial stiffness parameters is caused by the non-equality of these directions in the hexagonal lattice (figure 1*b*). While there are bonds oriented along circumferential directions, there are no bonds strictly oriented in longitudinal direction. A fixed level of relative circumferential extension of the sample causes the same extension of axial bonds, while longitudinal extension of the sample causes both rotations and extensions of intercellular bonds.

### 4.3. Residual stresses and opening angles

*In vivo*, the arterial tunica media, as well as tunica adventitia, are subjected to residual bending stresses [31]. The residual bending stresses have both longitudinal and circumferential components; thus, a ring fragment of tunica media with an axial cut transforms into a near-planar surface with negative curvature (figure 9*a* shows numerical simulation); the curvature of the resulting surface characterizes the level of residual bending stress [31]. Similarly, due to longitudinal residual bending, an initially planar strip of tunica media gets bent (figure 9*b*, for corresponding simulation). To mimic residual bending, we have intensified the attractive interactions in the outermost layer of SMC: this was achieved by increasing the *k _{BT}*

_{_c}and

*k*

_{BT}_{_x}constants for the outermost cells. While mimicking circumferential opening angle formation was possible with a moderate increase of circumferential spring stiffness (up by 20%,

*k*

_{BT}_{_c}= 9.38 × 10

^{−2}mN μm in the outer layer), mimicking longitudinal opening angle formation, on the contrary, required a very high change of axial spring stiffness (approximately by five times,

*k*

_{BT}_{_x}= 9.02 × 10

^{−3}mN μm in the outer layer). Such difference between two directions can be explained by their inequality in the dense hexagonal lattice. Moreover, in reality, axial opening angle formation can be caused by a low-pitch spiral package of smooth muscle spindles straightening when cut. However, the current model does not take into account this spiral build-up. Therefore, the question of modelling axial opening angle is still open and, in the future, may require changing the discrete arterial wall structure presented in figure 1

*c*: at bending a planar layer of close-packed cells into a tube, the edges of the layer should be shifted along each other.

### 4.4. Stent deployment simulation

Stent deployment simulation was performed as a preliminary step for in-stent restenosis simulations [1–6], in order to obtain the initial configuration of the vessel lumen. In the current vessel model, the vessel wall only contained the tunica media layer (figure 10*a*). Radial movement of stent struts was simulated by imposing a kinematic boundary condition on the radial displacements of the cells contacting with the struts (these cells are shown in red in figure 10*a*). The radial displacements of stent struts are equal to 0.1 of the vessel's inner radius which is 1.4 mm. The end sections of the vessel are stress-free. The movement of the vessel as a rigid body is eliminated by the viscous term in the motion system (2.9). The explicit time-stepping Runge–Kutta procedure mentioned in §2.5 was used for this test. In future studies, we plan to include tunica adventitia in the vessel wall model, for modelling the stent deployment stage in the in-stent restenosis analysis. For patients with atherosclerotic and non-atherosclerotic intimal thickening, the tunica intima also becomes a load-bearing layer [17]. Therefore, including the tunica intima into the mechanical vessel model could also be studied.

Circumferential projections of intercellular forces are shown in figure 10*b*,*c*, for two stages of stent deployment, for a radial displacement of stent struts equal to 5% and 10% of the inner radius of the artery respectively.

### 4.5. Generalization of the model to arbitrary cell size

The calibrated parameters presented in table 2 were obtained for a tissue model taking a cell radius of 36 µm. In general, these calibration results can be applied to arbitrary cell sizes using correction coefficients determined so that at a fixed level of relative stretch, the total force response from cells in the strip cross section preserves its value, for any cell size. The number of cells composing a cross section of the sample with a hexagonal close-packed lattice is (figure 1*b*) , where *b* is cross-section width, *h* is tunica media layer width. If the cell radius is scaled down by a factor of *n* (*R*_{new} = *R*/*n*), the number of cells composing a cross section is naturally scaled up by a factor of *n ^{2}.* The Hertz force acting on each cell at the fixed level of applied stretch is scaled down by

*n*

^{2}(see equation (2.1) and note that indentation

*d*in (2.1) is also scaled by

*n*at the fixed level of relative stretch of the sample). The total repulsive response of the sample preserves its value under the given level of stretch. Also, the wormlike chain force on each cell is scaled up by a factor of

*n*due to scaling of

*l*

_{max}(see equation (2.5)). To keep the total attractive response constant, we scale

*F*

_{wmlc}(

*z*) down by

*n*

^{3}.

Therefore, the calibrated parameters from table 2 can be used for arbitrary cell size with the following correction: *k _{BT}*

_{_c_corrected}=

*k*

_{BT}_{_c}/(36 µm/

*R*)

^{3},

*k*

_{BT}_{_x_corrected}=

*k*

_{BT}_{_x}/(36 µm/

*R*)

^{3}. The rest of model parameters do not depend on the cell size.

## 5. Conclusion

A cell-based anisotropic model of tunica media in human coronary artery has been proposed and numerically verified against the *in vitro* uni-axial stretch tests from [17]. The *in silico* stretching tests have shown that the mechanical behaviour of tunica media can be realistically mimicked with a cell-based model using wormlike chain potential (representing compliance of cytoskeleton network and ECM fibres) and neo-Hookean extension of Hertz contact potential (representing repulsive contact interactions between SMC). The JKR adhesive contact model produces stress–stretch curves which are qualitatively different from the curves obtained for elastomeric stretch patterns observed in arterial wall layers.

Regarding the compressive stiffness of tunica media, we have considered the neo-Hookean extension of the Hertz contact model as an appropriate model of intercellular repulsion; however, this assumption needs to be verified against experiments.

Our further research on tunica media mechanics will focus on carrying out a coupled biomechanical/biochemical study of neo-intima formation dynamics in order to assess interaction cut-off margins for synthetic SMC and correctly mimic SMC migration.

## Data accessibility

The source code is publicly available and hosted in the BitBucket directory at https://bitbucket.org/naunat/tunicamedia. The animations of the *in silico* tests (including stretching, open angle formation and stent deployment) are available at the Dryad repository [35].

## Authors' contributions

D.R.H. conceived the original model; N.B.M. improved it and implemented the current version of the model, designed and performed the simulations, analysed the results and drafted the manuscript. A.I.S. analysed the results and helped in drafting the manuscript. A.G.H. conceived of the study, designed the study, coordinated the study and helped in drafting the manuscript. All authors reviewed the manuscript and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

We acknowledge financial support by the Russian Science Foundation, Agreement no. 14-11-00826 (10.07.2014). We also acknowledge partial funding from the European Union Horizon 2020 research and innovation programme under grant agreement no. 671564, the ComPat project (http://www.compat-project.eu/) under grant agreement 675451, and the CompBioMed project (http://www.compbiomed.eu/).

## Acknowledgements

The authors would like to express their gratitude to Prof. Gerhard Holzapfel, Department of Computational Biomechanics, Graz University of Technology, for providing the digital data on the *in vitro* stretch tests of dissected human coronary artery layers.

- Received January 17, 2017.
- Accepted June 5, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.