Coronary Smooth Muscle Cell Calcium Dynamics: Effects of Bifurcation Angle on Atheroprone Conditions

This work investigates the effect of arterial bifurcation angulation on atherosclerosis development through in-silico simulations of coupled cell dynamics. The computational model presented here combines cellular pathways, fluid dynamics, and physiologically-realistic vessel geometries as observed in the human vasculature. The coupled cells model includes endothelial cells (ECs) and smooth muscle cells (SMCs) with ion dynamics, hetero and homotypic coupling, as well as electro-diffusive coupling. Three arterial bifurcation surface models were used in the coupled cells simulations. All three simulations showed propagating waves of Ca²⁺ in both the SMC and EC layers, following the introduction of a luminal agonist, in this case ATP. Immediately following the introduction of ATP concentration Ca²⁺ waves propagate from the area of high ATP toward the areas of low ATP concentration, forming complex patterns where waves interact with eachother, collide and fade. These dynamic phenomena are repeated with a series of waves of slower velocity. The underlying motivation of this research was to examine the macro-scale phenomena, given that the characteristic length scales of atherosclerotic plaques are much larger than a single cell. The micro-scale dynamics were modeled on macro-scale arterial bifurcation surfaces containing over one million cells. The results of the simulations presented here suggest that susceptibility to atherosclerosis development depends on the bifurcation angulation. In conjunction with findings reported in the literature, the simulation results demonstrate that arterial bifurcations containing wider angles have a more prominent influence on the coupled cells pathways associated with the development of atherosclerosis, by means of disturbed flow and lower SMC Ca²⁺ concentrations. The discussion of the results considers the findings of this research within the context of the potential link between information transport through frequency encoding of Ca²⁺ wave dynamics and development of atheroprone conditions.

1. Introduction

Over the past 30 years, since the pioneering research of Caro et al. (1971), a substantial body of work has been carried out in order to elucidate the relationship between vascular geometry and the growth of atherosclerotic plaques. A recent review by Heo et al. (2014) shows the current views of flow-mediated influence on the development of atherosclerosis. Much has been discovered—however little has been attempted in order to explain the role of spatially varying agonists on cellular dynamics in combination with the cellular pathways located in the intima of the vessel wall, mediated by luminal agonists, and disturbed fluid flow.

1.1. Atherogenisis and Vascular Geometry

Multiple studies agree that disturbed vascular flow is the potential trigger for both endothelial and smooth muscle cell dysfunction in human vasculature, and lead to the activation of pro-inflammatory phenotype and coronary plaque formation (Libby, 2002; Davignon and Ganz, 2004; Plank et al., 2006). Similarly, a study by David (2003) pinpoints wall shear stress (WSS) and disturbed flow in wide-angle bifurcations as a mechanical level trigger of atherosclerosis. However, the underlying biochemical mechanisms of atherogenesis, potentially triggered by spatially varying agonist concentrations in the boundary layer flow remain unclear.

Further, it has been shown that the localization of plaque formation is associated with properties of arterial geometries especially in or close to the coronary and carotid bifurcations (Thomas et al., 2005; Lee et al., 2008; Dong et al., 2014). The review by Chiu and Chien (2011) demonstrates that the majority of plaque formation cases occur close to bifurcations and, more specifically on the outer walls of those bifurcations. Lee et al. (2008) have conducted a number of computational fluid dynamics (CFD) simulations on carotid bifurcation geometries derived from magnetic resonance imaging (MRI) data and found significant relationships between disturbed flow, proximal area ratio and bifurcation tortuosity, whilst Dong et al. (2014) concluded that wider bifurcation angles, specifically in the left circumflex (LCX) bifurcation shoulder, are likely sites of disturbed flow.

1.2. The Frequency-Mediated Protein Phosporylation Hypothesis

In contrast to the view that [Ca²⁺] (calcium ion concentration) is the mediator for protein phosphorylation the work of Goldbeter et al. (1990) suggested that the [Ca²⁺] oscillatory frequency encodes for phosphorylation of proteins related to atherogenesis, where the phosphorylation kinetic rates were a function of [Ca²⁺] frequency. Hence one could hypothesize that the kinetic rates may be a marker for predisposition to kinase activation if the Ca²⁺ dynamics varied sufficiently in areas prone to plaque formation.

Multiple works indicate the possibility that cellular calcium dynamics, in terms of variations in frequency of oscillations, could enhance pathways that increase the probability of atherosclerotic plaque formation. For example, the mitogen-activated protein kinase (MAPK) pathway has been shown to be a crucial signaling mechanism in foam cell formation, endothelial cell (EC) activation, and vascular smooth muscle cell (SMC) proliferation (Hakala et al., 2006; Elliott and Eguchi, 2015). Our hypothesis is that Ca²⁺ dynamics could encode for MAPK pathway activation and that the resulting oscillatory [Ca²⁺] variations as a function of geometry may be a biomarker for the early initiation of plaque formation.

1.3. In-silico Studies of Ca²⁺ Dynamics

Variations in both blood flow, the associated wall shear stress and luminal surface agonist concentrations occur over length scales much larger than a single cell. A number of models were used to investigate the [Ca²⁺] dynamics phenomena occurring via coupled cells. However, none of them provided simulations which occur over the macro scale where characteristic lengths are centimeters rather than micrometers. Similarly, the studies reported to date included a limited number of pathways modeling hetero and homocellular communication. For example, Koenigsberger et al. (2006) modeled a vessel segment of 50 micrometers in length consisting of 80 SMCs and 30 ECs where coupling consisted of assuming a diffusive process for membrane potential. The experiments simulate pressure variation by including a stretch channel, but did not include any alternative pathways such as Nitric Oxide (NO). The work of Jacobsen et al. (2008) reported experiments with 45 coupled SMCs without an endothelial layer. The report by Koenigsberger et al. (2010) describes a model of a 2D tissue strip consisting of 130x3 SMCs.

A small number of studies report macro-scale simulations of arterial geometries which examine the coupling of intracellular chemical species. For example, the work of Shaikh et al. (2012) used a massively parallel one-dimensional (1D) model of coupled cells to investigate the propagation of [Ca²⁺] waves with spatially varying agonist concentrations. The results demonstrated that by varying the coupling coefficients, significantly different wave propagation profiles were obtained including propagation upstream with respect to blood flow. Further, the work of Boileau et al. (2015) provided a computational framework of cellular dynamics over large vessel segments with an emphasis on realistic arterial geometries generated from patient data. However, this work did not include a model of the endothelium, nor that of the myoendothelial gap junctions linking SMCs and ECs, otherwise known as heterotypic coupling. The experiments by Boileau et al. (2015) modeled membrane potential using a linear gradient term, rather than modeling physiologically realistic coupling which includes both Fickian and electro-diffusive terms.

1.4. Coupled Cells Dynamics Simulations in Arterial Bifurcations

The model presented here combines cellular chemical pathways, fluid dynamics, and physiologically realistic vessel geometries. The micro-scale dynamics, including chemical species concentrations within individual ECs and SMCs were modeled on arterial bifurcation surfaces containing over one million cells. The coupled cells model includes EC and SMC ion dynamics, hetero and homotypic coupling, as well as Fickian and electro-diffusive models. The underlying motivation of this research was to examine the macro-scale phenomena, given that atherosclerotic plaques have characteristic length scales much larger than a single cell in the arterial wall.

In addition, the inclusion of CFD-derived agonist concentrations as a function of vessel geometry provides a state-of-the-art investigation into the relationship between spatially varying agonist concentrations, ionic dynamics, and vessel geometry. The model and the in-silico experiments presented in this work offer the opportunity to further the understanding of the complex relationship between atherosclerosis and vessel geometry.

The remainder of this article discusses the enhancements to the coupled cells model and simulation framework in section 2, followed by the simulations results in section 3. Section 4 considers the findings of this research within the context of the potential link between information transport through frequency encoding of Ca²⁺ wave dynamics and development of atheroprone conditions.

2. Methods

The coupled cells experiments presented in this work bring together a series of methodological enhancements which embody a step change in the area of in-silico atheroprone conditions modeling. These enhancements include physiologically realistic arterial geometries, an IP₃ pathway and electro-diffusive coupling in the coupled cells model. While the foundation of the experimental framework is provided in the report by Zakkaroff et al. (2016), this section focuses on the details of the set of bifurcating arterial segments and the corresponding agonist maps along with the refined coupled cells model.

2.1. Bifurcating Arterial Surfaces

Three varying angles of bifurcation were chosen to represent the variation in human population: 50°, 80° and 110°, based on the physiological range of arterial bifurcation angle variation reported in the work of Dong et al. (2014). Figure 1 shows the three surface models along with the agonist map overlays, where each bifurcation model is approximately 18 mm in length and 2.5 mm in diameter. An arterial vessel surface of this scale may include over one million ECs and SMCs combined.

FIGURE 1

Figure 1. Arterial bifurcations and CFD-based ATP maps. The time-averaged ATP maps for the 50° (Top), 80° (Middle), and 110° (Bottom) arterial segments bifurcation generated using CFD. High ATP concentrations are observed at the apical regions, while the lateral seams of the arterial segments show low ATP concentrations. The large-angle bifurcations correspond to higher ATP concentrations in the areas of flow detachment over the bifurcation seams.

The method of creating EC and SMC surface meshes corresponding to endothelial and smooth muscle layers presented in Zakkaroff et al. (2016) was extended to enable the creation of physiologically accurate bifurcating arterial segments with varying angles between daughter branches.

In addition to the varying angles between branches, these surface meshes feature two additional geometric properties: first, they are non-planar—the vessels considered in this research belong to a coronary tree, which is naturally affected by curvature as the vessels wrap around the myocardial surface; second, the radii of the daughter branches are related by Murray's law:

$\begin{array}{ll}{r}_p^3=r_{d_1}^3+r_{d_2}^3,& (1)\end{array}$

which relates the radii of downstream branches to the main stem (Murray, 1926), where r_p is the radius of the parent branch while r_d1 and r_d2 are the radii of the respective daughter branches.

2.2. Agonist Maps

Time-averaged agonist maps for the arterial bifurcation surfaces are shown in Figure 1. As the agonist, adenosine triphosphate (ATP), concentration was used as an input to the cellular model via the endothelium, special care was given to deriving a physiologically realistic spatial distribution for the concentration maps. Luminal receptor agonist concentration maps (concentration per unit area) were derived from steady-state CFD simulations as described in the work of Zakkaroff et al. (2016), where the surface ATP concentrations were derived from flowing ATP concentrations and a reactive state at the luminal wall. Due to the low diffusivity of agonists in flowing blood, time-dependent profiles have a negligible difference to that of steady state.

The agonist maps show relatively high ATP concentrations at the bifurcation apex, low concentrations at the lateral bifurcation seams, and near-uniform concentrations in all other regions. The work of Wan et al. (2008) provided ATP concentrations as a function of shear stress and found ATP concentrations between 0.2 and 1 μM, while the results reported by Gorman et al. (2007) found the concentration of ATP in blood to be approximately 0.2 μM. Variations of ATP in this research show similar concentrations to that of both Gorman et al. (2007) and Wan et al. (2008).

The agonist maps presented here link the geometric properties of a given arterial segment, via spatially varying WSS on the endothelium as shown by David (2003), with the cellular dynamics described by our physiological model. We note that while ATP is treated in this model as an example of an agonist, other mechano-transduction pathways exist, but they are not included in this model.

2.3. The Coupled Cells Model

The coupled cells model presented in Zakkaroff et al. (2016) describes the binding of the agonist (ATP) to the luminal surface of the endothelium, which triggers a cascade of reactions and movements of ionic chemical species in both ECs and SMCs. Figure 2 provides an overview of the coupled cells model. However, the G-protein pathway responsible for IP₃ generation via ATP binding to the luminal endothelial surface was reduced to a linear relationship between ATP and IP₃ concentrations, which was a simplification in need of refinement. In contrast, this research considers the IP₃ pathway in greater detail by building on the two models of Lemon et al. (2003) and Bennett et al. (2005).

FIGURE 2

Figure 2. Coupled cells model. Schematic representation of mass transfer dynamics in a single coupled EC/SMC unit. (1) ATP binds to the purinergic P2Y receptors on the EC surface, activating the G-protein-coupled receptors (GPCR) which then activates a membrane-bound enzyme, phospholipase C (PLC), in the presence of Ca²⁺. PLC activation results in the phosphorylation of phosphatidylinositol 4,5-bisphosphate (PIP₂). (2) PIP₂ gives rise to IP₃ which is then released in the cytosol. IP₃ binds to IP₃ receptors on the endoplasmic reticulum (ER)/sarcoplasmic reticulum (SR) surface. (3) IP₃-bound IP₃-receptor induces release of Ca²⁺ from the ER/SR into the cytosol. (4) Ca²⁺ release from intracellular store sensitizes the IP₃-receptor further which releases more Ca²⁺, thus making a Ca²⁺ rich domain in the cytosol in both ECs and SMCs, which depolarizes the membrane potential. This is know as calcium-induced calcium release (CICR). (5) The sarco/endoplasmic reticulum Ca²⁺-ATPase (SERCA) pump requires ATP to move Ca²⁺ against the concentration gradient. Cytosolic Ca²⁺ encourages the replenishment of the intracellular stores via this pathway. (6) Ca²⁺ leaks from ER/SR consistently under a concentration gradient between cytosolic and ER/SR luminal Ca²⁺ and keeps the Ca²⁺ in equilibrium during a non-stimulated state of the cell. (7) The EC cytosolic Ca²⁺ favors the influx of extracellular Ca²⁺ from non-selective cation channels. (8) The SERCA pump pushes out cytosolic Ca²⁺ back into the SR. (9, 10) Activation of the KCa channel in the EC, upon binding to Ca²⁺ ions intracellularly at BKCa and SKCa, allows potassium ions (K⁺) to move out of the cytosol. This hyperpolarizes the membrane potential. (11) Residual current, predominantly consisting of monovalent ions, contributes to the membrane potential repolarization. (12) The IP₃ concentration increases in the SMC cytosol via transmission of IP₃ from coupled ECs. This IP₃ increase activates the downstream IP₃-induced Ca²⁺ release. (13) The IP₃ induced and CICR Ca²⁺ depolarizes the membrane potential. (14) The membrane depolarization results in the influx of Ca²⁺ from voltage-operated calcium channels (VOCCs) which will close upon repolarization in the following steps. (15) Ca²⁺, in addition to other pathways, is pushed out via sodium ion Na⁺/Ca²⁺ exchanger. (16) Binding of Ca²⁺ ions to K_Ca opens BKCa channels in the SMC causing K⁺ ionic efflux and membrane repolarization. (17) Influx of chloride ions (Cl⁻) adds to the repolarization. (18) The medium for intercellular communication via homocellular gap junctions can either be Ca²⁺, IP₃ or membrane potential gradients. (19) Heterocellular gap junctions can couple an EC to an SMC via Ca²⁺, IP₃ or membrane potential coupling. Hyperpolarized EC membrane potential can hyperpolarize the SMC plasma membrane and consequently close VOCCs.

2.3.1. IP₃ Pathway

The IP₃ dynamics in the coupled cells model used in this research, corresponding to the component 1 in Figure 2, were adapted from the work of Lemon et al. (2003) and are described as follows:

$\begin{array}{ll}\frac{dĨ}{dt}=ϵ·r_h·PI{P}_{2,tot}-J_{degrade}+\mathrm{\text{coupling terms}},& (2)\end{array}$

where $ϵ=u_c·{(N_a·V_{EC})}^{-1}$ is a constant to ensure dimensional equivalence; ${PI{P}_2}_{tot}=5.0\times 10^7$ which is used in place of the variable PIP₂ concentration, which typically shows little variance, as noted in the work of Johny et al. (2017); the coupling terms are the hetero and homoocellular diffusion rates of IP₃ between neighboring ECs and SMCs.

The hydrolysis rate of PIP₂ induced by fully activated PLC was modeled as follows:

$\begin{array}{ll}{r}_h=\alpha \left(\frac{\tilde{c}}{\tilde{c}+K_{Ca}}\right){\tilde{G}}_{prot},& (3)\end{array}$

where K_Ca is the dissociation constant for Ca²⁺ binding to PLC, $\tilde{c}$ is the cytosolic concentration of Ca²⁺ in the EC, and α is an effective signal gain parameter which was fitted to experimental data, and is treated as a constant.

The dynamics of activated G-protein are described as:

$\begin{array}{ll}\frac{d{\tilde{G}}_{prot}}{dt}=k_a·(\delta +{\rho }_{P2Y})·(G_{pro{t}_{tot}}-\widetilde{G_{prot}})-k_d·{\tilde{G}}_{prot},& (4)\end{array}$

where k_a and k_d are the activation and deactivation rates of G-protein respectively; G_prottot is the total number of G-protein molecules within the EC; and δ is the G-protein intrinsic activity parameter. Finally, the ratio of bound to total P2Y receptors (under the assumption of fast kinetics), taken from Bennett et al. (2005), is as follows:

$\begin{array}{ll}{\rho }_{P2Y}=\frac{\left[ATP\right]}{K_{ATP}+\left[ATP\right]},& (5)\end{array}$

where K_ATP is the effective Michaelis-Menten constant for ATP binding to a P2Y receptor and [ATP] is the lumenal ATP concentration. Complete definitions of the ODEs and fluxes in our model, including parameter values, are provided in the Supplementary Material.

2.3.2. Electro-Diffusive Coupling

The work of Koenigsberger et al. (2005) offered a simplified hetero and homocellular coupling model based on an assumption that membrane potential is diffused between adjacent cells. A more accurate coupled cell model needs to consider that membrane potential is a phenomenon that exists due to the imbalance of charge between two separated bodies resulting in an electrical gradient, and hence cannot be diffused. Instead, ionic species carrying charge are diffused, which alters the membrane potential through electro-diffusion, in addition to concentration gradient-based Fickian diffusion.

Jacobsen et al. (2007) described the existence of gap junction currents explained by both membrane potential and concentration gradients between neighboring cells. This was modeled in the presented research as follows:

$\begin{array}{cc}{I}_{gap,\varphi }=P_x\sigma A_{segment}F\left[\nabla [\varphi ]+\left(\frac{zF}{RT}\overline{[\varphi ]}\nabla V_m\right)\right].& (6)\end{array}$

Here $\overline{\left[\varphi \right]}$ = ([ϕ]_a+[ϕ]_b)/2, which is the average concentration of ion ϕ between two cells a and b; $\nabla V_m\simeq \frac{V_{m,a}-V_{m,b}}{\delta x}$, which is the membrane potential gradient; $\nabla \left[\varphi \right]\simeq \frac{{\left[\varphi \right]}_a-{\left[\varphi \right]}_b}{\delta x}$, which is the concentration gradient of ion ϕ; A_segment is the interface area of neighboring cells; and δx is the distance between cells. The concentration gradient here is only relevant for Ca²⁺ in our model, as intracellular concentrations of all other ions were assumed to remain constant on the considered time scale. The values of these constant concentrations are taken from Kapela et al. (2008). The ion ϕ has the valence z; F is Faraday's Constant; T is the absolute temperature; P_x is the permeability of ion ϕ; and R is the universal gas constant. The currents described by Equation (6) were added to the hetero and homocellular coupling terms for membrane potential. The definitions, values, and dimensions of these parameters are provided in the Supplementary Material. IP₃ coupling due to its zero valence is coupled using Fickian diffusion.

2.4. Coupling Coefficients

The macroscopic mean conductance G in membrane potential coupling has been estimated by Van Rijen et al. (1997) to be 30 nS, whereas Lidington et al. (2000) estimated the macroscopic homocellular electrical resistance in ECs to be 3 MΩ (or $G=\frac{1}{3M\Omega }=333nS$). As $g=\frac{G}{C_m}$, the results of Van Rijen et al. (1997) yield a conductance one order of magnitude less than achieved by the evaluation of G from Lidington et al. (2000). Therefore, in the interest of coherent parameter use, the value of G = 30 nS was used here, as chosen in the work of Koenigsberger et al. (2005).

The membrane capacitance value between ECs is of the order of 30 pF as given by Schuster et al. (2003) and for SMCs it is 10 pF from Parthimos et al. (1999). Thus, for homocellular electrical coupling between SMCs, Yamamoto et al. (2001) calculated the macroscopic gap junctional resistance to be of the order of 90 MΩ (or $G=\frac{1}{90M\Omega }=11nS$). Therefore, the value of homocellular electrical coupling coefficients g and $\tilde{g}$ used for SMC and EC respectively, was 1000 s⁻¹.

Experimental reports on the single channel conductance values of Cx40 and Cx43 gap junctions is scarce, and to the best of our knowledge the data on a cell's macroscopic or net conductance of intercellular Ca²⁺ and IP₃ transfer is not available. The work of Koenigsberger et al. (2006) treat these values as free parameters and set the coefficients such that homocellular Ca²⁺ coupling between SMCs is able to synchronize the Ca²⁺ oscillations of five coupled SMCs in conjunction with the electrical coupling. SMCs are weakly coupled, and hence homocellular coupling coefficient was set to 0.05. EC homocellular Ca²⁺ coupling coefficient was also set to 0.05. In the case of heterocellular coupling, ECs and SMCs are coupled via the same three media.

Macroscopic intercellular resistance across myoendothelial junction has been estimated by Yamamoto et al. (2001) to be 900 MΩ, with the net capacitance between the two cells being approximately 20 pF, therefore making the heterocellular electrical coupling coefficient 50 s⁻¹. For Ca²⁺ and IP₃ transport across the myoendothelial junction, information on permeability is scarce and Koenigsberger et al. (2005) set their coupling coefficient to 0.05s⁻¹ for IP₃. The authors in Kapela et al. (2009) reduced the heterocellular IP₃ coupling coefficient by assuming that myoendothelial coupling was weak, using a gap junction resistance of 13.5 GΩ, however this came from porcine coronary arteries. The variation in unitary conductance is large, judging from the current publications.

2.4.1. NO, EDH, and EET Considerations

Fast [Ca²⁺] changes in ECs and SMCs occur due to the exchange of Ca²⁺ in the ER/SR and the cytosol. On the one hand, the characteristic timescales for this are of the order of a second to several seconds. On the other hand, kinetic parameters that mediate cytosolic Ca²⁺ are of the order of 50 to 100 s, as noted in the work of Yang et al. (2005). We therefore note that the NO pathway is relatively slow and does not play an important part in wave propagation. Thus, although nitric oxide (NO) is a recognized antithrombic vasodilator, it is speculated here that due to the difference in the rates of reaction kinetics it does not influence the Ca²⁺ dynamics investigated in this work. A further examination of the NO-induced influences on cellular function other than cytosolic Ca²⁺ oscillations, is offered in section 4. In addition, it is recognized here that endothelial derived hyperpolarizing factor (EDH) may decrease Ca²⁺ levels, but experimental evidence on this pathway is still controversial.

The work of Hadfield et al. (2013) showed that epoxyeicosatrienoic acid (EET) production follows that of Ca²⁺. If one assumes that a similar pathway exists in the EC and SMC then the EET dynamics would induce potassium efflux, close the VOCC and thence reduce [Ca²⁺] with a possible reduction in [Ca²⁺] oscillations. Potentially, this may have the effect of inhibiting the MAPK pathway, and possibly inhibiting the initiation of plaque growth. However, the model of Hadfield et al. (2013) did not include a VOCC. The model presented in this research does not include an EET influence and it is still unknown as to the nominal value of EET produced in the coronary vasculature. Similarly, it is still unclear exactly how EET influences the BK channel. Thus, the role of EET in the Ca²⁺ dynamics is left out of this work, as it has yet to be fully investigated in the further work.

2.5. Parallel Implementation

To enable parallel simulations of the coupled cells dynamics the arterial bifurcation surfaces were subdivided into quadrilateral domains. The quadrilaterals were grouped into three 2D surface-meshes, each representing a tubular segment of a bifurcation. Within each segment, the quadrilaterals corresponded to computational domains that were mapped to distinct tasks/cores in the Message Passing Interface (MPI)-enabled parallel environment. A periodic boundary condition was imposed along the longitudinal edge of each 2D grid. In this case, the periodic boundary condition refers to the communication-specific mapping of MPI tasks/domains. Each domain included 208 ECs and 80 SMCs, coming to 288 cells per quadrilateral, with a total of 1,175 million cells, spread across 4080 quadrilateral domains. The parallel implementation of the coupled cells simulations presented here is based on the work of Zakkaroff et al. (2016).

The system of ODEs simulating the time-dependent transfer of ions and membrane potential within each quadrilateral domain was solved at a specified time-step. For the simulations presented here the inter-domain communication interval was set to 0.01 s. This interval was chosen on the basis of a number of tests to maintain the balance between MPI communication and computation times. At regular time intervals the state variables from the cells along the edges of the quadrilateral domains were passed to the adjacent cells in the neighboring domains. The inter-domain communication was implemented through the MPI library, the solution of which follows the mesh ghost cell communication pattern recommended for discretised domain decomposition (Gropp et al., 2014). The solution for the ODEs was obtained by using the numeric differential equation solver ARKode, from the Suite of Nonlinear and Differential/Algebraic Equation Solvers (SUNDIALS) (Hindmarsh et al., 2005).

The simulations were performed on an IBM POWER7 architecture, each simulation using 120 cores, with each core solving 34 quadrilateral domains. Approximately 20 h of wall clock time was required to complete each of the 500-physiological-second simulations. Values of state variables, such as Ca²⁺ concentrations and membrane potential, for all ECs and SMCs were recorded after every elapsed physiological second. During post-processing, the output data were mapped onto the surface meshes and visualized and analyzed in ParaView (Ayachit, 2015).

3. Results

Snapshots of Ca²⁺ concentrations in SMCs and ECs for the bifurcation angle of 110° are presented in Figures 3A,B for t = 500 physiological seconds. The results demonstrated that ATP-induced production of IP₃ in the endothelium and its subsequent diffusion into the SMC induced cytosolic Ca²⁺ oscillations for certain ATP concentrations, where the frequency of oscillation is a monotonic increasing function of the production rate of endothelial IP₃. Given that IP₃ is determined by the spatially varying ATP concentration at the endothelial wall, SMCs can oscillate as a function of space. However, a greater appreciation of the complexity of this phenomenon, observable in the form of slow-moving waves both up and downstream through the arterial wall, is gained by viewing the animated visualizations available on the University of Canterbury High Performance Computing YouTube channel¹, which show cellular dynamics over the whole of each simulation.

FIGURE 3

Figure 3. Simulation at t = 500 physiological seconds of Ca²⁺ dynamics in SMCs and ECs using a 110° geometry. In the SMCs (A), waves of Ca²⁺ with maximal amplitude of 1 μM propagated outwards from the apical region. These waves struggled to penetrate into the lateral seams where the SMC Ca²⁺ concentration remained below 0.1 μM. Similar behavior was observed in the ECs (B).

3.1. Wave Patterns

The waves of propagating SMC Ca²⁺ oscillations were observable in all simulations. Immediately following the introduction of the ATP concentration Ca²⁺ waves propagate from the area of high ATP concentration—the apex of the bifurcation. The initial wave front moves relatively quickly in a circumferential direction until it reaches the parent artery stem where it spreads out, reaching the upstream edge in less than 15 s. The wave moves to the outer wall where ATP concentration is low and forms a circular pattern, reducing in radius as time progresses. These dynamic phenomena are repeated with a series of waves slower in velocity. The period of the oscillations emanating from the apex of the bifurcation was approximately 20 s.

As time increases the uppermost parts of the circular waves in the outer wall become unstable and form complex patterns within areas of low ATP concentration, similar to the phenomenon of re-entrant waves. At approximately 350 physiological seconds into the simulation the low ATP concentration area has become highly complex as smooth Ca²⁺ waves move downstream along the inner wall of the daughter arteries, and the outer walls produce elaborate wave structures similar to that found in low ATP concentration regions. Waves also move upstream toward the inlet with a leading edge positioned at the top of the artery.

The Ca²⁺ waves produced at the apex were rapid and distinct, while those at the lateral bifurcation seams propagated more slowly, sporadically, and over significantly shorter distances. Larger bifurcation angles appeared to accentuate these observations, most notably the Ca²⁺ dynamics at the lateral seams of the 110° bifurcation geometry. Head-on collisions of Ca²⁺ waves result in the destruction of both waves.

Finally, the area affected by this complex Ca²⁺ behavior reduced in size as a function of wider bifurcation angles, seemingly as the low ATP concentration area (lateral seam as shown in Figure 1) reduces in size as the angulation increases. In both SMCs and ECs [Ca²⁺] oscillates with a slight phase lag.

3.2. Temporal Averaging and Surface Mesh Unwrapping

There exist significant differences in time-scales between the cardiac pulse, Ca²⁺ wave frequency, and atherosclerotic plaque formation. The temporal oscillatory behavior of EC and SMC Ca²⁺ concentrations and their influence on protein phosphorylation is complex and it is somewhat difficult, if not impossible, to describe it quantitatively in the absence of dynamic visualizations. Temporal averaging was used with the view of examining the influence of the dynamic behavior in the areas of high and low ATP concentration on MAPK pathway. Temporal averaging was applied for each cell x in the surface arterial mesh as follows:

$\begin{array}{ll}{\left[{\overline{Ca}}^{2+}\right]}_x=\frac{1}{(t_f-t_i)}{\displaystyle \sum _{t=t_i}^{t_f}}{\left[C{a}^{2+}\right]}_x(t).& (7)\end{array}$

Here ${\left[{\overline{Ca}}^{2+}\right]}_x$ is the averaged Ca²⁺ concentration in cell x between times t_i and t_f, and ${\left[C{a}^{2+}\right]}_x(t)$ is the Ca²⁺ concentration in cell x at time t. Thus, the maximum/minimum values of oscillations were averaged to provide a more interpretable account of cellular dynamics. For ease of comparison, the differently angled bifurcations were mapped to a 2D surface, to allow a single perspective to capture [Ca²⁺] dynamics over the entire arterial segment. The circumferential edges on these surfaces are periodic and thus, this re-mapping can be thought of as an unwrapping of the original bifurcations. The result of applying these post-processing techniques is shown in Figure 4 for SMC Ca²⁺ values. Temporal averaging was performed over the final 200 s of the simulation such that t_i = 300 s and t_f = 500 s. These values were chosen in order to exclude transient behavior, while averaging Ca²⁺ concentrations over a sufficiently long period.

FIGURE 4

Figure 4. Time-averaged SMC Ca²⁺ concentration over 200 physiological seconds using flat mesh re-mappings of 50° (Left), 80° (Center), and 110° (Right) bifurcation geometries. We observed high averaged Ca²⁺ concentrations at the apical regions and low averaged Ca²⁺ concentrations at the lateral seams in all three geometries. Further, greater disturbances in Ca²⁺ wave propagation occurred at the lateral seams of geometries with wider bifurcation angles.

The results confirmed the observation that greater disturbances in Ca²⁺ wave propagation occurred at the lateral seams of bifurcations with increasing bifurcation angle. Further, this temporal averaging disclosed that geometries with larger bifurcation angles were subject to lower Ca²⁺ concentrations at the bifurcation apex. The areas of averaged Ca²⁺ reduce in size as a function of increasing bifurcation angle. This is primarily due to the influence of fluid dynamics with the growing intensity of the rotational flow caused by the increase in curvature of the bifurcation Comerford et al. (2008) especially in the lateral region.

3.3. Protein Phosphorylation

Our hypothesis has been, on the basis of Goldbeter's work, that frquency encoding in the form of [Ca²⁺] oscillations caused by cellular dynamics, mediated by geometry and spatial profile of the agonist, may influence the MAPK pathway, and hence induce plaque formation. Using the same notation as Goldbeter et al. (1990), the rate of change of the fraction of phosporylated protein, W*, can be written as:

$\begin{array}{ll}\frac{d{W}^{*}}{dt}=\frac{v_p}{W_T}\left[\frac{v_k}{v_p}\frac{1-W^{*}}{K_1+1-W^{*}}-\frac{W^{*}}{K_2+W^{*}}\right],& (8)\end{array}$

where

$\begin{array}{ll}{v}_k=V_{mk}\frac{C{a}^{2+}(t)}{K_a+C{a}^{2+}(t)}& (9)\end{array}$

Here the components in the brackets correspond to the production and decay of the phosphorylated protein respectively. Both of these reactions are assumed to be defined by Michaelis-Menten kinetics. The maximum rate of phosphorylated protein production is a function of the cellular cytosolic Ca²⁺ concentration. As [Ca²⁺] oscillates, production of W* follows. The simulations of W* in this work use the same parameter values as Goldbeter, and can be found in the Supplementary Material. Figure 5 shows the time-dependent Ca²⁺ profile for a single EC and SMC taken from the center of the low ATP concentration in the lateral region, whilst Figure 6 shows the resulting phosphorylated protein W*using the SMC Ca²⁺ taken from the profile shown in Figure 5.

FIGURE 5

Figure 5. Time-dependent calcium profile for both SMC (blue) and EC (red) taken from the center of the lateral region. The SMC Ca²⁺ concentration spiked sporadicly, but for the majority of the time remained at a low concentration (approximately 0.2 μM). EC Ca²⁺ concentration oscillated within a small range approximately between 0.6 and 0.7 μM.

FIGURE 6

Figure 6. Time-dependent profile of the fraction of phosphorylated protein using the SMC Ca²⁺ profile shown in Figure 5. The phosphorylated protein dynamics follow the Ca²⁺ dynamics closely due to the relatively fast reaction rates.

The Ca²⁺ profile shows sporadic spiking with maximum amplitudes of up to 0.8 μM. However, most of the time the Ca²⁺ concentration is low: approximately 0.2 μM. Figure 6 shows that the phosphorylated protein dynamics follow the Ca²⁺ dynamics due to the relatively fast reaction rates, compared to the Ca²⁺ frequency, with maximum values of 0.95 μM. Most of time the phosphorylated protein concentration is low. On the basis of this result the lateral region, where one would expect plaque growth over time, the activation of MAPK pathway is low.

Since plaque formation is a chronic condition, and it may take months to form, a time-averaged value for the fraction of the phosphorylated protein is a more realistic parameter to investigate, which is defined as follows:

$\begin{array}{ll}\overline{W^{*}}=\frac{1}{T_{max}-T_{min}}{\displaystyle {\int }_{T_{min}}^{T_{max}}}{W}^{*}(s)ds,& (10)\end{array}$

where T_min is zero, and T_max to the time at the end of the simulation.

Figure 7 shows these average values as functions of cell type and region for the three angles of bifurcation. While in the apical region, both the time averaged EC and SMC $\overline{W^{*}}$values are increasing functions of bifurcation angle, that is not the case for the lateral region. Here EC-based $\overline{W^{*}}$increases with angle, but SMC-based $\overline{W^{*}}$ is a monotonic decreasing function and of much lower value. In the apical region the time-varying EC [Ca²⁺] values are near-constant and of relatively large value, of the order of 0.9 μM (data not shown), and the SMC [Ca²⁺] values are again near-constant, but lower, on average 0.4 μM (data not shown).

FIGURE 7

Figure 7. Protein phosphorylation for ECs and SMCs in geometries with varying angles. In the apical region, both the time averaged EC and SMC $\overline{W^{*}}$ values are increasing functions of bifurcation angle. Lateral seam ECs also follow this behavior, while lateral SMCs do not.

4. Discussion

All three simulations presented in this work showed propagating waves of Ca²⁺ in both the SMC and EC layers following the introduction of a luminal agonist, in this case ATP. The agonist binds to endothelial receptors, induces the production of IP₃, and consequently the release of Ca²⁺ in both ECs and SMCs. Both ECs and SMCs are excitable and their Ca²⁺ concentrations can oscillate under certain conditions (Kapela et al., 2008). The work of Koenigsberger et al. (2005) used a strip of coupled ECs and SMCs to show that cells could develop in-phase oscillations, otherwise known as vasomotion, however those simulations did not include spatial agonist variation, nor the resultant wave propagation. Since the cells in the experiments presented in this work did not produce in-phase [Ca²⁺] variations, the phenomenon produced is not vasomotion.

Complex behavior in SMC Ca²⁺ dynamics was observed at the lateral bifurcation seams. In particular, these areas exhibited relatively higher Ca²⁺ concentrations in ECs, while Ca²⁺ waves in SMCs propagated in a slower, non-uniform manner, and over significantly shorter distances. These sporadic Ca²⁺ oscillations in the SMC seem to be primarily due to the cell being in or close to a steady state and receiving a diffusive flux of Ca²⁺ from the EC via gap junctions, which perturbs them into an oscillatory state and, as a consequence, provides sufficient cytosolic Ca²⁺ to excite the cells into continuing the propagating wave. In agreement with previous simulations presented in the work of Zakkaroff et al. (2016), the Ca²⁺ waves propagate through regions that are otherwise devoid of oscillations.

Pinto and Campos (2016) suggested that areas under arterial bifurcations are susceptible to potential plaque development, due to observed low time-averaged wall shear stress (TAWSS) and high relative resident time (RRT) values. Dong et al. (2014) presented a more specific hypothesis—that it is these areas under bifurcations with wider angles that are at greater risk of plaque formation. We observed low time-averaged Ca²⁺ concentrations in SMCs, in regions that CFD experiments have identified as sites most likely to experience plaque development (Steinman et al., 2002; Peiffer et al., 2013; Pinto and Campos, 2016) but contrastingly not for ECs. Further, we note the low Ca²⁺ concentrations in these areas are more prominent in arterial geometries with wider bifurcation angles, in agreement with the findings of Dong et al. (2014).

4.1. NO Production and Atheroprone Conditions Considerations

There is considerable evidence that the production of NO is reduced in areas of low WSS (Kanai et al., 1995; Malek et al., 1999; Kavdia and Popel, 2003), which in turn can lead to endothelial dysfunction (Davignon and Ganz, 2004) and may result in an increased probability of atherogenesis. The eNOS-mediated production of NO is a function of cytosolic Ca²⁺ concentration. Work by Comerford et al. (2008) showed that areas of eNOS reduced as the angle of the bifurcation increased. However, in terms of Ca²⁺ oscillations nitric oxide has no influence on the oscillatory dynamics, since this is primarily driven by the relative fluxes of the CICR mechanism and the SERCA pump associated with the SR storage of Ca²⁺.

Although the NO pathway does not appear in the presented model, we can rely on the eNOS and NO data from previous work. For example, the work of Dormanns et al. (2016) utilized the same model for eNOS generation as described by Comerford et al. (2008); their data indicated that the EC and SMC concentrations of NO, as determined from the eNOS concentration, were similar in value. Thus, we can surmise that even for conditions of oscillating Ca²⁺ concentration, areas of low ATP and time-averaged cytosolic EC and SMC Ca²⁺ concentration would also produce low areas of NO concentration, thereby possibly enhancing the progression of EC dysfunction and subsequent atherogenesis.

On the basis of a coupled cellular model combined with fluid dynamics, the areas of low WSS, ATP, and low NO seem to correlate with areas susceptible to plaque formation. At the same time the simulations reported here indicate that time averaged EC [Ca²⁺] is considerably higher than that of SMCs in the lateral region, giving rise to NO concentrations, compatable with the apical region where, from experiment, the probablity of plaque growth is low. This conflicting result warrants further investigation, and hence the inclusion of the full NO pathway along with WSS-mediated NO production.

4.2. Myoendothelial Gap Junctions, Cx43, and WSS

It has been suggested that geometric and subsequent biochemical changes are related to atherosclerotic regions (Caro et al., 1971; Ku et al., 1985). Ca²⁺ ions move between ECs and SMCs in healthy cellular states, primarily through the Cx43 connexin protein subunit. A concentration change in one layer will therefore affect the other, and so Ca²⁺ dynamics in both ECs and SMCs should be considered when investigating the onset of vascular disease.

The review of cell coupling and atherosclerosis by Burnier et al. (2009) states that “lesions are typically observed in arterial vessels with high WSS such as bifurcations,” in contrast to most accepted theories that low WSS is associated with atherosclerosis. The review, however, also states that the Cx43 gap junction is thought to be sensitive to mechanic stimuli due to its overexpression during shear stress, while there is evidence suggesting that this particular gap junction has increased expression in the intima during the early phases of atherosclerosis. The work of Inai et al. (2004) provides some evidence that Cx43 is increased in areas of high WSS gradients, although this was only expressed in ECs of cardiac valves. At the same time the work of DePaola et al. (1999) showed that Cx43 expression increased with increasing WSS, as well as disorganization of the Cx43 protein.

4.3. Frequency Encoding and Ca²⁺ Dynamics Information Transport

Meyer and Stryer (1988) suggested that calcium-activated proteins, such as calmodulin, are highly suited for detecting spikes in ionic concentrations, where calcium could excite different effector systems. First Goldbeter and Koshland (1981), and later Goldbeter et al. (1990) proposed the idea that [Ca²⁺] oscillations might be encoded in terms of their frequency. They suggested that one possible scenario of information transport through frequency encoding may be based on protein phosphorylation, but only in the case of zero-order kinetics, where the substrate concentration is far in excess. This provides support to the premise that protein phosphorylation is highly sensitive to changes in the kinetics. Goldbeter's analysis indicates that the time variation of the phosphorylated protein is directly associated with the oscillations of cytosolic Ca²⁺ and to the subsequent variation in kinase activity, where a higher frequency of the SMC Ca²⁺ oscillation produces a higher mean value of the phosphorylated protein.

Research by Goldbeter and Koshland (1981) and Goldbeter et al. (1990) further state that when “the period of oscillations is longer, the protein undergoes significant dephosphorylation from one Ca²⁺ peak to the next”, resulting in larger values of dephosphorylated protein. The results presented here indicate that along the outer wall of the bifurcations in areas of low ATP, where cytosolic Ca²⁺ oscillations are uncoordinated, the encoding information is disrupted causing possible dysfunction to those protein pathways. For example, as noted before MAPK pathways have been implicated in foam cell formation, EC activation, and vascular smooth muscle migration and proliferation, and several receptor tyrosine kinases (RTKs) pathways in a variety of cells have been shown to be important for the development of atherosclerosis.

It emerges that no single marker exists for plaque formation and that a number of conditions have to be present in order to initiate atherosclerotic conditions. On the basis of Ca²⁺ dynamics alone, aligned with variation in geometry and luminal wall agonist concentration, we conclude that initiation of atherosclerotic plaque may not be due to frequency encoded Ca²⁺ phosphorylation of the protein alone. However, combinations of more complex cellular models and the introduction of other pathways such as NO, along with the use of numerical simulations could provide further insight into the relationship of vascular geometry and atherosclerotic plaque formation.

5. Conclusions

This research examined single-cell dynamics on a macro-scale, and the presented model was compared with other mathematical models of cellular physiology in the related literature. The underlying coupled cells model was extended to include a detailed IP₃ mediated Ca²⁺ pathway, and gap junction currents that replace simple linear gradients for membrane potential and Fickian diffusion of ions. The premise of this work was such that variations in the frequency and dynamics of the Ca²⁺ waves in areas of low WSS have the potential to influence the kinase reactions of proteins under the assumption that Ca²⁺ encodes for high sensitivity phosphorylation and possibly provides a marker for atherogenesis.

The macro-scale simulations in this work demonstrated oscillatory behavior in SMC and EC species' concentrations. Moreover, significantly lower Ca²⁺ concentrations were observed in areas under the lateral surface of bifurcations, in SMCs. The SMC Ca²⁺ dynamics in these areas were difficult to analyse at any given time point: waves propagated in a slow non-uniform manner, and over significantly shorter distances. Post-processing, in the form of time-averaging, was applied to gain insight into the emergent properties of these dynamics. This showed significantly lower Ca²⁺ concentrations over time at the lateral bifurcation seams compared to the other surface regions of the arterial segments. Ca²⁺ concentration frequencies and average values in the lateral domain for SMCs were relatively small. In contrast EC Ca²⁺ dynamics showed that average values were high in both apical and lateral regions. Since experiments have shown that plaque growth is rare in the apical region this would indicate that frequency encoding would not be a prominent mechanism in the very early stages of plaque growth.

In conjunction with findings reported in the literature, the results of this work suggest that arterial bifurcations containing wider angles may indeed influence the development of atherosclerosis, by means of disturbed flow and, consequently, lower SMC Ca²⁺ concentrations. The overal outcome of this work suggests that although Ca²⁺ dynamics are of considerable importance in cellular functions, other related pathways may combine in a non-linear fashion to produce an environment conducive to plaque formation.

This research provides a flexible simulation platform that allows for the investigation of cellular dynamics with geometrically varying, multi-scale arterial sections. The potential influence of bifurcation angulation on atherogenesis has been studied, and the investigation of many other phenomena is now possible. Finally, all additions and modifications made to this platform were designed and implemented to promote scalability and further model development. Future simulations of arterial physiology, including further detailed cellular chemistry and surface-meshes that span the entire coronary tree, should benefit from these extensions.

6. Limitations and Future Work

It is important to note that biophysical values and coefficients used in the presented model were collected from different configurations of experiments. For example, the coefficients of electro-diffusive coupling can vary widely depending on experimental setup. These considerations present a potential limitation for in-silico experiments; however, we hope that future research in physiology may help improve the homogenization of experimental findings for computational modeling.

For a given mesh, the ECs and SMCs on each quadrilateral were aligned orthogonally with respect to each other, where ECs align axially with the flow and the SMCs align circumferentially. These arrangements were generated such that they communicate with four neighboring cells homocellularly, when in fact an EC and SMC can typically communicate with up to six neighboring ECs and SMCs, respectively (Koenigsberger et al., 2005). Additionally, ECs tend to align in the direction of the dominant local WSS vector component whereas in the present case ECs are always aligned on an axial direction. The random nature of alignment for ECs at positions of low WSS has the potential to affect the homotypic and, possibly, heterotypic coupling coefficients as shown by Chiu and Chien (2011).

The NO pathway, including cyclic guanosine monophosphate (cGMP) in the EC, was not included in this model. This pathway describes the production of NO as a result of agonistic molecules within the lumen binding to EC receptors. NO has known atheroprotective properties, such as inhibiting the production of adhesion molecules that encourage the recruitment of monocytes. Following the implementation of this pathway in macro-scale simulations, it will be possible to examine NO dynamics across large sections of ECs, and, possibly, provide further insight into the process of atherosclerotic plaque development. Similarly, the introduction of stretch-activated channels and vasculature motion in the model may improve the fidelity of future simulation results.

Many geometric variations were considered in a study conducted by Lee et al. (2008), among which was bifurcation angulation. This property was chosen for analysis as a result of findings that bifurcations with wider angles are more likely to experience greater flow detachment (Dong et al., 2014). However, there are many other geometric properties (such as bifurcation tortuosity, for example) that are worthy of future exploration.

Shaikh et al. (2012) introduced multiple coupling cases to earlier Ca²⁺ dynamics simulations in order to model both healthy and pathological conditions in coupled cells behavior. These coupling cases were not the focus of the research presented here—the particular coupling cases modeling unhealthy cellular behavior already assume mid-to-late stages of atheroma development. However, future work should include the creation of new coupling cases that reflect initial stages of disease, before plaque formation is observable. Further, these coupling cases should vary both temporally and spatially as noted above due to the variation in EC alignment. Similarly, in their discussion of Ca²⁺ dynamics in the vascular endothelium, Taylor and Francis (2014) stated that it was important for future studies to move away from assumptions based on global Ca²⁺ changes and instead focus on spatially and temporally relevant aspects of real-time signaling.

The results of our simulations show abnormal cellular dynamics at specific locations on arterial geometries; therefore different coupling parameters should be employed at these locations to reflect potentially emergent dysfunctional cellular communication, but there seems to be little evidence to indicate the coupling parameters variation with respect to WSS. This is particularly important, given that simulations presented here span physiological time periods of minutes, while atherosclerotic plaques develop over time scales of months. Modification of localized communication channels as a function of atypical cellular dynamics, such as significantly lower Ca²⁺ concentrations, would allow us to simulate the progression of vascular disease in relatively short simulations.

Data Availability Statement

The computational model and required input files for this study can be found in the Coupled Cells Github release: https://github.com/BlueFern/CoupledCells/releases/tag/v-4.0.

Author Contributions

All authors contributed equally to the direction and ideas of the presented work. SD and CZ generated the new surface meshes, extended the physiological model, and ran parallel simulations. SM produced the input ATP maps using computational fluid dynamics methods. SD, CZ, and TD developed the manuscript.

Funding

The University of Canterbury provided a UC HPC Master's scholarship that partially funded this work.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank the members of the University of Canterbury High Performance Computing Center, Tony Dale, and François Bissey for their enthusiastic support and expert advice on the technical front. We also thank Allanah Kenny and Michael Plank for proofreading and providing valuable suggestions to the manuscript. The authors wish to acknowledge the contribution of New Zealand eScience Infrastructure (NeSI) to the results of this research. New Zealand's national compute and analytics services and team are supported by NeSI and funded jointly by NeSI's collaborator institutions and through the Ministry of Business, Innovation and Employment.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01528/full#supplementary-material

Footnotes

1. ^https://goo.gl/AQivPz—University of Canterbury High Performance Computing channel on YouTube.

References

Ayachit, U. (2015). The ParaView Guide: A Parallel Visualization Application. New York, NY: Kitware, Inc.

Bennett, M. R., Farnell, L., and Gibson, W. (2005). A quantitative model of purinergic junctional transmission of calcium waves in astrocyte networks. Biophys. J. 89, 2235–2250. doi: 10.1529/biophysj.105.062968

PubMed Abstract | CrossRef Full Text | Google Scholar

Boileau, E., Parthimos, D., and Nithiarasu, P. (2015). An extended computational framework to study arterial vasomotion and its links to vascular disease. Biomed. Technol. 11, 129–144. doi: 10.1007/978-3-319-10981-7_9

CrossRef Full Text | Google Scholar

Burnier, L., Fontana, P., Angelillo-Scherrer, A., and Kwak, B. R. (2009). Intercellular communication in atherosclerosis. Physiology 24, 36–44. doi: 10.1152/physiol.00036.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Caro, C. G., Fitz-Gerald, J. M., and Schroter, R. C. (1971). Atheroma and arterial wall shear observation, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis. Proc. R. Soc. Lond. B Biol. Sci. 177, 109–133.

PubMed Abstract | Google Scholar

Chiu, J. J., and Chien, S. (2011). Effects of dsisturbed flow on vascular endothelium: pathophysiological basis and clinical perspectives. Physiol. Rev. 91, 327–387. doi: 10.1152/physrev.00047.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Comerford, A., Plank, M. J., and David, T. (2008). Endothelial nitric oxide synthase and calcium production in arterial geometries: an integrated fluid mechanics/cell model. J. Biomech. Eng. 130:011010. doi: 10.1115/1.2838026

PubMed Abstract | CrossRef Full Text | Google Scholar

David, T. (2003). Wall shear stress modulation of ATP/ADP concentration at the endothelium. Ann. Biomed. Eng. 31, 1231–1237. doi: 10.1114/1.1615574

PubMed Abstract | CrossRef Full Text | Google Scholar

Davignon, J., and Ganz, P. (2004). Role of endothelial dysfunction in atherosclerosis. Circulation 109, 27–32. doi: 10.1161/01.CIR.0000131515.03336.f8

PubMed Abstract | CrossRef Full Text | Google Scholar

DePaola, N., Davies, P. F., Pritchard, W. F., Florez, L., Harbeck, N., and Polacek, D. C. (1999). Spatial and temporal regulation of gap junction connexin43 in vascular endothelial cells exposed to controlled disturbed flows in vitro. Proc. Natl. Acad. Sci. U.S.A. 96, 3154–3159. doi: 10.1073/pnas.96.6.3154

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, J., Sun, Z., Inthavong, K., and Tu, J. (2014). Fluid-structure interaction analysis of the left coronary artery with variable angulation. Comput. Methods Biomech. Biomed. Eng. 18, 1500–1508. doi: 10.1080/10255842.2014.921682

PubMed Abstract | CrossRef Full Text | Google Scholar

Dormanns, K., Brown, R. G., and David, T. (2016). The role of nitric oxide in neurovascular coupling. J. Theor. Biol. 394, 1–17. doi: 10.1016/j.jtbi.2016.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Elliott, K. J., and Eguchi, S. (2015). “Phosphorylation regulation by kinases and phosphatases in atherosclerosis,” in Atherosclerosis: Risks, Mechanisms, and Therapies, eds H. Wang and C. Patterson (Hoboken, NJ: John Wiley & Sons Inc.), 453–465.

Google Scholar

Goldbeter, A., Dupont, G., and Berridge, M. J. (1990). Minimal model for signal-induced Ca²⁺ oscillations and for their frequency encoding through protein phosphorylation. Proc. Natl. Acad. Sci. U.S.A. 87, 1461–1465.

PubMed Abstract | Google Scholar

Goldbeter, A., and Koshland, D. E. (1981). An amplified sensitivity arising from covalent modification in biological systems. Proc. Natl. Acad. Sci. U.S.A. 78, 6840–6844.

PubMed Abstract | Google Scholar

Gorman, M. W., Feigl, E. O., and Buffington, C. W. (2007). Human plasma atp concentration. Clin. Chem. 53, 318–325. doi: 10.1373/clinchem.2006.076364

PubMed Abstract | CrossRef Full Text | Google Scholar

Gropp, W., Hoefler, T., Thakur, R., and Lusk, E. (2014). Using Advanced MPI: Modern Features of the Message-Passing Interface. Cambridge, MA: MIT Press.

Google Scholar

Hadfield, J., Plank, M. J., and David, T. (2013). Modeling secondary messenger pathways in neurovascular coupling. Bull. Math. Biol. 75, 428–443. doi: 10.1007/s11538-013-9813-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hakala, J. K., Lindstedt, K. A., Kovanen, P. T., and Pentikainen, M. O. (2006). Low-density lipoprotein modified by macrophage-derived lysosomal hydrolases induces expression and secretion of il-8 via p38 mapk and nf-κb by human monocyte-derived macrophages. Arterioscler. Thromb. Vasc. Biol. 26, 2504–2509. doi: 10.1161/01.ATV.0000245796.97133.ad

PubMed Abstract | CrossRef Full Text | Google Scholar

Heo, K. S., Fujiwara, K., and Abe, J. (2014). Shear stress and atherosclerosis. Mol. Cells 37, 435–440. doi: 10.1016/S0021-9290(06)84109-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban, R., Shumaker, D. E., et al. (2005). SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Trans. Math. Softw. 31, 363–396. doi: 10.1145/1089014.1089020

CrossRef Full Text | Google Scholar

Inai, T., Mancuso, M. R., McDonald, D. M., Kobayashi, J., Nakamura, K., and Shibata, Y. (2004). Shear stress-induced upregulation of connexin 43 expression in endothelial cells on upstream surfaces of rat cardiac valves. Histochem. Cell Biol. 122, 477–483. doi: 10.1007/s00418-004-0717-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobsen, J. C., Aalkjaer, C., Matchkov, V. V., Nilsson, H., Freiberg, J. J., and Holstein-Rathlou, N. H. (2008). Heterogeneity and weak coupling may explain the synchronization characteristics of cells in the arterial wall. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 366, 3483–3502. doi: 10.1098/rsta.2008.0105

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobsen, J. C., Aalkjaer, C., Nilsson, H., Matchkov, V. V., Freiberg, J., and Holstein-Rathlou, N. H. (2007). A model of smooth muscle cell synchronization in the arterial wall. Am. J. Physiol. Heart Circ. Physiol. 293, H229–H237. doi: 10.1152/ajpheart.00727.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Johny, J. P., Plank, M. J., and David, T. (2017). Importance of altered levels of SERCA, IP3R, and RyR in vascular smooth muscle cell. Biophys. J. 112, 265–287. doi: 10.1016/j.bpj.2016.11.3206

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanai, A. J., Strauss, H. C., Truskey, G. A., Crews, A. L., Grunfeld, S., and Malinski, T. (1995). Shear stress induces atp-independent transient nitric oxide release from vascular endothelial cells, measured directly with a porphyrinic microsensor. Circ. Res. 77, 284–293.

PubMed Abstract | Google Scholar

Kapela, A., Bezerianos, A., and Tsoukias, N. M. (2008). A mathematical model of Ca2+ dynamics in rat mesenteric smooth muscle cell: agonist and NO stimulation. J. Theor. Biol. 253, 238–260. doi: 10.1016/j.jtbi.2008.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Kapela, A., Nagaraja, S., and Tsoukias, N. M. (2009). A mathematical model of vasoreactivity in rat mesenteric arterioles. ii. conducted vasoreactivity. Am. J. Physiol. Heart Circ. Physiol. 298, H52–H65. doi: 10.1152/ajpheart.00546.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Kavdia, M., and Popel, A. S. (2003). Wall shear stress differentially affects no level in arterioles for volume expanders and hb-based o₂ carriers. Microvasc. Res. 66, 49–58. doi: 10.1016/S0026-2862(03)00008-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Koenigsberger, M., Sauser, R., Bény, J. L., and Meister, J. J. (2005). Role of the endothelium on arterial vasomotion. Biophys. J. 88, 3845–3854. doi: 10.1529/biophysj.104.054965

PubMed Abstract | CrossRef Full Text | Google Scholar

Koenigsberger, M., Sauser, R., Bény, J. L., and Meister, J. J. (2006). Effects of arterial wall stress on vasomotion. Biophys. J. 91, 1663–1674. doi: 10.1529/biophysj.106.083311

PubMed Abstract | CrossRef Full Text | Google Scholar

Koenigsberger, M., Seppey, D., Bény, J. L., and Meister, J. J. (2010). Mechanisms of propagation of intercellular calcium waves in arterial smooth muscle cells. Biophy. J. 99, 333–343. doi: 10.1016/j.bpj.2010.04.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Ku, D. N., Giddens, D. P., Zarins, C. K., and Glagov, S. (1985). Pulsatile flow and atherosclerosis in the human carotid bifurcation. Positive correlation between plaque location and low oscillating shear stress. Arteriosclerosis 5, 293–302.

PubMed Abstract | Google Scholar

Lee, S. W., Antiga, L., Spence, J. D., and Steinman, D. A. (2008). Geometry of the carotid bifurcation predicts its exposure to disturbed flow. Stroke 39, 2341–2347. doi: 10.1161/STROKEAHA.107.510644

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemon, G., Gibson, W. G., and Bennett, M. R. (2003). Metabotropic receptor activation, desensitization and sequestration - I: modelling calcium and inositol 1,4,5-trisphosphate dynamics following receptor activation. J. Theor. Biol. 223, 93–111. doi: 10.1016/S0022-5193(03)00079-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Libby, P. (2002). Inflammation in atherosclerosis. Nature 420, 868–874. doi: 10.1038/nature01323

PubMed Abstract | CrossRef Full Text | Google Scholar

Lidington, D., Ouellette, Y., and Tyml, K. (2000). Endotoxin increases intercellular resistance in microvascular endothelial cells by a tyrosine kinase pathway. J. Cell. Physiol. 185, 117–125. doi: 10.1002/1097-4652(200010)185:1<117::AID-JCP11>3.0.CO;2-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Malek, A. M., Alper, S. L., and Izumo, S. (1999). Hemodynamic shear stress and its role in atherosclerosis. J. Am. Med. Assoc. 282, 2035–2042.

PubMed Abstract | Google Scholar

Meyer, T., and Stryer, L. (1988). Molecular model for receptor-stimulated calcium spiking. Proc. Natl. Acad. Sci. U.S.A. 85, 5051–5055.

PubMed Abstract | Google Scholar

Murray, C. D. (1926). The physiological principle of minimum work: a reply. J. Gen. Physiol. 12:445.

Google Scholar

Parthimos, D., Edwards, D. H., and Griffith, T. (1999). Minimal model of arterial chaos generated by coupled intracellular and membrane ca2+ oscillators. Am. J. Physiol. Heart Circ. Physiol. 277, H1119–H1144.

PubMed Abstract | Google Scholar

Peiffer, V., Sherwin, S. J., and Weinberg, P. D. (2013). Does low and oscillatory wall shear stress correlate spatially with early atherosclerosis? A systematic review. Cardiovasc. Res. 99, 242–250. doi: 10.1093/cvr/cvt044

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinto, I. S., and Campos, J. B. (2016). Numerical study of wall shear stress-based descriptors in the human left coronary artery. Comput. Methods Biomech. Biomed. Eng. 19, 1443–1455. doi: 10.1080/10255842.2016.1149575

PubMed Abstract | CrossRef Full Text | Google Scholar

Plank, M. J., Wall, D. J. N., and David, T. (2006). Atherosclerosis and calcium signalling in endothelial cells. Prog. Biophys. Mol. Biol. 91, 287–313. doi: 10.1016/j.pbiomolbio.2005.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Schuster, A., Bény, J. L., and Meister, J. J. (2003). Modelling the electrophysiological endothelial cell response to bradykinin. Eur. Biophys. J. 32, 370–380. doi: 10.1007/s00249-003-0279-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaikh, M. A., Wall, D. J., and David, T. (2012). Macro-scale phenomena of arterial coupled cells: a massively parallel simulation. J. R. Soc. Interface 9, 972–987. doi: 10.1098/rsif.2011.0453

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinman, D. A., Thomas, J. B., Ladak, H. M., Milner, J. S., Rutt, B. K., and David Spence, J. (2002). Reconstruction of carotid bifurcation hemodynamics and wall thickness using computational fluid dynamics and MRI. Magnet. Reson. Med. 47, 149–159. doi: 10.1002/mrm.10025

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, M. S., and Francis, M. (2014). Decoding dynamic Ca2+ signaling in the vascular endothelium. Front. Physiol. 5:447. doi: 10.3389/fphys.2014.00447

CrossRef Full Text | Google Scholar

Thomas, J. B., Antiga, L., Che, S. L., Milner, J. S., Steinman, D. A. H., Spence, J. D., et al. (2005). Variation in the carotid bifurcation geometry of young versus older adults: Implications for geometric risk of atherosclerosis. Stroke 36, 2450–2456. doi: 10.1161/01.STR.0000185679.62634.0a

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Rijen, H., Van Kempen, M. J., Analbers, L. J., Rook, M. B., Van Ginneken, A. C., Gros, D., et al. (1997). Gap junctions in human umbilical cord endothelial cells contain multiple connexins. Am. J. Physiol. Cell Physiol. 272, C117–C130.

PubMed Abstract | Google Scholar

Wan, J., Ristenpart, W. D., and Stone, H. A. (2008). Dynamics of shear-induced atp release from red blood cells. Proc. Natl. Acad. Sci. U.S.A. 105, 16432–16437. doi: 10.1073/pnas.0805779105

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamamoto, Y., Klemm, M. F., Edwards, F. R., and Suzuki, H. (2001). Intercellular electrical communication among smooth muscle and endothelial cells in guinea-pig mesenteric arterioles. J. Physiol. 535, 181–195. doi: 10.1111/j.1469-7793.2001.00181.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Clark, J. W., Bryan, R. M., and Robertson, C. S. (2005). Mathematical modeling of the nitric oxide/cgmp pathway in the vascular smooth muscle cell. Am. J. Physiol. Heart Circ. Physiol. 289, H886–H897. doi: 10.1152/ajpheart.00216.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Zakkaroff, C., Moore, S., Dowding, S., and David, T. (2016). 3D time-varying simulations of Ca(2+) dynamics in arterial coupled cells: a massively parallel implementation. Int. J. Numer. Methods Biomed. Eng. 33:e02786. doi: 10.1002/cnm.2786

CrossRef Full Text | Google Scholar