3D time‐varying simulations of Ca2+ dynamics in arterial coupled cells: A massively parallel implementation

Summary Preferential locations of atherosclerotic plaque are strongly associated with the areas of low wall shear stress and disturbed haemodynamic characteristics such as flow detachment, flow recirculation and oscillatory flow. The areas of low wall shear stress are also associated with the reduced production of adenosine triphosphate in the endothelial layer, as well as the resulting reduced production of inositol trisphosphate (IP3). The subsequent variation in Ca2+ signalling and nitric oxide synthesis could lead to the impairment of the atheroprotective function played by nitric oxide. In previous studies, it has been suggested that the reduced IP3 and Ca2+ signalling can explain the correlation of atherosclerosis with induced low WSS and disturbed flow characteristics. The massively parallel implementation described in this article provides insight into the dynamics of coupled smooth muscle cells and endothelial cells mapped onto the surface of an idealised arterial bifurcation. We show that variations in coupling parameters, which model normal and pathological conditions, provide vastly different smooth muscle cell Ca2+ dynamics and wave propagation profiles. The extensibility of the coupled cells model and scalability of the implementation provide a solid framework for in silico investigations of the interaction between complex cellular chemistry and the macro‐scale processes determined by fluid dynamics. © 2016 The Authors. International Journal for Numerical Methods in Biomedical Engineering published by John Wiley & Sons Ltd.


INTRODUCTION
The recently coined term in silico research refers to computer simulations of complex biological systems dynamics. This research approach offers the potential of increasing the speed of knowledge discovery [1,2]. In silico experiments have the potential to provide insight into the observations obtained by the experimental science. In the context of biological systems dynamics, the in silico simulations enable the rapid pruning of the parameter search space for the refinement and integration of cellular-level models into biologically realistic macro-scale models. The large-scale physiological simulations described here were designed to provide insight into the effects of the luminal concentration variations on adenosine triphosphate (ATP)-dependent dynamics in the coupled endothelial cells (ECs) and smooth muscle cells (SMCs) making up an arterial wall. The simulations of this nature provide a unique opportunity to perform experiments, which would never be possible in the in vivo and in vitro settings. For example, various specific pathological conditions, as described further in the text, can be simulated by changing the homocellular and heterocellular coupling parameters. Numerical simulations of this nature have never been attempted before at the scale of millions of coupled cells.
The presented work is an initial investigation into the relationship between micro-scale cellular dynamics and the resulting 'emergent behaviour' that occurs on a macro-scale exemplified by the size of atherosclerotic plaques. There is a clear correlation between the dynamics and the increase in intimal hyperplasia as found in in vitro and in vivo (cadaver) studies. Indeed, there is considerable evidence in the work of [12] and [21] demonstrating that eNOS and the production of NO is reduced in areas of low WSS, which in turn leads to endothelial dysfunction [22], resulting in the increased probability of atherosclerotic plaque formation. The eNOS-mediated production is a function of cytosolic calcium, which is available from the stores through the IP 3 -mediated channel. IP 3 concentration is determined by the P2Y receptor mediated by agonists such as ATP. Finally, we note that as early as 1988 [23], there was found a clear relationship between WSS/calcium-induced potassium currents; these biochemical changes were hypothesised to be related to atherosclerotic regions [11,24]. In our work, we explore the relationship of calcium dynamics and arterial geometry as a way of determining how spatially and time-varying agonist concentration can produce calcium transients over much larger scales than a single cell, giving credence to the idea that plaques grow both upstream and downstream. This work offers the possibility of understanding further the process of the initiation and growth of atherosclerotic plaques.

Resistance arterioles versus coronary artery simulations
The majority of experiments showing electrical signalling are conducted using small diameter vessels either denuded of endothelium or vessels that contain only ECs ( [25]), pressurised vessels without fluid flow ( [26]) or that the current through the gap junction remains constant ( [25]). These experiments also do not show oscillatory behaviour as a function of space because this is difficult if not impossible to replicate under laboratory conditions. These constraints do not occur in our simulations, thus allowing for a more physiologically correct simulation.
The bi-directional relationship between Ca 2C and WSS is a lot more prominent in small diameter resistance arteries. Vasomotion is not considered in this work because this phenomenon occurs more readily in resistance vessels rather than, say, the coronary tree. The focus of the paper is on much larger atherosclerosis prone vessels such as the coronary where contractility does not vary as much as that seen in smaller resistance arterioles. Fluid dynamics calculations show that WSS varies little in terms of the dilation or contraction of the artery (unless it is a small resistance vessel) compared with the base WSS induced by the vascular geometry of a bifurcation.

Parallel coupled cells simulations
The macro-scale simulations reported in this work were implemented on the Blue Gene parallel architecture where a number of coupled EC/SMC units are grouped into a single task/domain in the parallel computational environment. These tasks interact with each other by exchanging state variable values for the neighbouring cells along the edge of the computational domain at every time step during the simulation process. The interaction between the groups of coupled EC/SMC units was implemented through the MPI. Full details of the computational solution are provided in Section 2.3.
In this article, Section 2 briefly reviews the model for a single pair of coupled ECs/SMCs, which is based on the work of [18]. The mass transfer dynamics are modelled by systems of non-linear ordinary differential equations (ODEs). The implementation of this coupled-cell model along with initial simulations has been previously described in the work of [27]. The core of the computation model described here is a single-coupled EC/SMC unit modelled by a system of nine ODEs with four and five equations for an EC and SMC, respectively.
The simulations reported in this work are divided into two groups. The first group is based on a synthetic ATP agonist map generated with a sigmoid function parametrised with the distance to the origin of the vessel segment, while the second group is based on the CFD-based agonist map generated from a model of physiologically-realistic mass transfer simulation of blood flow in a vessel. The results for both simulation groups are provided in Section 3, followed by discussion and conclusions in Section 4 and Section 5, respectively.

METHODS
The experimental methods described here extend the work reported previously by [27]. An overview of the coupled interactions between SMCs and ECs is provided in Section 2.1. The coupled cells physiological model remains unchanged; a brief review of the model is provided in Appendix A.1 for completeness. The cell coupling parameter variations in which were used in the simulations reported here are provided in Section 2.2. The computational solution, however, has been significantly extended to allow the modelling of arterial segments with complex geometric features such as branching and bifurcations. The critical difference from the computational solution presented by [27] is the 2D Cartesian communication between the groups of coupled cells comprising the discrete domains of the vessel surface. This Cartesian model of inter-domain communication supports coupled cells simulations with physiologically realistic arterial segments; the enhanced arterial surface generation is based on a method for parametric modelling of bifurcating structures. The computational solution description is provided in Section 2.3, including the Cartesian surface modelling and macro-scale inter-cellular communication. Section 2.4 describes the generation of synthetic ATP agonist maps, while Section 2.5 describes the CFD-based modelling of the ATP simulation inputs. Figure 1. Schematic representation of mass transfer dynamics in a single endothelial cell (EC)/smooth muscle cell (SMC) unit; this model accounts for the essential mechanisms of IP 3 -induced cytsolic Ca 2C release and the cascade of events following it. The diagram is adapted from the previous work by [27]. The numbers refer to the pathways described in the text.

Overview of a coupled cells unit
The adhesion of agonist (ATP) to the endothelial surface starts a complex cascade of reactions in the ECs and SMCs as shown in Figure 1. The following is a brief description of participating pathways: (1) Agonist binds to the purinergic (P2Y) receptors on the EC surface, activating the GPCR which then activates membrane-bound phospholipase C (PLC). PLC activation allows phosphorylation of PIP2. (2) PIP2 gives rise to IP 3 that is then released in the intracellular space. This nascent IP 3 binds to IP 3 receptor (IP 3 ) on the ER/SR surface. (3) IP 3 -bound IP 3 induces release of Ca 2C ions from the ER/SR into the cytosol. (4) The Ca 2C release from intracellular store sensitises the IP 3 further, which releases more Ca 2C , thus making a Ca 2C -rich domain in the cytosol in both EC and SMC. The excess of intracellular Ca 2C depolarises the membrane potential. (5) ER/SR has low-affinity binding sites on the cytosolic side of a channel, which constitutes a pump called the (SERCA) pump. Cytosolic Ca 2C encourages the replenishment of the intracellular stores via this pathway. (6) Ca 2C leaks from ER/SR consistently under a concentration gradient between cytosolic and ER/SR luminal Ca 2C and keeps the Ca 2C in equilibrium during a non-stimulated state of the cell. (7) In an EC, the cytosolic Ca 2C favours the influx of extracellular Ca 2C from non-selective cation channels. (8) The SERCA pump pushes out cytosolic Ca 2C to extracellular space. (9,10) In ECs, activation of KCa, upon binding to Ca 2C ions intracellularly at BKCa and SKCa, let K C ions move out of the cytosol. This hyperpolarises the membrane potential. (11) Although K C ions efflux is the main repolarising current, residual current (mainly consisting of monovalent ions) also contributes to membrane potential repolarisation. (12) The IP 3 concentration increases in SMC cytosol via transmission of IP 3 from coupled EC.
This IP 3 increase activates the downstream IP 3 -induced Ca 2C release. (13) The IP 3 induced, and calcium-induced calcium release (CICR) Ca 2C depolarises the membrane potential. (14) The membrane depolarisation results in the influx of Ca 2C from VOCCs, which will close upon repolarisation in the following steps. The ATP to J PLC conversion models reported by [28] and [29] show that the relationship between ATP and IP 3 is monotonic and increasing; thus it can be approximated by a linear function over the range of ATP. Therefore, a linear function is utilised to simulate the resulting relationship between ATP and the flux of IP 3 into the EC, that is J PLC . The range of ATP concentration is represented by the equivalent values of the J PLC ; the flux of IP 3 into the EC is determined by the concentration of ATP in the vessel lumen. The whole range of J PLC consists of three bands, which correspond to the ATP-dependent Ca 2C dynamics in an SMC. As shown in Figure 2, low and high bands of J PLC result in steady state Ca 2C concentrations, while the middle band shows the presence of Ca 2C oscillations. The oscillation envelope for this band indicates the relative high and low amplitudes for lower and higher J PLC concentrations.

Cell coupling variations
The variations of coupling parameters described here model the connexin gap junction permeability used in the simulations in this work. The Cx37 and Cx40 down-regulation along with Cx43 upregulation in early atheroma has been reported in the work of [30]. [31] have shown that Cx43 gap junctions have a conductance change of 50 nS for each mV change in membrane potential; this value is used for modelling heterocellular electrical coupling. Homocellular coupling parameter values were obtained from the work of [32] to set electrical resistance at 33 M (equivalent to 30 nS) as chosen in the work by [18]. With a membrane capacitance between ECs of 30 pF whilst for SMCs this is 10 pF and using the results of [33]. The macroscopic gap junctional resistance of 90 M (equivalent to 11 nS) [33] is the value of the homocellular coupling for both EC and SMC set as 1000 s 1 . The Ca 2C and IP 3 homocellular coupling parameters are treated as free in a similar way to the work of [18], assuming they are relatively weak; the value of is set to 0.05 s 1 . The coupling of membrane potential is a combination of simple ion diffusion across the gap junction and drift given by the product of ion concentration and spatial gradient of the membrane potential. For this particular model, the electro-diffusion drift is ignored as it has been shown to have negligible effect. Unpublished work with coupled cellular oscillatory models shows this to be the case. Table I summarises the coupling variations that have been used in the current work; the coupling parameter sets reproduce the parameters used in [27] in order to produce a comparison with the original 1D model with a 2D surface simulation. Cases 1 and 2 reflect 'healthy' (non-pathological) coupling cases, where Case 2 has the added heterocellular Ca 2C coupling. The remaining two coupling configurations are designed to simulate pathological cases. In Case 3, IP 3 coupling is added to the default homocellular EC coupling to simulate upregulation of Cx43 in lesion-prone areas described in the work of [30]. Case 4 has the homocellular EC and heterocellular membrane potential Ca 2C coupling disabled to simulate the coupling in progressive atherosclerosis when Cx43 connexin is the dominant mechanism of intercellular communication because of the down regulation of Cx37 and Cx40 connexins [5]. It is recognised that these four cases represent simple approximations to a complex physiological phenomena and should not be seen as a complete description of the mechanisms involved. However, they are put forward an initial framework from which further more complex pathways can be simulated. In particular, the role of IP 3 diffusion will need careful consideration due in part to its assumed smaller connexin channel permeability. Figure 2 shows the bifurcation diagrams (maximum/minimum Ca 2C concentration vs the IP 3 flux parameter J PLC ); for case 1, the range of IP 3 where oscillation exists is substantial whereas case 4 has a reduced domain, and the lower bifurcation points are shifted from approximately 3.7 to 3.0. Cases 2 and 3 are identical in their bifurcation points and maximum/minimum values of Ca 2C oscillations, because of the identical heterocellular coupling parameters. The bifurcation diagrams were generated for a single EC/SMC unit without homocellular coupling.

Computational solution methodology
Simulations of cellular dynamics even for small vessel segments require a substantial amount of computational resources because of the density and number of coupled cells comprising the endothelial and smooth muscle surfaces making up the arterial wall. The primary requirement in the development of the computation solution for the large-scale simulations was the decomposition of the coupled EC and SMC layers into a number of ODEs systems, each corresponding to a number of 'neighbourhoods' of coupled cells grouped into quadrilateral surface elements later mapped to distinct computational domains. With a suitably chosen granularity of the surface-to-domain decomposition, the entire bifurcation surface can be mapped onto individual computational nodes in a parallel architecture.

Bifurcation meshes.
A method for generic parameterisation of bifurcating structures published by [34] offers an elegant solution for generating high-quality surface meshes from centre lines of segmented arteries. For the presented cases, a simple symmetric bifurcation is used and is defined by three centre lines forming a Y-shaped bifurcation. For the present model, we do not utilise a patient-specific geometry but assume that the topological equivalent of an arterial bifurcation provides a sufficient surface. Each Y bifurcation can be decomposed into three semi-tubular surface segments by using the solution to a biharmonic equation, originally developed as a computer-aided design tool by Bloor and Wilson [35]. The biharmonic equation has the operator form given by where is the space vector .x; y;´/ and f .u; v/ can be either zero or some a priori defined function.
For the present, we use f .u; v/ D 0, where u and v are orthogonal coordinates, which provide axial and circumferential directions. The Neumann and Dirichlet boundary conditions for Eq. (2.1) define the rate at which the surface moves away from the boundary and its direction, respectively. Eq. (2.1) is discretised into m (u direction) and n (v direction) points on the boundary. With this method, the surface of each segment is unwrapped into a 2D m n grid in the u and v directions, respectively; in this case, m and n represent the longitudinal and circumferential number of quadrilaterals in each semi-tubular segment accordingly. The solution now becomes one of a 2D solution, rather than a 3D solution. To obtain a complete 3D surface, the semi-tubular segments can be joined together at the boundaries.  full bifurcation geometry. Each segment surface here is defined by the solution of the biharmonic equation. Figure 3(b) shows the Dirichlet and Neumann boundary conditions used as input for a single biharmonic solution.
The mapping of bifurcations to 2D-quadrilateral meshes enables the decomposition of the whole computational domain into local domains, where each domain corresponds to a single computational task in a parallel environment. Furthermore, this type of decomposition enables the use of MPI-specific virtual Cartesian topologies for efficient inter-processor communication. Each tubular segment in the computational solution presented here is mapped to a 2D Cartesian grid of size m n with a periodic boundary condition imposed along the longitudinal edge. In this case, the periodic boundary condition refers to the MPI-enabled communication-specific mapping of MPI tasks. It must be noted that the bifurcation surfaces are not constrained to symmetric or circular geometries as the solution of the biharmonic equation can be defined for any (continuous) boundary in .x; y;´/ space.  Figure 4(a) shows a bifurcation surface mesh consisting of Q D 4080 quadrilateral elements, where each element represents a cell neighbourhood on the bifurcation surface. With this type of surface decomposition, each quadrilateral domain corresponds to a parallel computational task running on a single core. The SMCs and ECs within each domain can then be placed into the quadrilaterals with little complexity because the .u; v/ directions provide physiologically correct axial and circumferential alignments for ECs and SMCs, respectively. The SMC and EC cells are aligned orthogonally with respect to each other where ECs align axially with the flow, whilst SMCs align circumferentially. Figure 4(b) shows ECs and SMCs arranged in two layers within a single quadrilateral on the bifurcation surface. With the approximate physical dimensions of SMCs and ECs, 50 5 m and 65 10 m, respectively, the SMCs and ECs are arranged into 'fundamental' surface units, where the area covered by five ECs corresponds to the area of 13 SMCs; this arrangement corresponds to the layer arrangement used by [27]. The surface units are placed into the quadrilaterals  in a 4 4 grid giving 288 cells per a single quadrilateral/domain with 208 and 80 SMCs and ECs, respectively. An implementation of the biharmonic equation solver in Fortran is freely available through online resources ‡ . The Fortran code was integrated into a C++ 'wrapper' as a part of a small library for generating surfaces from vessel centre lines; the library is using the relevant functionality of the (VTK) library [36] such as input, conversion of point clouds returned by the dbihar solver into surface meshes and output. The wrapper was implemented as a VTK-style filter to be used as a part of a pipeline for generating meshes from vessel centrelines.
Two synthetic arterial bifurcation surface meshes were generated for the simulations reported in this work. The geometric configuration properties, such as proportional length and angle of the branches, were similar in both cases; however, the second mesh was generated for the vessel surface approximately twice the size of the first mesh. The physical dimensions of the cells are only relevant when generating the surface from a centreline with a specified vessel diameter; the knowledge of the physical dimensions of cells enables the generation of meshes with the correct cell density for a given centreline and vessel diameter. For example, a surface mesh for a bifurcation vessel segment with a parent and two branch arteries of 8.8 mm length and 2.6 mm diameter each is covered by about 1.2 million cells, where there are 326 400 ECs, 848 640 SMCs, with a total of 1 175 040 cells. If each EC and SMC are modelled by four and five ODEs, respectively, then the total number of ODEs for this mesh is 10 575 360, while the larger mesh requires 20 901 888 ODEs. These numbers of ODEs need to be solved thousands of times even for a short temporal domain. Table II provides the details of the bifurcation meshes.

Parallel simulation implementation.
The simulations were performed on the IBM Blue-Gene P architecture at the University of Canterbury High Performance Computing Centre.
The implementation also has been tested on the BlueGene Q architecture at the Argonne National Laboratory.
For each quadrilateral domain, the system of ODEs was solved using an explicit Runge Kutta (4, 5) scheme initially implemented in Fortran as rksuite_90 component of the (NAG) Fortran Library [37]; in the work reported here, a version of rksuite_90 converted to C++ was used instead. 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 neighbouring domains. The inter-domain communication was implemented through the MPI library. The solution for inter-domain communication and the corresponding boundary conditions enforcement within a single MPI domain was implemented following the mesh ghost-cell communication pattern recommended for discretised domain decomposition [38]. For the simulations presented here, the inter-domain communication interval was set to 0.01 s; this interval was chosen to maintain the balance between MPI communication and computation times.
A parallel simulation requires as input the basic geometric configuration the surface for a given bifurcation in terms of the number of quads along axial and circumferential dimension; in addition, each simulation requires the configuration of the coupled EC and SMC in a quad. The quads are grouped into 2D Cartesian meshes, each representing a tubular segment/branch of a bifurcation. Within each 2D Cartesian mesh, the quads correspond to computational domains, which are mapped to distinct MPI tasks. In addition, each simulation requires an explicit agonist (ATP) map, where each EC cell is assigned a surface layer ATP value. The agonist maps were generated from the corresponding EC surfaces in order to preserve the implicit topology of the cells. The types of agonist maps used in this work are described in Sections 2.4 and 2.5.
The state variables' values for ECs and SMCs were written to HDF5 format after every elapsed physiological second. In the post-processing step, the HDF5 output data for every second was mapped onto the surface/cell mesh and written out in VTK's VTU format, which is an XML-based format with binary data for unstructured grid data. The visualisation and analysis were performed with the ParaView [39].

Synthetic agonist map
Work by [14] has shown that even in a time-dependent CFD-based solution, the ATP concentration can be adequately modelled by a time-averaged profile in areas exhibiting flow separation and low WSS. The areas characterised by large Péclet number (Pe) values (defined by the low diffusion coefficient of ATP) and thick concentration boundary layer of mass transfer are not significantly affected by the pulsatile flow. The ATP concentration in these areas can be represented by functions of a linear form. In addition, these areas are often preceded (upstream) and followed (downstream) by a constant WSS profile [14].
Synthetic agonist maps for the simulations presented here were generated with a sigmoid-like distribution of J PLC values in the endothelial layer defined as a function of axial distance only. As noted previously, ATP concentration is linearly related to the rate of IP 3 production given as J PLC . Figure 5 shows plot of a sigmoid function for a surface mesh in Figure 4(a). The properties of a synthetic agonist map generated in this manner enable a simulation for an arterial bifurcation where there exist both constant and linearly-varying ATP concentrations. Figure 6(a) shows the corresponding synthetic ATP map for number of domains Q D 4080.

CFD-based (physiological) agonist map
The CFD simulations for generating physiologically realistic ATP maps follow the algorithm described in the work of [14]. Flow and mass transfer fields were solved with the open-source CFD tool, OpenFOAM [40]. The flow field is solved iteratively via the continuity and momentum  Here, u represents the blood velocity vector field u D OEu; v; w, p the kinematic pressure field and the kinematic viscosity. The species mass transport was solved simultaneously via the conservation of species equation given by  where is the species ATP concentration and D is the isotropic diffusion coefficient. The walls of the artery are assumed to be rigid and stationary. Due to the very high Péclet number, the timeaveraged nucleotide concentration does not differ significantly from the steady state concentration [14]. The governing equations were solved iteratively via the Semi-Implicit Method for Pressure-Linked Equations algorithm for pressure-velocity coupling; the convective terms were handled with second-order upwinding for momentum and third-order upwinding Monotonic Upstream-Centred Scheme for Conservation Laws scheme for chemical species. A parabolic Dirichlet velocity profile was implemented on the inlet (the peak velocity was determined for the actual inlet area and the specified flow-rate) and zero velocity gradient Neumann boundary conditions on the outlets. The pressure field values were specified as zero-pressure gradient Neumann boundary conditions for the inlet and zero-pressure Dirichlet boundary condition on the outlets.
For the species boundary conditions, a constant Dirichlet value of 0.1 M was specified on the inlet and zero gradient Neumann boundary conditions on the outlets. As noted in Section 2.1, the relationship between ATP and IP 3 production rate is monotonic and increasing, which can be approximated by a linear function. Therefore, a linear relationship between ATP and the IP 3 flux into the EC (J PLC ) is deemed sufficient. Following the work of [17], a boundary condition is defined as follows: Here, S.x/ represents the production of ATP as a function of WSS in the form of: where the K and m values as reported in the work of [17]. In this case, the boundary condition on the endothelial surface is set to Robin (reactive) boundary condition specified on the arterial walls. This boundary condition is implemented iteratively by specifying a Neumann boundary condition for the ATP gradient, where the gradient itself is a function of ATP concentration. For every iteration, the boundary condition is evaluated based on the current ATP concentration, then the ATP gradient is specified on the wall, and the scalar transport equation is solved. The process is repeated until convergence is obtained. All other parameters for mass transport were taken from the work of [17], and the WSS is calculated from the Navier-Stokes flow-field. Michaels-Menten formulation used to determine the IP 3 concentration in the EC. Figure 11(a) shows the ATP concentration, resulting from CFD mass transport simulations for a mesh with the number of domains Q D 4080.

RESULTS
The simulations were performed with both synthetic and CFD-based ATP agonist surface concentration profiles. In both cases, the simulations were carried out for all four cell-coupling cases.   Figure 6 shows the agonist map for the Q D 4080 mesh along with the propagating Ca 2C oscillations in the SMC layer at two simulation steps of t=60 and 120 s, respectively. The oscillations originate in the downstream branches of the bifurcation at the upper end of the J PLC -dependent Ca 2C oscillatory domain § . As the oscillations propagate upstream, the wave amplitude increases, while the wave period decreases. This behaviour continues with the time progress of the simulation. § In the simulation snapshots the 'upstream' direction is defined from right to left (in opposition to blood flow) and 'downstream' to be from left to right (in the same direction as physiological blood flow).   Figure 7 zooms in on the upstream wave front, where the oscillations are dampened while propagating into the non-oscillatory domain at t=150 s. The phenomenon of wave propagation past the oscillatory domain is known to exist in excitable media. In this case, the Ca 2C waves travel past the oscillatory domain because the system is excitable for a range of low J PLC in the non-oscillatory domain close to the Hopf bifurcation. This excitability is found in similar Ca 2C SMC models such as the model by [41]. Spiral waves can form in excitable systems from open wave ends [42], which we see in Figure 7. Figure 8 shows the comparative visualisations for all coupling cases at the later stage of the simulations (t=600 s). This illustrates the differences in the oscillatory behaviour determined by the differing coupling configuration.

Synthetic agonist map
In order to generate 2D space/time plots for visualising and comparing temporal wave propagation patterns, a subset of SMCs was extracted from all time steps for a given simulation. This subset was defined by plotting time-dependent SMC concentrations along a line extending along the outer surface of the bifurcation from inlet to outlet. Figure 9 presents the comparison of the Ca 2C oscillation dynamics for the two mesh sizes. The 2D wave plots show that the overall simulation results show significant similarities for both bifurcation meshes. These differences are explained by the following: (1) The Q D 8064 mesh was not an exact scaling of the Q D 4080 mesh.   Figure 10 shows the zoomed-in versions of the space time plots for all four coupling cases. Significant differences exist for the wave profiles in each of the cases. The visualisations for all four coupling cases with the synthetic ATP map are available online ¶ .

Physiological agonist map
CFD-based agonist maps were generated only for the Q D 4080 bifurcation surface mesh described in Section 2.3.1. Figure 11 shows the agonist map for this surface along with the propagating Ca 2C oscillations in the SMC layer at time t D 300 s. The oscillations originate largely in the branches of the bifurcation where the J PLC concentration is high compared with the rest of the vessel segment. Figure 12 presents a comparison of simulation results with a CFD-based agonist map for all four coupling cases at t D 300 s. One commonality between all cases is that similar to the synthetic agonist map the oscillations propagate along the direction of the decreasing ATP gradient. Although not shown in the static figures, waves propagating on the surface of the inlet cylinder seem to do so in a circumferential direction. Variations in Ca 2C are determined from the small neighbourhood of the lower bound of the J PLC -dependent Ca 2C oscillatory range, as shown in Figure 2. The wave trains seen in the downstream branches move from right to left, along the gradient from high to low ATP   Figure 13 shows weak scaling results for four meshes of increasing size. Each simulation was performed for 1000 physiological seconds, which on average took 6.5 h. The results show a remarkable parallel scalability of the simulations with the increasing total problem size, while keeping the individual core workload constant. Figure 14 shows the breakdown of core hours between intradomain ODE solver time, inter-domain MPI communication time and IO-specific MPI aggregation time. This figure shows that with the given configuration of simulations (the number of cells per MPI domain, inter-domain communication interval and frequency of HDF5 IO) the IO-specific MPI aggregation time is optimal for the balance between computation and communication for this specific type of simulations.

DISCUSSION
Although the analysis of the model is qualitative in nature, the simulations results allow us to investigate the differences among coupling cases for a particular region in the ATP map. The visualisations (figures in the text as well as supplementary videos) are designed to demonstrate the emergent cell behaviour between healthy and pathological conditions. As it has been noted earlier, the four coupling cases represent simple approximations to a complex physiological phenomena and should not be seen as a complete description of the mechanisms involved. However, they are put forward as an initial framework from which further more complex pathways can be simulated. In particular, the role of IP 3 diffusion will require a close review due in part to its relatively small homotypic connexin channel permeability even though Ca 2C buffering is taken into account.

Synthetic agonist map
We first turn to discussion and analysis of Ca 2C waves for the synthetic agonist concentration maps. Figure 6 for Case 1 shows clearly the propagation direction of waves and that the wave length is ||   an increasing function of axial distance (length decreases as the wave moves from right to left). Initial analysis of the dynamics of the SMC indicates that this is due to reduced cytosolic Ca 2C whilst remaining in the oscillatory regime. Overall, waves do not propagate into areas of high ATP/J PLC /IP 3 concentrations defined by the upper bound of the J PLC -dependent Ca 2C oscillatory range shown in Figure 2. Figure 7 indicates three characteristic features of the wave propagation phenomenon. Firstly, that of variations in the width of the propagating wave as a function of circumferential direction. This is due to diffusion in both directions (axial and circumferential) where as stated before the width is a function of the flux of Ca 2C moving into neighbouring cells and being either taken up by the SR or diffused out to other nearby cells whose Ca 2C concentration is of a size that will allow ions to move across the cell boundary through gap junctions. The second feature is that of how waves terminate: as the concentration of ions in the cells decreases due to diffusion, the waves become weak and narrow in width and eventually are extinguished by seemingly curving back upon themselves. Spiral waves are known to occur in a number of biological conditions as reported in the work of [42]. Finally, the results show that the Ca 2C is dispersed in wave packets. The blocks of Ca 2C concentration are not single cells (because these are much smaller than a single mesh element). The initial examinations indicate that these packets of cells have similar cytosolic Ca 2C concentrations and are in phase in terms of their oscillatory dynamics. However, it is as yet unclear as to the exact reason for this. Figure 8 indicates the variations in wave propagation phenomena as a function of the four cases of coupling. Case 1 (Figure 8(a)) has no heterocellular Ca 2C coupling. Thus, any oscillations produced in the ECs do not acquire fully propagated to the SMCs, although because of the heterocellular IP 3 diffusion, cytosolic Ca 2C is propagated from cell to cell in only small amounts. It should be noted that if the sensitivity of the IP 3 receptors on the SR is increased then this may very well produce profiles similar to Case 2 because of increased cytosolic Ca 2C effluxed from the SR.
For Case 1, waves occur for a substantial distance downstream (left to right) in the daughter artery whereas for Cases 2, 3 and 4, the daughter arteries are relatively free of oscillation. For Case 2, compared with Case 1, the dynamics are substantially different due to the inclusion of heterocellular Ca 2C coupling. The agonist-induced increase in IP 3 and subsequent Ca 2C efflux from the ER is now diffused to the SMC via gap junctions. The increase in SMC Ca 2C causes oscillations (via the CICR mechanism) to occur, and further diffusion allows cells not normally participating in an oscillatory regime to become activated. The curvature of the waves seen in (Figure 8(b)) is a result of a complex dynamic. Initially, the waves propagate from the bifurcation upstream in an axisymmetric fashion, similar to Case 1. When the waves are weakened close to being extinguished, a reflected wave allows Ca 2C waves to move upstream. The curvature of the waves changes from positive to negative as time progresses. This can be seen in the video associated with the supplementary videos provided online. Case 3 (Figure 8(c)) is similar to Case 2, as expected because the difference between Cases 2 and 3 is that EC homocellular IP 3 coupling is removed in Case 3. It indicates once again that Ca 2C coupling is the dominant mechanism of wave propagation.
Case 4 coupling parameters results are shown in Figure 8(d); in this case, the wave dynamics are more chaotic and unstructured compared with the other coupling parameter cases. At the first stages of activation, the waves propagate upstream but become very short and rapidly disappear at a position just upstream from the Hopf bifurcation point denoted in Figure 2. As time progresses, the waves continue to move upstream and become extinguished further and further upstream from the original Hopf bifurcation point. At approximately t D 100 s, the upstream Ca 2C dynamics initiate very quick 'flashes' of intense concentration just downstream of the upstream wave domain border; these 'flashes' are evident for approximately 5 to 6 s and leave a band of inactive Ca 2C as shown in Figure 8(d) at t D 600 s. These can be seen in the supplementary video material provided online.
In the space/time concentration plots shown in Figure 9, the Ca 2C waves along the outer edge from inlet to upper daughter outlet indicate a strong similarity between the Q D 4080 (Figure 9(a)) simulation and the Q D 8096 simulation (Figure 9(b)). The smaller mesh is not an exact physical scaling from 4080 to 8064 mesh elements; thus, the J PLC maps are slightly different for the two mesh sizes. However, in both meshes the waves initially stop at the Hopf bifurcation point. After approximately 200 s, the waves move upstream with varying wave speed. The figures show that waves stop and start suddenly. The analysis of 1D waves (data not shown) indicates that this is caused by the diffusion of Ca 2C to the extent where the cytosolic Ca 2C is too low to cause the CICR mechanism to initiate oscillations. A finite time is needed for diffusion to build up the cell Ca 2C again to renew oscillatory behaviour. Figure 10 shows the space/time Ca 2C concentration plots for all four coupling cases. Case 1 coupling simulation results have been discussed earlier. Case 2 coupling plot indicates the effect of 'wave packets' propagating in both upstream and downstream directions. This points to distinct 'communities' of cells whose Ca 2C oscillations are in phase and hence do not provide sufficient diffusion between them. Only when the boundaries of these packets allow significant Ca 2C into the community does the cell packet breaks up. This phenomena need further analysis to confirm our initial explanation. Case 3, as mentioned previously is similar to Case 2. Case 4, however, presents a combination of phenomena where packets of cells exist throughout specific positions on the surface. At the later simulation stages, the waves propagate further upstream in a well-defined manner. For Cases 2, 3 and 4, the propagation beyond the bifurcation point starts immediately in contrast to Case 1, where it takes a substantial time to produce upstream wave propagation past into the nonoscillatory region. It is clear that the diffusion of Ca 2C through heterocellular gap junctions is the dominant mechanism for all wave profiles, although gap junction transmission of IP 3 does have an effect through the action of IP 3 receptors on the SR of an SMC. However, this may be a function of the IP 3 permeability value being relatively small.

CFD-based agonist map
We now turn to the case for the CFD-based ATP agonist map. Figure 11 shows the ATP surface profile derived from the CFD solution described earlier. Here, we see a substantially different ATP concentration from the synthetic map. Figure 11(b) shows the Ca 2C concentration at time t D 300 s. Although this ATP map is substantially different from the synthetic condition, the waves still propagate in a direction towards low ATP (or J PLC rate) concentration and, for this particular, ATP concentration profile waves move downstream into the daughter artery. Low ATP concentration occurs upstream in the parent artery segment. Here, complex waves and patterns occur, but this is due to the low levels of ATP.

Coupling and the transition from healthy to diseased conditions
A comparison with the simulation results for the CFD-based ATP map in Figure 12 demonstrates that these cases have markedly different ATP-defined oscillatory responses in the areas that are most likely to develop atherosclerotic plaque, that is the outer bends of a bifurcation surface used for the simulations in this work. The outer bends of the bifurcations under the assumed non-pathological coupling conditions (Cases 1 and 2) demonstrate lack of oscillatory behaviour. Under the pathological coupling conditions (specifically Case 4) the cells in the outer bends demonstrate a markedly different behaviour. While the exact transition path from healthy to diseased conditions remains to be established, it is reasonable to hypothesise that in an in vivo domain localised cell-coupling parameters could change in response to the haemodynamic-dependent agonist concentrations and stress on the cyto-skeleton over long periods of time. This agonist-induced change in coupling conditions could be one of the triggers linked with the transition from healthy to diseased state. This hypothesis requires further investigation with more complex coupling models. It may be suggested that the transition from healthy to diseased state is a multi-stage process where a response to the changing haemodynamic conditions results in coupling condition changes followed by the corresponding change in the oscillatory behaviour. Figure 12 illustrates this hypothesis. For example, the outer bend localised coupling could change over the course of time from healthy to pathological coupling, resulting in a different response to the agonist. For Case 1 (Figure 12), no oscillations occur in the outer bend (at any time during the simulation) whilst for Case 4, waves occur in the outer bend area for the majority of the simulation time, starting at approximately 100 s into the simulation. Further investigation of the response may help reveal the mechanics of arterial plaques growth.

Limitations and future work
Limitations of the model include only one-cell-deep layer of SMCs. Future work will extend the simulations to include a model of SMC layer of variable thickness in terms of the number of cells. Additionally, Figure 6(c) indicates that there is an apparent dependence of the wave speed (an increase) on the size of the quadrilaterals in the neighbourhood of the bifurcation. This is an artefact of the tessellation of the bifurcation surface. Around bifurcations, the area and dimensions of the surface quads are accurate only to the level of magnitude because of the complexity of the surface geometry. Some of the quads acquire stretched out along the first dimension and shortened along the second dimension. Thus, in this neighbourhood, the density of cells varies, because for each quadrilateral tessellating the arterial surface, there are a specific number of SMCs and ECs along each dimension. The simulation uses lumped parameter equations where the ion concentrations are constant throughout the cell. Consequently, waves appear to travel faster along the neck of the bifurcation. Although this is not physiologically accurate, diffusion inside the cell is ignored in the current model because of the assumptions of a lumped model. Future work will account for the improved distribution in cell size across the quadrilateral domains while allowing variations in ionic cytosolic concentrations. Most importantly, future projects will investigate atherosclerotic disease states transitions along with dynamic state transitions within homo-typic and hetero-typic coupling parameter space. This will involve explicit modelling of the effects of WSS on the endothelium as well as spatially varying coupling coefficients simulating the alignment of ECs in straight arteries compared with the random orientation and distribution of ECs found close to bifurcations and stagnation points defined by the blood fluid mechanics. Another branch of future work will extend the model to predict the sites of atherosclerosis initiation and compare the results of simulations to image-based clinical data.

CONCLUSIONS
We have developed a massively parallel framework for performing 3D simulations of coupled EC/SMC dynamics in a bifurcating artery. The simulations reported in this work include up to 2.3 million cells modelled by more than 20 million ODEs. The results of the time-dependent simulations show a radically varying range of Ca 2C wave propagation profiles determined by the coupling configurations. The simulations demonstrate that heterocellular diffusion of both IP 3 and Ca 2C through gap junctions plays an important role in information transfer over large spatial scales. Zero membrane potential coupling does not eliminate Ca 2C oscillations but provides a significant difference in the oscillatory dynamics over scales much larger than a single cell. Indeed the simulations show that IP 3 is an important transport mechanism for oscillations to occur and one that should not be neglected over and above that of Ca 2C .
The results of the reported simulations show that for both the synthetic and CFD-based agonist maps the phenomena of macro-scale Ca 2C oscillation propagation is strongly dependent on the ATP gradient along the EC layer and the Ca 2C cellular dynamics. Although the bifurcation meshes used in this work are symmetric, and the use of symmetric boundary conditions could have been used to reduce the computation time, the simulations were performed with full bifurcation meshes in preparation for large-scale simulations with meshes modelling realistic arterial segments where daughter arteries are of different diameters, and the bifurcation angles change depending on the location of the arterial segment in the human vasculature. We believe that the coupled EC/SMC dynamics modelling framework is an invaluable tool capable of providing deep insights into the cellular dynamics and the development of vascular disorders in connection with the local properties of arterial geometries. We would like to emphasise that the relationship between Ca 2C dynamics (in terms of cellular oscillations and wave propagation) and vascular disease is tentative at present, while considerable further research is required to advance the understanding of the interplay between vascular geometry, fluid dynamics and cellular function. It is our intention to extend the parallel implementation of the model to simulate patient-derived geometry at a realistic scale. This work will involve the use of substantially larger computational resource. The presented framework is a foundation step towards simulating a full coronary tree, which involves the solution of billions of ODEs, and is considered to be a classic exascale problem.

A.1. Coupled cells physiological model
A single pair EC/SMC model from the work of [18] accounts for the essential mechanisms of IP 3induced cytosolic calcium release and the cascade of the triggered processes. The work of [27] extends the single SMC and EC Koenigsberger model with the homocellular and heterocellular coupling. Similar to the intracellular ion concentrations and membrane voltage potential, the coupling terms are defined as ion fluxes; thus, the intercellular fluxes are added to the corresponding conservation equations for the intracellular concentrations and membrane voltage potential. The parameters of the EC/SMC model are provided in [27].
The following equations ((A.1)-(A.5)) specify the dynamics for cytosolic calcium concentration (c), SR calcium concentration (s), membrane potential (v), cytosolic IP 3 concentration (I) and the open state probability of the calcium-activated potassium channels (!) for the SMC layer.
ƒ‚ … SM C membrane potential terms Coupling membrane potential terms The following equations ((A.6)-(A.9)) specify the dynamics for cytosolic calcium concentration ( Q c), SR calcium concentration (Q s), membrane potential ( Q v) and cytosolic IP 3 concentration ( Q I) for the ECs layer in a similar manner to those for the SMC.
Coupling membrane potential terms The definitions of the intracellular fluxes for SMCs and ECs are provided in Appendix A.2 and Appendix A.3, respectively, while the definitions for the homocellular and heterocellular coupling terms are provided in Appendix A.4 and Appendix A.5. The complete list of algebraic constants in the SMCs and ECs equations can be found in the work of [18].
The terms specific to a single SMC or EC (intracellular dynamics) are detailed in the work of [18], while the coupling terms (intercellular dynamics) originate from the work of [27]. The complete list of algebraic constants in the SMCs and ECs equations can be found in the work of [18].

A.2. Intracellular SMC fluxes
(A.10) Ca 2C -induced Ca 2C release from SR into cytosol: Ca 2C transfer back into the SR via the Ca 2C -activated SERCA pump: Efflux of Ca 2C via the Na C /Ca 2C exchanger: Removal of intracellular Ca 2C through the plasma membrane Ca 2C -ATPase (PMCA) pump: Function of the Ca 2C -activated K C channels determined by the cytosolic Ca 2C concentration: Potassium efflux through the plasma membrane-bound BK Ca channels, where ! is the open channel probability of the Ca 2C -activated channels expressed as a function of K activation : Influx of chloride ions caused by the plasma membrane depolarisation: Constant rate efflux of potassium through the Na C /K C pump: Linear IP 3 concentration degradation:

A.3. Intracellular EC fluxes
Ca 2C -induced Ca 2C release from ER into cytosol: Ca 2C transfer back into the ER via the Ca 2C -activated SERCA pump:  Table A.I. SMC model parameters as described by [18].  38) In the above-mentioned equations, the k and l subscripts are the indices of the adjacent SMCs and ECs, respectively.

A.5. Heterocellular coupling
The work of [27] provides the following definitions for the heterocellular coupling components of the conservation equations for an SMC and its adjacent ECs include Ca 2C (A.39) and IP 3 In the equations above, the m and n subscripts are the indices of the ECs adjacent to a given SMC and SMCs adjacent to a given EC, respectively.