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

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.

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 arterial wall consists of three layers of tissue: tunica intima, tunica media and tunica adventita. The intima consists of a single layer of ECs; the media consists of the SMCs, and the adventita is made up of connective tissue. The endothelial layer plays the role of a mediator in the transfer of chemical species to and from the surrounding tissue. The mass transport phenomenon where chemical species are actively or passively diffused within large populations of cells is made possible through cellular communication channels. These channels are known as gap junctions; they are composed of connexin protein subunits and allow direct intercellular communication. For example, the Cx37 homotypic connexins are highly selective and allow only the transfer of monovalent cationic currents; thus, they have the most influence on the membrane potential 3. The Cx40 connexin protein subunit favours divalent cations, and thus, it is the major contributor to a gradient‐driven Ca2+ concentration current. The Cx43 connexin is the least selective, and it allows the passage of a range of large and small molecules such as ATP and inositol triphosphate (IP3) 4. The Cx43 connexin is normally present only in ECs although it has been reported to appear at the shoulder regions of the atherosclerotic lesions 5.

Recent publications on the physiology of vascular system dynamics hypothesise that cellular Ca2+ concentration oscillations (spreading both downstream and upstream) triggered by low WSS play an important role in atherosclerotic plaque formation 6, 7, 8, 9. The outcomes of the published in vivo and in vitro experiments provide descriptive models that use cellular‐level observations, but these observations fall short of explaining the macro‐scale phenomena of plaque formation and development 10. Atherosclerotic plaque lesions have characteristic‐length scales much larger than a single cell in the arterial wall. Hence, the mechanisms of plaque formation must be studied over considerable distances, such as hundreds to thousands of cell lengths.

Emergent macro‐scale behaviour

On the macro‐scale level, local geometric properties of arterial trees such as curvature, branching and bifurcations determine the local haemodynamics properties such as pressure, flow velocity and WSS. Among those properties, WSS carries a special significance as it influences the localisation of atherosclerotic lesions 11. Recent studies have established a link between atherosclerosis and low WSS, resulting from disturbed flow 12, 13. The magnitude of WSS determines the ATP concentration in the luminal boundary layer, which is in constant contact with the endothelium 14, 15. ATP is often described as the chemical energy transporter in cellular metabolism, and in the case of the vascular dynamics, it plays the role of a vasoconstrictor. In other studies, it has been proven that reduced ATP concentration is correlated with low WSS 16. Attempting to explain the onset of atherosclerosis, a study by 17 has proposed that atheroprotective mechanisms in a vascular wall are inhibited by reduced IP3 and Ca2+ signalling induced by low WSS.

The work of 9 presents a comprehensive example of cellular dynamics modelling, focusing on Ca2+ wave propagation in the SMC layer of arterial wall. This example of cellular dynamics modelling is aligned with the work presented here because 9 take into consideration the properties of WSS 8. However, the research presented in this article differs in two major ways from the approach of 9. Firstly, the model reported here is a coupled cells model, which includes both SMC and EC layers, and is therefore physiologically correct and able to interpret the transduction properties of the endothelial layer. Secondly, the macro‐scale phenomena of cell coupling are modelled on the micro‐scale where the ECs and SMCs in the arterial wall have physiologically accurate dimensions.

This study investigates temporal and spatial cellular dynamics of large populations of coupled ECs and SMCs with the view of gaining insight into the Ca2+ signalling in arterial wall tissue. The basic EC/SMC coupled model 18 has been validated through in vivo experiments in the work of 19. The experiments were performed with SMC arterial tissue strips of 750 × 50μm, which were equivalent to populations of approximately 140 SMCs. The in silico experiments reported in the work of 20 were performed with the population of about 400 SMCs. Our simulations included the populations of 850 000 SMCs and 330 000 ECs, which total is approximately 1.2 millioncoupled cells. These simulations provide new evidence over much larger scales equivalent to that found for atherosclerotic plaques.

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 IP3‐mediated channel. IP3 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 Ca2+ 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.

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:

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 IP3‐induced cytsolic Ca2+ 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.

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.

PIP2 gives rise to IP3 that is then released in the intracellular space. This nascent IP3 binds to IP3 receptor (IP3) on the ER/SR surface.

IP3‐bound IP3 induces release of Ca2+ ions from the ER/SR into the cytosol.

The Ca2+ release from intracellular store sensitises the IP3 further, which releases more Ca2+, thus making a Ca2+‐rich domain in the cytosol in both EC and SMC. The excess of intracellular Ca2+ depolarises the membrane potential.

ER/SR has low‐affinity binding sites on the cytosolic side of a channel, which constitutes a pump called the (SERCA) pump. Cytosolic Ca2+ encourages the replenishment of the intracellular stores via this pathway.

Ca2+ leaks from ER/SR consistently under a concentration gradient between cytosolic and ER/SR luminal Ca2+ and keeps the Ca2+ in equilibrium during a non‐stimulated state of the cell.

In an EC, the cytosolic Ca2+ favours the influx of extracellular Ca2+ from non‐selective cation channels.

The SERCA pump pushes out cytosolic Ca2+ to extracellular space.

10) In ECs, activation of KCa, upon binding to Ca2+ ions intracellularly at BKCa and SKCa, let K+ ions move out of the cytosol. This hyperpolarises the membrane potential.

Although K+ ions efflux is the main repolarising current, residual current (mainly consisting of monovalent ions) also contributes to membrane potential repolarisation.

The IP3 concentration increases in SMC cytosol via transmission of IP3 from coupled EC. This IP3 increase activates the downstream IP3‐induced Ca2+ release.

The IP3 induced, and calcium‐induced calcium release (CICR) Ca2+ depolarises the membrane potential.

The membrane depolarisation results in the influx of Ca2+ from VOCCs, which will close upon repolarisation in the following steps.

Ca2+, in addition to other pathways, is pushed out via Na/Ca exchanger.

Binding of Ca2+ ions to KCa opens BKCa channels in SMC causing K ions efflux and membrane repolarisation.

Influx of Cl2 ions add to the repolarisation.

Medium for intercellular communication via homocellular gap junctions can either be Ca2+, IP3 or membrane potential coupling.

Heterocellular gap junctions can couple an EC and an SMC via Ca2+, IP3 or membrane potential coupling. Hyperpolarised EC membrane potential can hyperpolarise SMC plasma membrane and consequently close VOCCs.

The ATP to J conversion models reported by 28 and 29 show that the relationship between ATP and IP3 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 IP3 into the EC, that is J. The range of ATP concentration is represented by the equivalent values of the J; the flux of IP3 into the EC is determined by the concentration of ATP in the vessel lumen. The whole range of J consists of three bands, which correspond to the ATP‐dependent Ca2+ dynamics in an SMC. As shown in Figure 2, low and high bands of J result in steady state Ca2+ concentrations, while the middle band shows the presence of Ca2+ oscillations. The oscillation envelope for this band indicates the relative high and low amplitudes for lower and higher J concentrations.

Figure 2

Bifurcation diagram for cytosolic smooth muscle cell (SMC) Ca2+ concentration in a single coupled endothelial cell/smooth muscle cell (SMC) unit as function of increasing J at the endothelial cell surface for coupling cases 1–4.

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 up‐regulation 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 Ca2+ and IP3 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 1 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 Ca2+ coupling. The remaining two coupling configurations are designed to simulate pathological cases. In Case 3, I P 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 Ca2+ 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 IP3 diffusion will need careful consideration due in part to its assumed smaller connexin channel permeability.

Table 1

Coupling cases defining healthy and pathological modes of intercellular communication.

	Homocellular		Heterocellular
Number	SMC	EC	SMC ↔ EC	Case Description
1	V, Ca2+, I P 3	V, Ca2+	V, I P 3	Non‐pathological coupling.
2	V, Ca2+, I P 3	V, Ca2+	V, I P 3, Ca2+	Non‐pathological coupling case with the
				added heterocellular Ca2+ coupling.
3	V, Ca2+, I P 3	V, Ca2+, I P 3	V, I P 3, Ca2+	Pathological case simulating lesion‐
				prone areas with the added I P 3 coupling
				in ECs.
4	V, Ca2+, I P 3	I P 3	I P 3	Pathological case simulating progressive
				atherosclerosis by removing homocellular
				(ECs) and heterocellular (SMCs
				and ECs) membrane potential and Ca2+
				coupling.

Figure 2 shows the bifurcation diagrams (maximum/minimum Ca2+ concentration vs the IP3 flux parameter J); for case 1, the range of IP3 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 Ca2+ 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,z) and f(u,v) can be either zero or some a priori defined function. For the present, we use f(u,v) = 0, where u and v are orthogonal coordinates, which provide axial and circumferential directions. The Neumann and Dirichlet boundary conditions for Eq. (1) define the rate at which the surface moves away from the boundary and its direction, respectively. Eq. (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. Figure 3(a) shows the boundaries of the semi‐tubular segments comprising the 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.

Figure 3

Composition of a bifurcations surface from semi‐tubular segments (top) with a single patch boundary conditions for a semi‐tubular segment (bottom).

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,z) space.

Figure 4(a) shows a bifurcation surface mesh consisting of Q = 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.

Figure 4

An example of a bifurcation surface mesh decomposed into Q = 4080 computational domain (top); arrangement of smooth muscle cells (SMCs) and endothelial cells (ECs) in a single quadrilateral mapped to a computational task in parallel environment (bottom). ODE, ordinary differential equation.

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.

Bjorstad, P. E. (1980). Numerical Solution of the Biharmonic Equation. Retrieved from http://http://www.netlib.no/netlib/bihar/dbihar.f.

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 326400 ECs, 848640 SMCs, with a total of 1175040 cells. If each EC and SMC are modelled by four and five ODEs, respectively, then the total number of ODEs for this mesh is 10575360, while the larger mesh requires 20901888 ODEs. These numbers of ODEs need to be solved thousands of times even for a short temporal domain. Table 2 provides the details of the bifurcation meshes.

Table 2

Parameters of the arterial surface meshes

Number	L B	L Max	D	α B	Q	ECs	SMCs	Total
1	8.84	17.68	2.55	30	4,080	326,400	848,640	1,175,040
2	10.92	21.84	4.08	30	8,064	645,120	1,677,312	2,322,432

Arterial surface meshes include branch length L (mm), maximum length L (mm), diameter D (mm), branch centreline angle α (degrees), number of domains Q, number of ECs, number of SMCs and the total number of cells in a vessel bifurcation segment; each domain for both mesh variants included a total of 288 cells, with 208 and 80 SMCs and ECs, respectively.

Parallel simulation implementation

The simulations were performed on the IBM BlueGene 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 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 IP3 production given as J. 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 = 4080.

Figure 5

Sigmoid function used in the modelling of the adenosine triphosphate distribution for a synthetic agonist map; adenosine triphosphate concentration is linearly related to the rate of inositol trisphosphate (IP3) production given as J; the function is centred on the bifurcation point; the negative endothelial cell positions are decreasing towards the inlet of the bifurcation segment, while the positive numbers are increasing towards the outlets. EC, endothelial cell.

Figure 6

Synthetic adenosine triphosphate (ATP) concentration map and a visualisation of Ca2+ waves for the Q = 4080 mesh under Case 1 coupling conditions.

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 Eqs. (2) and (3), on the assumption that the blood flow is steady state, incompressible and of Newtonian rheology:

Here, u represents the blood velocity vector field u = [u,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 time‐averaged 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 IP3 production rate is monotonic and increasing, which can be approximated by a linear function. Therefore, a linear relationship between ATP and the IP3 flux into the EC (J ) 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 τ 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 IP3 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 = 4080.

Figure 11

CFD‐based adenosine triphosphate concentration map and a visualisation of Ca2+ waves.

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. All simulation runs were executed for the Q = 4080 mesh (corresponding to 4080 cores of Blue Gene/P). For the Q = 8064 mesh (corresponding to 8064 cores of Blue Gene/P), only Case 1 coupling with a synthetic agonist profile was performed. This simulation was designed to verify that a larger population of cells did not change the dynamics of the simulations in a significant way. All simulations were run for 1000 physiological seconds, where output was recorded after every physiological second. The small mesh bifurcation simulations with Q = 4080 cores on Blue Gene/L on average took approximately 20.3 h to complete, while the large mesh bifurcation simulations with Q = 8064 cores on Blue Gene/P took approximately 15.1 h. The visualisations for all simulations are available on the UC HPC YouTube Channel.

Figure 6 shows the agonist map for the Q = 4080 mesh along with the propagating Ca2+ 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‐dependent Ca2+ 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 Ca2+ waves travel past the oscillatory domain because the system is excitable for a range of low J in the non‐oscillatory domain close to the Hopf bifurcation. This excitability is found in similar Ca2+ 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 7

Propagating wave‐front at t = 150 s for coupling Case 1 with the 4K mesh. SMC, smooth muscle cell.

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.

Figure 8

Wave plots for the synthetic agonist map for all coupling cases at t = 600. SMC, smooth muscle cell.

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 Ca2+ 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:

Figure 9

Comparison of the Ca2+ (μ M) waves for Case 1 coupling in the SMC, smooth muscle cell layer along the extracted SMC, smooth muscle cells line; the red line indicates the lower bound of the J‐dependent Ca2+ oscillatory range, where the oscillations are ‘extinguished’ at the start of the simulation, as the simulation progresses Ca2+ waves tend to spread into the non‐oscillatory domain.

The Q = 8064 mesh was not an exact scaling of the Q = 4080 mesh.

The Hopf bifurcation point (lower bound of oscillations for a single SMC) in the ATP profile for the Q = 8064 bifurcation is at a different location in the parent artery segment compared with the Q = 4080 mesh.

The Ca2+ oscillations in the SMC under the Case 1 coupling conditions are initially ‘extinguished’ at the lower border of the oscillatory J concentration range. However, in both cases the oscillations propagate upstream past the lower boundary (red horizontal line in the figure) at approximately 200 and 250 s for the Q = 4080 and Q = 8064 meshes, respectively.

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 .

Figure 10

Comparison of zoomed‐in plots showing the Ca2+ (μ M) waves in the SMC, smooth muscle cell layer along the extracted SMCs line; the red line indicates the lower bound of the J‐dependent Ca2+ oscillatory range, shown originally in Figure 2.

Coupling Case 1—https://www.youtube.com/watch?v=ax8I6rZG-6Q,Coupling Case 2—https://www.youtube.com/watch?v=HfMHsZrM8bw,Coupling Case 3—https://www.youtube.com/watch?v=lx-Pza_{Z}HVk,Coupling Case 4—https://www.youtube.com/watch?v=5JpfDcs5GTA.

Physiological agonist map

CFD‐based agonist maps were generated only for the Q = 4080 bifurcation surface mesh described in Section 2.3.1. Figure 11 shows the agonist map for this surface along with the propagating Ca2+ oscillations in the SMC layer at time t = 300 s. The oscillations originate largely in the branches of the bifurcation where the J 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 = 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 Ca2+ are determined from the small neighbourhood of the lower bound of the J‐dependent Ca2+ 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 concentration. For Case 1, oscillations occur at the inner surface of the bifurcation point whilst for all of the other three cases, this area is devoid of oscillations and holds a constant Ca2+ concentration. The visualisations for all four coupling cases with the CFD‐based ATP map are available online .

Figure 12

Visualisation of Ca2+ oscillations for all four coupling cases for the computational fluid dynamics‐based agonist map at time step t = 300.

Coupling Case 1—https://www.youtube.com/watch?v=UPZZdbF8De0,Coupling Case 2—https://www.youtube.com/watch?v=QH9gyH5s1fo,Coupling Case 3—https://www.youtube.com/watch?v=yj_{q}1QIbrdg,Coupling Case 4—https://www.youtube.com/watch?v=8q6hrFsUy4s.

Parallel scalability

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 intra‐domain 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.

Figure 13

Weak scaling results for four meshes of increasing size (1000 physiological seconds simulations).

Figure 14

Simulation time breakdown for variable number of IO‐dedicated cores (1000 physiological seconds simulations); HDF5 output time is not shown due to negligible values. MPI, message passing interface; ODE, ordinary differential equation.

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 IP3 diffusion will require a close review due in part to its relatively small homotypic connexin channel permeability even though Ca2+ buffering is taken into account.

We first turn to discussion and analysis of Ca2+ 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 Ca2+ whilst remaining in the oscillatory regime. Overall, waves do not propagate into areas of high ATP/J/IP3 concentrations defined by the upper bound of the J‐dependent Ca2+ 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 Ca2+ moving into neighbouring cells and being either taken up by the SR or diffused out to other nearby cells whose Ca2+ 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 Ca2+ is dispersed in wave packets. The blocks of Ca2+ 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 Ca2+ 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 Ca2+ coupling. Thus, any oscillations produced in the ECs do not acquire fully propagated to the SMCs, although because of the heterocellular IP3 diffusion, cytosolic Ca2+ is propagated from cell to cell in only small amounts. It should be noted that if the sensitivity of the I P 3 receptors on the SR is increased then this may very well produce profiles similar to Case 2 because of increased cytosolic Ca2+ 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 Ca2+ coupling. The agonist‐induced increase in IP3 and subsequent Ca2+ efflux from the ER is now diffused to the SMC via gap junctions. The increase in SMC Ca2+ 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 Ca2+ 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 IP3 coupling is removed in Case 3. It indicates once again that Ca2+ 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 = 100 s, the upstream Ca2+ 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 Ca2+ as shown in Figure 8(d) at t = 600 s. These can be seen in the supplementary video material provided online.

In the space/time concentration plots shown in Figure 9, the Ca2+ waves along the outer edge from inlet to upper daughter outlet indicate a strong similarity between the Q = 4080 (Figure 9(a)) simulation and the Q = 8096 simulation (Figure 9(b)). The smaller mesh is not an exact physical scaling from 4080 to 8064 mesh elements; thus, the J 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 Ca2+ to the extent where the cytosolic Ca2+ is too low to cause the CICR mechanism to initiate oscillations. A finite time is needed for diffusion to build up the cell Ca2+ again to renew oscillatory behaviour.

Figure 10 shows the space/time Ca2+ 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 Ca2+ oscillations are in phase and hence do not provide sufficient diffusion between them. Only when the boundaries of these packets allow significant Ca2+ 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 non‐oscillatory region. It is clear that the diffusion of Ca2+ through heterocellular gap junctions is the dominant mechanism for all wave profiles, although gap junction transmission of IP3 does have an effect through the action of IP3 receptors on the SR of an SMC. However, this may be a function of the I P3 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 Ca2+ concentration at time t = 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 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 Ca2+ wave propagation profiles determined by the coupling configurations. The simulations demonstrate that heterocellular diffusion of both IP3 and Ca2+ through gap junctions plays an important role in information transfer over large spatial scales. Zero membrane potential coupling does not eliminate Ca2+ oscillations but provides a significant difference in the oscillatory dynamics over scales much larger than a single cell. Indeed the simulations show that IP3 is an important transport mechanism for oscillations to occur and one that should not be neglected over and above that of Ca2+.

The results of the reported simulations show that for both the synthetic and CFD‐based agonist maps the phenomena of macro‐scale Ca2+ oscillation propagation is strongly dependent on the ATP gradient along the EC layer and the Ca2+ 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 Ca2+ 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.

Acknowledgements

The authors would like to thank the members of the University of Canterbury High Performance Computing Centre, Tony Dale and François Bissey for their enthusiastic support and expert advice on the technical front. We also thank Allanah Kenny for the dynamical system behaviour analysis and the final manuscript review.

The authors wish to acknowledge the contribution of NeSI to the results of this research. New Zealand's national compute and analytics services and team are supported by the New Zealand eScience Infrastructure (NeSI) and funded jointly by NeSI's collaborator institutions and through the Ministry of Business, Innovation and Employment.

Supporting info item

Click here for additional data file.

Supporting info item

Click here for additional data file.

Table A.1

SMC model parameters as described by 18.

F	Maximal rate of activation‐dependent calcium influx	0.23 μ Ms−1
K r	Half saturation constant for agonist‐dependent calcium entry	1 μ M
G Ca	Whole cell conductance for VOCCs	0.00129 μ MmV−1s−1
v Ca1	Reversal potential for VOCCs	100 mV
v Ca2	Half‐point of the VOCC activation sigmoidal	−24 mV
R Ca	Maximum slope of the VOCC activation sigmoidal	8.5 mV
G Na/Ca	Whole cell conductance for Na+/Ca2+ exchange	0.00316 μ MmV−1s−1
c Na/Ca	Half‐point for activation of Na+/Ca2+ exchange by Ca2+	0.5 μ M
v Na/Ca	Reversal potential for the Na+/Ca2+ exchanger	−30 mV
B	SR uptake rate constant	2.025 μ Ms−1
c b	Half‐point of the SERCA activation sigmoidal	1.0 μ M
C	CICR rate constant	55 μ Ms−1
s c	Half‐point of the CICR Ca2+ efflux sigmoidal	2.0 μ M
c c	Half‐point of the CICR activation sigmoidal	0.9 μ M
D	Rate constant for Ca2+ extrusion by the ATPase pump	0.24 s−1
v d	Intercept of voltage dependence of extrusion ATPase	−100 mV
R d	Slope of voltage dependence of extrusion ATPase	250 mV
L	Leak from SR rate constant	0.025 s−1
g	Scaling factor relating net movement of ion fluxes to the membrane potential (inversely related to cell capacitance)	1970 mV−1 μ M
F Na/K	Net whole cell flux via the Na+–K+–ATPase	0.0432 μ Ms−1
G Cl	Whole cell conductance for Cl− current	0.00134 μ MmV−1s−1
v Cl	Reversal potential for Cl− channels	−25 mV
G K	Whole cell conductance for K+ efflux	0.00446 μ MmV−1s−1
v K	Reversal potential for K+	−94 mV
λ	Rate constant for net KCa channel opening	45
c ω	Translation factor for Ca2+ dependence of KCa channel activation sigmoidal	0 μ M
β	Translation factor for membrane potential dependence of KCa channel activation sigmoidal	0.13 μ M 2
v_Ca_3	Half‐point for the KCa channel activation sigmoidal	−27 mV
R K	Maximum slope of the KCa activation sigmoidal	12.0 mV
k	Rate constant of IP3 degradation	0.1 s−1

Table A2

EC model parameters as described by 18.

F~	Maximal rate of activation‐dependent calcium influx	0.23 μ Ms−1
K~_r	half saturation constant for agonist‐dependent calcium entry	1 μ M
B~	ER uptake rate constant	0.5 μ Ms−1
c~_b	Half‐point of the SERCA activation sigmoidal	1.0 μ M
C~	CICR rate constant	5 μ Ms−1
s~_c	Half‐point of the CICR Ca2+ efflux sigmoidal	2.0 μ M
c~_c	Half‐point of the CICR activation sigmoidal	0.9 μ M
D~	Rate constant for Ca2+ extrusion by the ATPase pump	0.24 s−1
L~	Leak from SR rate constant	0.025 s−1
k~	Rate constant of IP3 degradation	0.1 s−1
G~_cat	Whole cell cation channel conductivity	0.66 μ M μ V −1s−1
E~_Ca	Ca2+ equilibrium potential	50 μ V
m~_3_cat		−0.18 μ M
m~_4_cat		0.37 μ M
J~_0	Constant calcium influx	0.029 μ Ms−1
C~_m	Membrane capacitance	25.8 pF
G~_tot	Total potassium channel conductivity	6927 pS
v~_K	K+ equilibrium potential	−80 μ V
ã_K1		53.3 μ M u V
ã_K2		53.3 μ V μ M −1
b~_K		−80.8 μ V
c~_K		−0.4 μ M
m~_3b		1.32 ×103 μ M μ V −1
m~_4b		0.30 μ M μ V −1
m~_3s		−0.28 μ M
m~_4s		0.389 μ M
G~_R	Residual current conductivity	955 pS
v~_rest	Membrane resting potential	−31.1 μ V

Table A.3

Coupling parameters for four intercellular communication configurations.

		Case 1	Case 2	Case 3	Case 4
Homocellular	g	1000	1000	1000	1000
	g~	1000	1000	1000	0
	p Ca	0.05	0.05	0.05	0.05
	p~_Ca	0.05	0.05	0.05	0
	p_IP_3	0.05	0.05	0.05	0.05
	p~_IP_3	0	0	0.05	0.05
Heterocellular	G	50.0	50.0	50.0	0
	G~	50.0	50.0	50.0	0
	P Ca	0	0.05	0.05	0
	P~_Ca	0	0.05	0.05	0
	P_IP_3	0.05	0.05	0.05	0.05
	P~_IP_3	0.05	0.05	0.05	0.05

Cases 1 and 2 simulate healthy paracrine communication, while Cases 3 and 4 simulate pathological cases, with early and progressive atherosclerotic lesions, respectively.