## Abstract

Arterial wall dynamics arise from the synergy of passive mechano-elastic properties of the vascular tissue and the active contractile behaviour of smooth muscle cells (SMCs) that form the media layer of vessels. We have developed a computational framework that incorporates both these components to account for vascular responses to mechanical and pharmacological stimuli. To validate the proposed framework and demonstrate its potential for testing hypotheses on the pathogenesis of vascular disease, we have employed a number of pharmacological probes that modulate the arterial wall contractile machinery by selectively inhibiting a range of intracellular signalling pathways. Experimental probes used on ring segments from the rabbit central ear artery are: phenylephrine, a selective *α*1-adrenergic receptor agonist that induces vasoconstriction; cyclopiazonic acid (CPA), a specific inhibitor of sarcoplasmic/endoplasmic reticulum Ca^{2+}-ATPase; and ryanodine, a diterpenoid that modulates Ca^{2+} release from the sarcoplasmic reticulum. These interventions were able to delineate the role of membrane versus intracellular signalling, previously identified as main factors in smooth muscle contraction and the generation of vessel tone. Each SMC was modelled by a system of nonlinear differential equations that account for intracellular ionic signalling, and in particular Ca^{2+} dynamics. Cytosolic Ca^{2+} concentrations formed the catalytic input to a cross-bridge kinetics model. Contractile output from these cellular components forms the input to the finite-element model of the arterial rings under isometric conditions that reproduces the experimental conditions. The model does not account for the role of the endothelium, as the nitric oxide production was suppressed by the action of L-NAME, and also due to the absence of shear stress on the arterial ring, as the experimental set-up did not involve flow. Simulations generated by the integrated model closely matched experimental observations qualitatively, as well as quantitatively within a range of physiological parametric values. The model also illustrated how increased intercellular coupling led to smooth muscle coordination and the genesis of vascular tone.

## 1. Introduction

The arterial wall can be viewed as a control system able to regulate blood flow in order to satisfy local tissue oxygenation and nutrition requirements. This function is performed, in part, by spontaneous fluctuations in vascular tone and diameter, known as vasomotion, and is facilitated by the contractile apparatus located within the smooth muscle layer of the arterial wall [1]. Arteries form an anisotropic structure composed of three primary layers that perform distinct functions. The outermost layer, adventitia, is made of a tissue with fibre dispersion which confines the inner arterial structures. Active vascular contractility is governed by the smooth muscle cells (SMCs) located in the second layer, called media. The endothelium, a cellular mono-layer lining the inner surface of the blood vessels, operates as an active interface between the blood flow and the medial layer. Between these structures lie layers of artery-specific connective tissue called elastic laminae [2], which largely convey the mechano-elastic properties of the vascular wall. The active contractile machinery of SMCs, driven by the phosphorylation of the actin–myosin motors, is catalysed by intracellular Ca^{2+}. Under physiological conditions the intracellular Ca^{2+} concentration exhibits modest variations, however, when operative conditions deviate far from normal, such as in the case of injury or under pharmacological interventions, significant fluctuations in intracellular Ca^{2+} concentrations may occur that are reflected in the pronounced dynamical variation in vascular tone [3–5]. The complex structure and function of the arterial wall, even in the absence of blood flow, suggest that a mathematical/computational multi-physics approach is necessary in order to elucidate the dynamical intricacy of the underlying biological system and to address questions that so far evade experimental investigations. To address this problem, a considerable number of multiscale/multi-component models for the arterial wall have been proposed in recent years [6–9]. The development of models accounting for the elastic behaviour of the vascular wall has been based on extensive experimentation on the mechanical properties of vascular tissue under a variety of stress–strain conditions [10–14]. Based on these findings, several methodologies have been proposed in the last decade for simulating the smooth muscle contractility [15–20]. In spite of considerable advances, the active component of vascular contractility, centred on the cellular Ca^{2+} dynamics of the smooth muscle has not been incorporated in a systematic way. This is particularly important as, according to classification by Fischer [21], the smooth muscle responsible for arterial vasomotion can be considered of the ‘fast type’ from a mechano-elastic point of view, and is therefore markedly sensitive to the cellular wall dynamics. Detailed modelling of vasomotion as an expression of multi-channel ionic signalling that regulates arterial smooth muscle Ca^{2+} dynamics, has been proposed in [22,23]. This work was extensively validated against a broad range of pharmacological interventions that specifically inhibit individual transport mechanisms. Extended cellular arrays of coupled SMCs were subsequently used to study the emergence of large-scale synchronization [24]. In this work, we have employed a hybrid model, based on [22,23], to integrate the active contractile behaviour of the media smooth muscle layer with the structural response of the arterial wall. The computational model developed incorporates two distinct scales: the cellular, where cytosolic Ca^{2+} catalyses cross-bridge (CB) kinetics, and the continuous where the contractile units (CUs), and therefore the arterial tissue, exhibit deformation and stress. Evaluation of the CB kinetics at cellular level relies on a modified version of the Hai & Murphy model [25,26]. The cellular network finite-element design follows the anatomical morphology of the tissue considered. For the mechano-elastic characterization of the arterial wall, we follow the work of [17,19,20]. From the structural point of view, the tissue is assumed to be a fibre-reinforced hyper-elastic material [11,27] and incompressibility is enforced by means of a standard penalty method. Simulation results were validated against experimental recordings of vascular tone obtained from rabbit ear arteries under isometric conditions. In the absence of fluid flow, the relationship between stress and deformation of the vascular structure becomes the only stress/deformation-generator mechanism. In addition, we were able to minimize the influence of the endothelium on the smooth muscle contractile apparatus by inhibiting production of nitric oxide (NO), the dominant endothelial control mechanism in the central rabbit ear artery [28]. The performance of the modelling framework has been tested against a number of cellular Ca^{2+} dynamics scenarios, induced by drug interventions able to specifically modulate the SMCs contractile machinery.

## 2. Modelling methodology

The multiscale modelling methodology that accounts for the contractile dynamics of the arterial wall is presented in figure 1, where the flow of information is highlighted. The model developed is based on the coupling of the smooth-muscle contractile apparatus with the mechano-elastic properties of the arterial wall. We outline here the mathematical formulation of the model components and their interface as a requirement for emergent coupled behaviour. At the cellular level, Ca^{2+} dynamics is described by a system of three coupled nonlinear ordinary differential equations which express cytosolic Ca^{2+} concentration (*χ*), Ca^{2+} concentration within the sarcoplasmic reticulum (SR) (*ζ*) and the membrane potential (*η*). Physiological regulation of these variables is due to the coordinated effect of multiple ionic currents that have been incorporated in the formulation of the model and are described in detail in the following section. Cytosolic Ca^{2+} is the primary catalyst for the phosphorylation of CB formation that occurs within each SMC, and which facilitates the contractile response of the cellular cytoskeleton. At the level of the CU, four distinct states can represent the CB kinetics: dephosphorylated myosin (*n*_{M}), phosphorylated myosin (*n*_{Mp}), attached actin–myosin filaments (*n*_{AM}) and attached phosphorylated actin–myosin filaments (*n*_{AMp}). The model assumes a spatial uniformity of CB within the cytosol. The state describing the fractions of attached actin–myosin filaments (represented as *n*_{AM}+*n*_{AMp}), along with variable *λ* that describes the physical deformation of the contractile elements (the stretch ratio in the tissue structure model), serve as inputs to the CU model, where the active free energy function is computed. The calculations performed in the CU model depend on the internal variable , which represents the current relative sliding between the filaments. The derivatives of the active free energy function with respect to the fourth invariant (, ), calculated at the CU subsystem, are then passed directly to the tissue structural mechanics model (with displacement *u* and pressure *p* as variables), where the total (passive and active) tissue stress and deformation is evaluated. This accounts for the arterial wall structural response. Several stages of this physiological process may be probed through targeted pharmacological interventions (e.g. specific ion channel blockers) to reproduce and investigate the origin of vascular disease states. All the model parameters are listed in tables 1, 2 and 4.

### 2.1. Cellular level

#### 2.1.1. Cytosolic Ca^{2+} homeostasis

The fundamental processes involved in the regulation of cytosolic Ca^{2+} concentrations are outlined in figure 2. All the ionic pathway contributions depicted have been incorporated in a model that combines the approaches previously proposed in [22,23]. In particular, we employed a system of three ordinary differential equations, as in [23], but maintained the predominance of ryanodine-sensitive over InsP3-sensitive intracellular stores, as appropriate for the rabbit ear arteries used in the experimental study [22]. Moreover, the model constructed maintains the well-established view that smooth muscle Ca^{2+} dynamics are due to the interplay of a membrane and an intracellular oscillator [30]. The mechanism underlying the intracellular Ca^{2+} oscillator is Ca^{2+} induced-Ca^{2+} release (CICR) via the ryanodine receptor. Once released into the cytosol, Ca^{2+} is re-sequestered by the stores through the sarcoplasmic/endoplasmic reticulum ATPase (SERCA), thus completing the cycle. By contrast, the membrane oscillator is centred on cellular polarization and voltage-operated Ca^{2+} channels (VOCCs). Membrane potential is mainly influenced by the balance between K^{+} and Cl^{−} gradients, and by Na^{+}-Ca^{2+} exchange (NCX), depending on whether the exchanger operates in forward or reverse mode. The current model also includes the store operated Ca^{2+} channel (SOCC), which links directly store content with cytosolic Ca^{2+} influx [31]. This addition allows for testing theoretically the role of this channel in specific pharmacological intervention scenarios. An explicit delineation of the modelling components is presented below.

##### 2.1.1.1. Ionic channels

The extracellular influx is the sum of the concentration currents *Φ*_{A}, *Φ*_{S}, *Φ*_{V} and *Φ*_{N}, which represent the Ca^{2+} permeable non-selective cation channels (NSCC), SOCC, voltage-operated cation channels (VOCC) and reverse mode Na^{+}-Ca^{2+} exchange (NCX), respectively. In this study, the first current is assumed to be constant, while the others depends on *χ* and *η* via the following equations:
2.1
2.2
2.3where *A*_{S}, *ζ*_{S}, *E*_{Ca}, *z*_{Ca1}, *R*_{Ca}, *E*_{NCX}, *x*_{NCX} and *z*_{NCX} are cellular model parameters. The SR acts as an inner store, uptaking cytosolic Ca^{2+} by means of the SERCA pump (*Φ*_{B}(*χ*)) and releasing it into the cytosol through ryanodine-sensitive SR-Ca^{2+} release channel (*Φ*_{C}(*χ*, *ζ*)). Store leakage (*Φ*_{L}(*χ*, *ζ*)) is also accounted for. These fluxes are described through the following expressions:
2.4
2.5
2.6where *B*_{SR}, *n*_{SR}, *x*_{SR}, *C*_{Ry}, *p*_{Ry}, *m*_{Ry}, *x*_{Ry}, *y*_{Ry} and *L*_{SR} are cellular model constants. The Ca^{2+} extrusion from the cytosol by ATPase pump (*Φ*_{D}(*χ*, *η*)) is modelled as follows:
2.7where *D*_{EX}, *k*_{EX}, *z*_{EX} and *R*_{EX} are cellular model parameters. The chloride (*Φ*_{Cl}(*χ*, *η*)) and potassium (*Φ*_{K}(*χ*, *η*)) ion membrane fluxes are modelled as follows:
2.8and
2.9where *E*_{Cl}, *x*_{Cl}, *z*_{Cl}, *E*_{K}, *z*_{K}, *β*_{K} and *R*_{K} are cellular model constants.

##### 2.1.1.2. Intercellular communication

In a cellular cluster, each cell is able to communicate with its neighbours by exchanging Ca^{2+} ion (*J*_{Ca}) and voltage (*J*_{V}) currents. Assuming that the cellular size (ratio between the volume and surface) is uniform along the grid, we can define the net fluxes and exchanged by the *i*th cell as follows:
2.10and
2.11where nNeigh is the number of adjacent elements, and and are the intercellular diffusion coefficients.

##### 2.1.1.3. Global cellular balance

The system variables *χ*, *ζ* and *η* evolve in time according to the following system of nonlinear ordinary differential equations:
2.12
2.13
2.14where *γ* is a scaling factor relating the net movement of ion fluxes to the membrane potential.

#### 2.1.2. Cross-bridges kinetics

The kinetics is described through four different states representing the fraction of CB that are: (i) attached and dephosphorylated (*n*_{AM}), (ii) attached and phosphorylated (*n*_{AMp}), (iii) detached and dephosphorylated (*n*_{A}) and (iv) detached and phosphorylated (*n*_{Ap}). The temporal evolution of the four-state kinetics shown in figure 3 is described through the following system of ordinary differential equations:
2.15
2.16
2.17
2.18where *τ*_{1}, *τ*_{2}, *τ*_{3}, *τ*_{4} and *τ*_{5} are the kinetic rates. Among these, *τ*_{1}, depends on the intracellular Ca^{2+} concentration (*χ*) via [22]:
2.19where *τ*_{0} and *χ*_{0} are the material constants. As *n*_{M} + *n*_{Mp} + *n*_{AMp} + *n*_{AM}=1, one of the states (*n*_{M} in this work) can be assumed dependent on the other three states, leading to a three independent variables system.

### 2.2. Continuous-level model

#### 2.2.1. Contractile unit mechanics

To model the force development within the CU, the work of Murtada *et al.* [17,19,20] is followed. Both phosphorylated and de-phosphorylated attached CB (*n*_{AMp}, *n*_{AM}) are considered elastic with the same mechanical stiffness. As mentioned previously, *λ*, the current stretch of the CU, can also be defined as the ratio between the current and the reference CU length. The average elastic elongation of the attached CB () can thus be calculated as follows:
2.20We note that both (representing the relative fibre sliding) and are normalized with respect to the reference CU length and are taken to be negative for contraction. The relative actin–myosin filament sliding in the CU can be driven either by the myosin power-stroke or the external force/deformation. This internal variable is thus decomposed into a chemical () and mechanical () component. The temporal evolution of the chemical component can be derived from the following force balance:
2.21where *P*_{c} is the stress associated with the driving force from the cross-bridges while *β*_{a} and *α*_{a} are fitting parameters. The internal driving stress *P*_{c} depends on the contraction/relaxation state of the CU, ie,

2.22

where *κ*_{AMp} is a parameter related to the force of a power-stroke of a single CB and *κ*_{AM} is related to the force-bearing capacity of a dephosphorylated CB during muscle extension. The energy stored in the CU is related to the filament sliding resistance from the surrounding matrix (*P*_{a}), which can also be seen as the (averaged) first Piola–Kirchhoff stress over the CU:
2.23where *μ*_{a} behaves like an active shear modulus and defines the relative filament overlap as a parabolic function of :
2.24where and are material parameters.

#### 2.2.2. Tissue structure

From the structural point of view, the medial tissue is considered as a hyper-elastic fibre-reinforced material [11,27], with the fibres aligned along the circumferential direction [20]. Various representations of the arterial tissue have been previously proposed, including models of transverse fibre distribution at a diagonal angle [11,27,32]. These approaches were specifically proposed to accurately represent the adventitia layer of the vessel wall. As we are primarily focusing on the media layer, in this study, we have followed the work proposed in [20] where the media layer fibres are aligned only in the circumferential direction. The free energy function (*Ψ*) is split into volumetric (*Ψ*_{vol}) and isochoric components; the latter is then decomposed into active (), accounting for the CU chemo-mechanics, and passive parts (), i.e.
2.25The volumetric part is assumed to be proportional to the energy potential as
2.26where *κ* and *J* are, respectively, the penalty number and the determinant of the deformation tensor **F** (*J* = det **F**). Both and depend on the CU stretch (*λ*), that can also be related to the fourth invariant as
2.27where is the deviatoric part of the right Cauchy deformation tensor and **a**_{0} is the direction of the media unstressed fibre. The active component depends directly on the CB attached (*n*_{AMp}, *n*_{AM}), i.e.
2.28The passive part is modelled as a classical anisotropic material with one fibre aligned with the SMCs as
2.29where *μ*_{p}, *c*_{p1} and *c*_{p2} are material constants. More details regarding the solid mechanics formulation are provided in electronic supplementary material, S1.

### 2.3. Multiscale coupling and solution procedure

#### 2.3.1. Space discretization

Each level of the framework is discretized by a spatial grid. The cellular network grid reflects the morphology of the tissue, so that each element represents a ‘real’ cell. For the continuous level, a finite-element discretization of the domain was carried out. By assigning each SMC to each finite-element the connectivity of the cellular grid and the mesh coincide. The variables computed at the cellular level (i.e. Ca^{2+} dynamics and CB kinetics) may be considered as internal variables in the finite-element framework. The CU variables (, *λ*, etc.) are evaluated at the Gauss integration points of the finite-element, in order to take into account the spatial diversity of deformation over the element. For the structural problem, staggered finite-elements are used in which the displacement field is interpolated linearly while the pressure and dilation coefficient are constant over each element. The nonlinear problem is solved via a classical Newton–Raphson procedure.

#### 2.3.2. Time integration

The models/subsystems constituting the framework are solved in a block segregated fashion, as depicted in figure 1. As there is no feedback between the subsystems, it is possible to employ a different and optimal time integration strategy for each of them. Thus, the Ca^{2+} dynamics is solved by an explicit and adaptive scheme (Runge–Kutta Merson), whereas the time-dependent equations for both the CB kinetics and CU mechanics are solved by the forward Euler method (as in [20]). Note that for computing the CU mechanics the deformed configuration (expressed in terms of *λ*) of the previous time step is used. By comparison, the tissue mechanics problem is solved in an implicit manner. The methodology proposed is valid for either quasi-static or dynamic problems, depending on whether the case considered is under isometric or non-isometric conditions. The same time step was employed for all subsystems. The solution procedure is presented step by step in table 3. The complete model has been implemented into an in-house C++ code. Visualization software ParaView [33] was used for post-processing analysis. The finite-element mesh has been realized by means of the three-dimensional mesh generator software Gmsh [34].

## 3. Experimental study

The pharmacological studies employed in the validation of the modelling methodology along with the model set-up developed to reproduce the experimental settings are reported below.

### 3.1. Experimental protocol

Isolated rabbit ears were obtained as described previously [1], the central ear artery was removed and cleaned of adherent fat and connective tissue. To measure force, 2 mm wide rings were mounted on 0.25 mm diameter steel hooks in a myograph (model 610M, Danish Myotechnology, Aarhus, Denmark) containing oxygenated (95% O_{2}; 5% CO_{2}) Holman's buffer (composition in mM: NaCl 120, KCl 5, NaH_{2}PO_{4} 1.3, NaHCO_{3} 25, CaCl_{2} 2.5, glucose 11 and sucrose 10) at 37.0°C. Prior to any pharmacological interventions, the rings were maintained at a resting tension of 1 mN over a 60 min equilibration period, with frequent readjustments in baseline tension to correct for stress relaxation. The average inner and outer diameters of the annular segments were approximately 0.7 and 0.8 mm, respectively (figure 4*a*). Following the loading phase, the arterial rings obtained the deformed configuration shown in figure 5. By considering a Cartesian reference system, the loading is applied along the *x*-direction. Preparations were incubated for 30 min with both the endothelial nitric oxide synthase inhibitor N(G)-Nitro-l-arginine methyl ester (l-NAME, 300 μM) and the cyclooxygenase inhibitor indomethacin (10 μM) to inhibit prostanoid formation. Rings were then constricted with phenylephrine (PE, 1 μM) and, once constrictor responses had reached a stable plateau, cumulative concentration–response curves to 10 and 30 μM CPA, and 10 and 30 μM ryanodine were obtained.

### 3.2. Model settings

From a structural point of view, each ring is assumed to deform symmetrically with respect to the *x*-plane, while no translations along *y* and in the longitudinal direction (*z*) are expected. The system can thus be reduced to one-eighth of the ring. A finite-element mesh consisting of 5750 linear hexahedral elements is used in the calculations. The global force *F* applied to the hooks is evaluated by integrating the traction contribution of the reduced arterial system over the global contact surface. For the cellular cluster, each cell was associated with one element. At subcellular level, we associate a CU to each Gauss integration point. To account for cellular variability (e.g. size, rates of Ca^{2+} uptake/extrusion) the term corresponding to the influx via Ca^{2+} permeable non-selective cation channels (*Φ*_{A}) is randomized with a normal distribution (mean value given in table 4 for the relevant pharmacological simulations, and s.d. = 0.1 μM s^{− 1}). Randomization of *Φ*_{A} creates a population of cells exhibiting variable oscillatory frequency. This allows us to study the effect of increasing strengths of intercellular coupling as a factor of smooth muscle synchronization and generation of vascular tone. For all simulations, the initial values for *χ*, *ζ* and *η* are set equal to 0.1 μM, 0.2 μM and −0.02 V, respectively. The reference set of kinetic rates (*τ*_{0}, *τ*_{2}, *τ*_{3}, *τ*_{4}, *τ*_{5}) necessary for solving the CB dynamics are taken from Parthimos *et al*. [22].

## 4. Results

### 4.1. Cellular coupling conditions

The proposed analysis is carried out for varying levels of cellular coupling in order to establish the dependency between the diffusion coefficients *α*_{C}, *α*_{V} and the global Ca^{2+} dynamics. These simulations are carried out for an unloaded ring configuration. The variables associated with the Ca^{2+} dynamics are monitored for three different cells: cell1, cell2 and cell3 (figure 4*b*). In figure 6, the time evolution of *χ* for (*α*_{C},*α*_{V}) equal to 0.0, 0.1, 1.0 s^{−1} is shown. The plot shows clearly that coupling does not affect significantly the amplitude of *χ* signal. The cellular coupling tends to synchronize the *χ* beating pattern along the cluster.

The spatial distribution of *χ* for two coupling levels (weakly coupled: (*α*_{C}, *α*_{V}) = 0.1 s^{−1}, strongly coupled: (*α*_{C}, *α*_{V}) = 1.0 s^{−1}) at two different time instants (*t* = 0.1 s, *t* = 6.0 s) are shown in figure 7. It is evident that coupling promotes the formation of travelling waves along the annular domain.

### 4.2. Framework validation

A range of pharmacological interventions associated with the activation/inhibition of specific cellular mechanisms was employed for the validation of the model. To mimic the diffusion of pharmacological agents, parametric changes were applied in a graded fashion to all the cells constituting the network. The interventions selected (i.e. phenylephrine, CPA and ryanodine) are associated with the modulation of the cellular Ca^{2+} homeostasis which is reflected in the contractile state of the smooth muscle. The actions of these pharmacological probes were simulated by the gradual variations of the associated model parameters, reported in table 4. The actin and myosin filaments are assumed to be detached for *t* = 0 s (*n*_{M} = 0.5 and *n*_{Mp} = 0.5), while the initial is set equal to 0 for each CU. Cells are assumed to be strongly coupled with (*α*_{C},*α*_{V}) set equal to 1.0 s^{−1}. The simulated drug interventions were performed only after the Ca^{2+} and CB variables reached stationary conditions. For each plotted result the drug intervention occurred at *t* = 0 s unless otherwise stated.

As the ring is clamped at the hooks, the resultant sum between internal and external forces at the nodes in contact with the hooks must be null. In this study, we are mainly interested in the global force developed at the hooks, which can be seen as the sum of nodal contributions along the contact surface. In addition to this, we compute also the force developed locally at cell1, which is assumed to be proportional to the sum *n*_{AMp}+*n*_{AM} [22]. For a more complete view of how the contractile force emerges from Ca^{2+} dynamics at the cellular level, see electronic supplementary material, S2.

#### 4.2.1. Phenylephrine intervention

Phenylephrine is a selective agonist of *α*1-adrenergic receptor, associated with vasoconstriction. The action of phenylephrine at cellular level was modelled by increasing the cytosolic Ca^{2+} influx via NSCCs. A more complete picture of the action of phenylephrine involves an initial Ca^{2+} release from the SR, via inositol 1,4,5-trisphosphate-sensitive Ca^{2+} release channels (IP3R channels), followed by sustained Ca^{2+} influx into the cytosol through NSCCs. To simulate the action of phenylephrine, the coefficient *Φ*_{A} of each individual cell was increased by a factor of 7. Therefore, the mean value of the *Φ*_{A} distribution increased from 0.6 to 4.2 μM s^{−1}. Two different drug dilution times (Δ*t*_{dil} = 50 and 100 s) were used, as shown in figure 8. By increasing cytosolic Ca^{2+} uptake, *χ* in cell1 starts to oscillate periodically. This pattern is reflected also in the store Ca^{2+} concentration variations (figure 8). As expected, higher Δ*t*_{dil} involves longer transient before reaching stationary conditions.

The forces generated by this intervention, normalized with respect to the force developed at the beginning of the intervention, are presented in figure 9, where three responses are shown for different dilution times (Δ*t*_{dil}) and different cellular coupling (*α*_{C}, *α*_{V}). By comparison to the simulated results, it appears that the drug was able to activate the muscle tissue at a dilution time Δ*t*_{dil} ∼ 50 s. We also observe that the final magnitude of the simulated and measured force is very similar.

#### 4.2.2. Cyclopiazonic acid intervention

Cyclopiazonic acid (CPA) is an inhibitor of the SERCA pump, preventing refilling of the store and is thus associated with Ca^{2+} store depletion. The effect of CPA, in terms of modulating the intracellular Ca^{2+} oscillator, is highly concentration-dependent as shown previously [35]. CPA was administered through concentration increases from 10 to 30 μM. In a few experiments CPA concentrations of 100 μM were used, which accelerated the tapering effect already in evidence (see electronic supplementary material, S2). For simulation purposes, the action of CPA was reproduced by linearly decreasing coefficient *B*_{SR} from 400 to 350 μM in 2000 s (while the randomized coefficient *Φ*_{A} was maintained constant for each cell, with a mean distribution value of 0.8 μM s^{−1}). In both the experiment and simulations CPA caused a small reduction in the oscillatory amplitude, while a regular waveform was maintained throughout the simulated intervention (figure 10). A comparison between the experimental and simulated time series is shown in figures 11*a*,*b*. Note that the oscillatory amplitude of a SMC weakly coupled to its neighbours is reduced, when compared with the same cell under strongly coupled conditions (figure 11*a*(ii)(iii)). This effect of weak cell–cell coupling is greatly enhanced in the case of force development, as seen in figure 11*b*. This is due to the fact that the force contributions of individual cells remain largely uncoordinated. Global force development is restored under strong intercellular coupling conditions (figure 11*b*(iii)).

#### 4.2.3. Ryanodine intervention

Ryanodine was administrated at the arterial sample at concentrations of 10 and 30 μM. As discussed previously, the action of ryanodine can be simulated in different ways depending on the concentration [22,23]. This is due to the complex multistage configuration of the ryanodine receptor tetramer. For the concentrations of ryanodine used in this study, it is accepted that the compound will block Ca^{2+} release from the SR in a concentration-related fashion. The action of ryanodine was simulated according to [22], by linearly decreasing the coefficient *C*_{Ry} at *t* = 100 s from 1250 μM s^{−1} down to 312.5 μM s^{−1}. The temporal evolution of variables *χ* and *ζ* within the reference cell1 is shown in figure 12 for different cellular coupling conditions ((*α*_{C}, *α*_{V}) = 0.1 s^{−1} and (*α*_{C}, *α*_{V}) = 1.0 s^{−1}). Gradual attenuation of Ca^{2+} release via the ryanodine receptor channels is associated with a decrease in frequency.

The forces obtained from the model are compared against the experimental values in figure 13*a*,*b*. In this case, coefficient *θ* was equal to 0.5. The pattern of the forces follows Ca^{2+} concentrations, with the same gradual decreasing magnitude. To simulate this aspect of the experimental traces, we needed to employ the term in equation (2.1) associated with store-operated Ca^{2+} entry, which is triggered in response to levels of SR Ca^{2+}. Note that this mechanism only had a minor effect in the simulation of the phenylephrine and CPA responses. Similar to the CPA simulations, weakly coupled SMCs remain uncoordinated in the presence of ryanodine, producing weak global contractile force (figure 13*b*). A consistent vascular tone is restored when strong intercellular coupling is maintained.

## 5. Discussion

We have developed a multi-component mathematical modelling framework that accounts for the structural response of the arterial wall under the active contraction of the media smooth muscle layer. The methodology was validated against a set of pharmacological interventions that modulate the contractile apparatus at the cellular level. The multiscale modelling approach combines dynamics and mechanics occurring at different levels, each requiring a specific solution strategy. With regard to the mechano-elastic component, the study was performed under isometric conditions, and thus the inertial force is neglected, allowing us to deal with a simplified system in which the evaluation of a number of dynamic parameters (such as density) was not necessary. Moreover, the ability to choose a finite-element discretization that conforms to the cellular grid, eliminated the need for spatial and temporal interpolation between the two subsystems. As a consequence, it was possible to adopt the same time step for all elements of the model. This is not a limitation, as different integration steps can be selected for each model level if required by the specific problem. This strategy can be implemented in conjunction with interpolation techniques that allow information transmission between the various contributing systems. Although necessary in many cases, this approach can result in loss of accuracy. All experiments in this study were performed following administration of l-NAME to eliminate the inhibitory effect of endothelium-derived nitric oxide on the smooth muscle contractile apparatus [36]. The involvement of secondary endothelium produced electrochemical factors in the contractile activity of the arterial wall has not been included, and should form the basis for further elaboration of the current model [37]. Several arterial ring experiments under isometric conditions have been carried out, each involving four rings ≈ 2 mm in length excised from the same animal. In spite of efforts to maintain identical conditions, variability was present, particularly between individual animals. With respect to the four arterial rings from the same arterial sample (e.g. figure 9), experimental and biological variability that could not be eliminated was due to differences in length, diameter (all rings were from adjacent sections, but there is, inevitably, a reduction in diameter as you descend the vessel), and purely initial condition considerations, such as Ca^{2+} and ionic uptake levels, which have been shown both theoretically and experimentally to greatly affect the oscillatory response of the arterial wall [22,23,30]. Owing to these factors, each arterial ring is at a different initial contractile state, reflected in variable unloaded geometry. We therefore elected to normalize the force traces against the pre-intervention steady state (e.g. the pre phenylephrine force in figure 9) which would allow us to reproduce the percentile contractile response during pharmacological probing. This approach successfully facilitated the main focus of this work which was the integration of the smooth muscle-based contractile apparatus and the mechano-elastic properties of the arterial wall. To investigate this fundamental interaction in the genesis of arterial tone, we have employed a series of experimental studies that probe distinct aspects of these mechanisms [22,23]. By simulating experimental findings with the vasoconstrictor phenylephrine, we were able to match the response time of the drug, and demonstrate the direct link between cellular Ca^{2+} dynamics and the development of force at the tissue level. Indeed, the cellular events occur at the same timescale as the global tissue contraction. This finding is correlated to the increased cytosolic Ca^{2+} concentrations associated with the administration of phenylephrine. This observation highlights the central role of cytosolic Ca^{2+} influx in the generation of vascular tone and the onset of oscillatory activity observed in both experiments and simulations. Indeed, we have previously shown that the levels of Ca^{2+} influx can modulate other ionic signalling pathways in a specific, clinically relevant manner [22]. As noted earlier, phenylephrine mediates Ca^{2+} release through the IP3R channel pathway. We have previously employed different approaches to account for this effect. In [22], we have considered two independent CICR and InsP3-induced Ca^{2+} (IICR) release pathways and assumed that the IICR results in a constant Ca^{2+} release current, which is essentially lumped in the overall Ca^{2+} uptake by the cell. Both effects are therefore accounted for by the same constant *Φ*_{A}. We have shown that such an assumption still allowed the mathematical model of vasomotion to accurately reproduce an extensive range of pharmacological interventions [22]. Further work to distinguish the CICR from the IICR mechanism reproduced the same results with no increased accuracy. We have therefore assumed here the simplified scenario proposed in [22]. The potential of Ca^{2+} influx to determine the natural frequency of the cell's contractile apparatus was the main reason that parameter *Φ*_{A} was selected for randomization across the cellular population. By comparison, cellular dynamics are considerably less sensitive to variations of other system parameters, associated with alternative ionic fluxes [29]. Blockage of the SERCA pump with CPA resulted in a modest pattern variation that highlights the robustness of the store-refilling mechanism. The resilience of Ca^{2+} dynamics under CPA was reproduced by the simulations. A noteworthy aspect of the function of CPA is its ability to transform arterial vasomotion in a controlled manner that follows hallmark transition routes out of chaotic behaviour [35,38]. Evidence of this behaviour is shown in figure 11*a*,*b* although detailed investigation of nonlinear oscillatory transitions was not an aim of this work. Ryanodine receptor dysregulation is implicated in a range of neuromuscular disorders and arrhythmogenesis in cardiovascular diseases [3]. This is mainly due to the complex inter- and intra-subunit interactions within the ryanodine receptor homotetramer [4,5]. Considering the intricate multistage dynamics of the ryanodine channel, computational simulation work can elucidate some of the dominant components of the mechanism. Although studied here in isolation, dysregulation of the ryanodine channel has been associated with upregulation of the SERCA pump protein, to allow for the maintenance of a level of sustainable homeostasis [39,40]. Such findings support the proposed methodology as a testing ground for hypotheses on the pathogenesis of vascular disease.

## Data accessibility

Further experimental data and numerical results supporting this article have been uploaded as part of electronic supplementary material, S2.

## Authors' contributions

A.C. implemented and carried out all the calculations. A.A. provided support to A.C. on implementation. D.H.E. designed and performed all experiments. P.N. and D.P. advised on model improvement, problem suggestions and interpretation of results. All authors developed the ideas and contributed to the writing and editing of the manuscript. All the authors gave their final approval for publication.

## Competing interests

We have no competing interests.

## Funding

The authors acknowledge the financial support provided by the S*ê*r Cymru National Research Network in Advanced Engineering and Materials and the EPSRC (grant number: EP/P018912/1).

## Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.3988731.v1.

- Received October 4, 2017.
- Accepted January 12, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.