A Graphical Approach for Hybrid Simulation of 3D Diffusion Bio-Models via Coloured Hybrid Petri Nets

,ree-dimensional modelling of biological systems is imperative to study the behaviour of dynamic systems that require the analysis of how their components interact in space. However, there are only a few formal tools that offer a convenient modelling of such systems. ,e traditional approach to construct and simulate 3D models is to build a system of partial differential equations (PDEs). Although this approach may be computationally efficient and has been employed by many researchers over the years, it is not always intuitive since it does not provide a visual depiction of the modelled systems. Indeed, a visual modelling can help to conceive a mental image which eventually contributes to the understanding of the problem under study. Coloured Hybrid Petri Nets (HPN ) are a high-level representation of classical Petri nets that offer hybrid as well as spatial modelling of biological systems. In addition to their graphical representations, HPN C models are also scalable.,is paper shows how HPN C can be used to construct and simulate systems that require three-dimensional as well as hybrid (stochastic/continuous) modelling. We use calcium diffusion in three dimensions to illustrate our main ideas. More specifically, we show that creating 3D models using HPN C can yield more flexible models as the structure can be easily scaled up and down by just modifying a few parameters. ,is advantage of convenient model configuration facilitates the design of different experiments without the need to alter the model structure.


Introduction
Petri nets, as intuitive graphical tools, play an important role in modelling and simulating biochemical reaction networks as well as many other kinds of systems (see, e.g., [1][2][3]). To address the varying aspects of different modelling scenarios, many extensions of the standard place/transition nets have been proposed in the literature, most notably various quantitative net classes, such as Stochastic Petri Nets (SPN ) [4] and Continuous Petri Nets (CPN ) [5].
In addition, Hybrid Petri Nets (HPN ), that combine all features of SPN and CPN were proposed in [1,5]. HPN provide a flexible tool to integrate smoothly continuous and discrete states. Furthermore, coloured Petri nets (PN C ) [6,7] have been introduced to provide a high-level representation of standard and quantitative Petri nets by associating colours with places and transitions.
Finally, based on HPN and PN C , coloured Hybrid Petri Nets HPN C have been proposed in [8][9][10] combining all the merits of HPN and PN C into one net class. e different notions of firing rates of HPN facilitate the modelling of multitemporal scales, while the use of colours allows a parametrised modelling of large systems with repeated components. In addition, PN C can be converted algorithmically into equivalent low-level PN via an unfolding process [11]. Indeed, HPN C render it possible to model complex biological systems using a graphical and flexible graphical language.
A good example of complex biological systems that require substantial efforts to understand their interacting components is the dynamics of intracellular calcium. Intracellular calcium (Ca 2+ ) dynamics [12][13][14][15] is an important biological pathway that has been recently widely studied. Intracellular Ca 2+ works as an essential controller in organising some basic cellular functions such as gene expressions [12,13,16]. Due to its complexities, modelling and simulation plays an important role in the understanding of the pathway of Ca 2+ dynamics. So far, many models were introduced to analyse Ca 2+ dynamics in two-dimensional space (see, e.g., [14][15][16][17][18][19]), but only a few models for the threedimensional space (see, e.g., [20][21][22][23]). However, constructing an in silico Ca 2+ model cannot be achieved by using a single modelling paradigm (i.e., deterministic or stochastic).
On one hand, the Ca 2+ pathway comprises processes that have a huge number of variables and states (e.g., diffusion). Modelling such processes using a stochastic approach will require longer simulation time (and in some cases, it is even not possible to obtain results). On the other hand, certain processes are crucial for the successful reconstruction of clinical and "wet-lab" results. ese processes (e.g., channel opening and closing) completely rely on variables (channel states) with a very few number of states. Simulating these variables deterministically will not capture the correct model behaviour.
us, 3D hybrid modelling of this case study and similar biological systems can be of paramount importance to better understand and get more insights into the dynamic behaviour. A tool that provides the hybridisation of both time and space, whereby deterministic and stochastic processes can interplay with each other, will considerably simplify the construction of 3D models.
is paper utilises HPN C [10] to construct and simulate biological 3D models. More specifically, we show how continuous and stochastic components can be combined to create 3D models that can be employed to experiment with different behaviours. We apply this framework to model intracellular Ca 2+ behaviour. e resulting model is flexible since it permits adjusting the model parameters to allow experimenting with different simulation configurations. In our modelling approach, HPN capture the continuous/ discrete semantics, while colours encode channels, clusters, and diffusion grids (see below). e constructed model is validated by conducting a few experiments and comparing their outcome with facts known from the literature. e remainder of the paper is organised as follows: Section 2 recalls Petri nets and coloured (hybrid) Petri nets. Section 3 reviews the required background about Ca 2+ dynamics as well as the utilisation of HPN C for modelling biological processes. Section 4 shows how HPN C can be used to construct and execute biological models illustrated by the Ca 2+ case study. Simulation results of the Ca 2+ model are presented in Section 4. We conclude the paper by some closing remarks and an outline of possible future work.

Background
In this section, we summarise the use of Petri nets for modelling biological applications including standard Petri nets and their extension and hybrid Petri nets and coloured hybrid Petri nets.

Petri Nets.
Petri Nets (PN ) are an excellent modelling formalism for studying the behaviour of complex systems made of concurrent components. ey have been applied in many areas and applications including biological systems [24]. In terms of discrete mathematics, PN are directed, weighted, bipartite graphs consisting of transitions, places, and arcs. Arcs are directed from places to transitions, or vice versa. In systems biology, places may represent species, while transitions might be used to represent any kind of chemical reactions or biological processes.
SPN extend PN by associating each transition with a random firing rate. e semantics of an SPN model is a Continuous Time Markov Chain (CTMC); therefore, it can be simulated using a Stochastic Simulation Algorithm (SSA) [26]. Besides, extended stochastic Petri nets (XSPN ) were proposed by considering different transitions types and different arcs types [27,28]. Moreover, in order to model continuous concentration changes of biological processes, CPN were proposed in [1]. e semantics of a CPN model corresponds to a system of Ordinary Differential Equations (ODEs) [29]. Furthermore, in order to incorporate both discrete and continuous modelling capabilities in one model, HPN were proposed in [1] and further augmented to Generalised HPN (GHPN ) in [25]. e motivation of GHPN , compared to HPN , is to support more timing features, including stochastic, discrete deterministic, and continuous deterministic firing of transitions [25]. GHPN simulation approaches combine stochastic simulation algorithms with ODE numerical integrators, employing appropriate mechanisms to synchronise the stochastic and deterministic firings of biological events [25,30]. More specifically, in GHPN , places are further classified into two categories: discrete and continuous. Discrete places hold nonnegative integer numbers, called tokens, while continuous places hold nonnegative real numbers that are often referred to as place marking in order to distinguish them from the discrete counterpart. Moreover, transitions are classified as continuous, stochastic, immediate, deterministically delayed, and scheduled ones [25]. Immediate transitions fire instantly with zero delay and have higher priority than any other transition types; they are useful to model certain system semantics. Deterministically delayed transitions fire after a deterministically specified value. Scheduled transitions are a special type of deterministically delayed transitions where the start and end time of the repeated firing is given in advance.
GHPN provide such special arcs as read, inhibitor, equal, modifier, and reset arcs, which always go from a place to a transition. eir purpose is to either impose certain conditions on the transition firings or allow the use of a preplace in the definition of the transition rate function. For example, if a transition is connected with a preplace via a read arc, then this transition cannot be fired unless the number of tokens in the preplace is greater than or equal to the arc weight. On the contrary, a transition connected with a preplace via an inhibitor arc cannot be fired unless the number of tokens in the preplace is less than the arc weight. Moreover, equal arcs require the number of tokens to be exactly the same as the arc weight. However, modifier arcs do not impose any constraint on transition firing; they are used to allow a place to occur in the transition rate function. Finally, reset arcs reset the preplace marking to zero when the corresponding transition fires. e precise semantics of these arcs as well as the syntactic rules that govern their use are detailed in [25]. Please note that, during this paper, we use HPN and GHPN interchangeably unless it is explicitly stated.

Coloured Hybrid Petri
Nets. Coloured Petri Nets (PN C ) extend low-level PN with a set of finite colour sets [6,7]. PN C offer a parameterised approach to construct large biological systems; there are many coloured net classes, among them coloured continuous Petri nets (CPN C ) [7] and coloured stochastic Petri nets (SPN C ) [31].
By augmenting GHPN with colours, we obtain a new class of hybrid Petri nets called Coloured Generalised Hybrid Petri Nets (GHPN C ) [8]. In GHPN C , places are annotated with colour sets and place markings are distinguished by colours. Similarly, transitions are extended by assigning guards to them. Guards control transition enabling in addition to the traditional enabling rules of uncoloured Petri nets. Arcs are also augmented with colour expressions which define the set of arc instances over the colour sets of the corresponding places.
GHPN C permit the modelling of systems that exhibit multiple temporal scales [25] (by supporting stochastic and continuous semantics) and also spatial scales [32] (by the help of colours). Figure 1 provides a simple example, based on the general framework presented in [32], illustrating the use of CPN C to model diffusion in 3D space.  ese six neighbours are calculated by (l ± 1, m, k), (l, m ± 1, k) or (l, m, z ± 1), while Figure 1(b) describes a component at the index (l, m, k) and its 26 neighbours at different indices. ese 26 neighbours are calculated by (l ± 1, m ± 1, k ± 1), (l ± 1, m ± 1, k ± 1), (l ± 1, m ± 1, k ± 1) or (l ± 1, m ± 1, k ± 1). To encode these components in 3D, we need a framework having the ability to define the space as in Figures 1(a) and 1(b). HPN C is a framework that has the ability to do that. To encode the low-level components of the model, our approach does not only rely on the structure of the HPN C , but also employs their colour expressions. To present a component we use the graph structure of HPN C by creating only one place connected with only one transition as shown in Figure 1(c), while the space and the other components of the system are defined by a neighbourhood function.
For instance, in Figure 1, we define 3D space with colour sets. x-dimension, y-dimension, and z-dimension are defined by the colour sets xdim, ydim, and zdim, respectively. Based on the previously defined colour sets to encode twodimensional space [33], we define, in this paper, a threedimensional compound colour set, Grid 3D, representing a cube form. Moreover, we define six variables, a, x of type xdim, b, y of type ydim, and c, z of type zdim, as well as three-dimensional coordinates, (x, y, z) to address every component in space. Each of these components has a set of neighbours which are defined by the function IsNeighbour(x, y, z, a, b, c): bool IsNeighbour (xdim x, ydim y, zdim z, xdim a, ydim b, zdim c) where xdim, ydim, and zdim are the same as in Figure 1(d).
e function IsNeighbour () returns either true or false depending on the evaluation of the function body. e list of variables, constants, and function is given in Figure 1(d). Furthermore, in Figure 1, we represent the components of the system by the place S. We connect the place S with the transition Diffc modelling the diffusion to all positions fulfilling the guard IsNeighbour.
In our developed HPN C model, we make use of the component presented in Figure 1 as a building block for modelling calcium diffusion in 3D. is component is indeed a compact representation of the diffusion process, which, when unfolded, will reproduce the exact model dynamics.
In summary, the key modelling features of GHPN C that can be helpful for constructing 3D models of biological systems are colours and the neighbouring function. Colours define component instances, while the neighbouring functions define how these instances are related to each other. Moreover, the discrete and continuous components of GHPN C permit to combine different timescales.

Intracellular Calcium Dynamics.
is section provides an example of how GHPN C can be employed to construct and simulate 3D models. e case study is about the 3D modelling of intracellular calcium dynamics. is specific example is selected because it is intricate to model it via existing modelling approaches. Moreover, this case study is Modelling and Simulation in Engineering highly suitable to demonstrate all required features to construct and simulate 3D systems that require hybrid and spatial characteristics. In what follows, a brief overview of the calcium dynamics is presented. Afterwards, our modelling approach is applied to this case study. We commence with an overview of the required background of Ca 2+ dynamics, followed by a detailed discussion of the GHPN C model construction.

Calcium
Release. Intracellular calcium is organised by the release of Ca 2+ ions from the cell membrane or internal stores such as Endoplasmic Reticulum (ER). Releasing of Ca 2+ from ER occurs via the opening of inositol-4,5-triphosphate receptor (IP 3 R) channels (among different types, such as RyR and L-type calcium channels [13]). ese channels are randomly grouped into clusters [13]. A stochastic model to study the random opening and closing of IP 3 R channels is presented in [17]. e stochastic behaviour of channel opening leads to a release event in the form of calcium pulses that elapse in a few milliseconds [16]. Figure 2 depicts the different components involved in Ca 2+ release from the ER to cytosol.
As it is detailed in [17], (IP 3 R) channels consist of four identical and independent subunits which can be combined with the deterministic part to study the Ca 2+ dynamics (see, e.g., [15]). Every subunit of the channel can be in any of eight different states and has three distinct binding sites: one for activation of IP 3 , one for activating Ca 2+ , and a third for Ca 2+ inhibiting. Each of these sites is either busy or idle by their    Figure 2: A schematic diagram describing the different processes and components involved in the Ca 2+ model according to [19,20,38]. When some channels are open, Ca 2+ flows from the ER to the cytosol. Ca 2+ returns back to the ER via a process called pump. Channel opening is affected by the Ca 2+ concentration surrounding the channel. Ca 2+ in cytosol can be associated/dissociated with different types of buffer. Ca 2+ in cytosol, ER, and Buffers rapidly diffuses to the neighbouring cells. particular regulatory agonist (Ca 2+ or IP 3 ) [17,34]. Building on the experimental studies in [35], Marchant et al. observed that in the absence of bound Ca 2+ , IP 3 does not bind to its activating site. erefore, the two states of unbound IP 3 can be neglected. Hence, the eight-state model can be reduced to six states without changing the mechanism of the channel model. Furthermore, in [36], the eight states of each channel are reduced to just four states, based on the faster time scale and the condition of equilibrium; each two states of the eight states can be grouped to only one state. In more details, there exist three distinct time scales: binding of IP 3 (the fastest), binding by Ca 2+ activation, and binding by Ca 2+ inhibition (the slowest). As IP 3 dissociation and association is faster than the competing processes, transitions between states of unbound IP 3 and bound IP 3 are the first to reach the equilibrium state. erefore, the unbound IP 3 state is grouped with the bound IP 3 state to form only one state.
In [33], we have adopted a model for the stochastic part that accounts for the eight states of each subunit. However, for the purpose of constructing a 3D model, we employ the reduced system for the interactions of the calcium subunits from [37], which exhibits less subunit states. e transition which connects two subunit states can be naturally modelled via a simple Markovian process, making it very suitable to use stochastic transitions in the Petri net notations. In [38], another state has been added to the four states of each subunit to denote the conformational change of a subunit (i.e., to denote the number of active subunits). Opening a channel occurs when three or more subunits are in this 5 th state.
Finally, the opening of an IP 3 R channel is influenced by the neighbouring Ca 2+ concentration in the cytosol. e chance of opening a channel is increased by opening neighbouring channels [16].
2.6. Calcium Diffusion. We aim to model and simulate the calcium dynamics graphically in the form of a cube. Mathematically, the dynamics of calcium can be described by the diffusion process via a system of partial differential equation. e diffusion process plays an important role in the movement of calcium in ER and cytosol. In ER, it affects the opening of channels as well as calcium releasing, while it increases cytosolic calcium concentration.
Calcium diffusion was extensively studied in [15,19,38] to describe the spatio-temporal behaviour of Ca 2+ dynamics in two-dimensional space, while in [20][21][22][23] Ca 2+ diffusion is studied in three-dimensional space. Furthermore, a hybrid model was introduced in [16] to study the Ca 2+ behaviour under the presence of only one channel on the membrane, while in [19] this model has been extended to support tens of channels grouped into clusters. is paper constructs an HPN C model, where the deterministic part depends on the ones proposed in [16, 19-21, 23, 38], while channel subunits are modelled based on the reduction of 8 states to 4 states as proposed in [36,39]. Our model is scalable since it employs colours to denote repeated components (as outlined in the following sections).

Materials and Methods
is section illustrates in more detail the construction of an HPN C model for intracellular calcium dynamics in threedimensional space. e model shows how the calcium dynamics between the ER and the cytosol as well as the corresponding diffusion process can be constructed via HPN C . Figure 3 provides the complete HPN C model. e construction of the HPN C model in Figure 3 can be divided into four parts. In part 1, the deterministic regime is described by specifying calcium reactions, its diffusion and association/dissociation with buffers. In part 2, the stochastic subregime is described as well as the modelling of clusters and their IP 3 R channels with their subunit states. Part 3 connects the stochastic reactions with the continuous ones; vice versa, part 4 connects the continuous regime with the stochastic one. In the following subsections, each part will be discussed in more details.

Modelling 3D Calcium Release.
We first present the deterministic part of the HPN C model in 3D space. is component models the inflow of calcium from ER into cytosol and vice versa. It also models reactions that account for the association/dissociation with buffers. Buffers are used to define the passive and active Ca 2+ as shown in Figure 3, where the place CaER gives the calcium concentration in ER, CaCy gives the calcium concentration in the cytosol, and the places CaCm and CaCs give the calcium concentration in mobile and stationary buffers, respectively. Seven (coloured) continuous transitions represent the interactions between Ca 2+ in cytosol and ER as well as buffers. Calcium inflow from the CaER into CaCy occurs through the transitions Leak and ER2Cy. On the contrary, calcium can return back to the ER from the cytosol by a transition called Pump. e interaction between the buffers and cytosol is controlled by the transitions Kbm − p, Kbm − m, Kbs − p, and Kbs − m. Since our main objective is modelling and simulating cytosolic Ca 2+ inside a cube, we associate all places with the colour set Grid 3D, which is a compound colour set of the dimension xdim × ydim × zdim (see Figure 3 and the declarations given in Table 1). e diffusion process plays a vital role in the calcium dynamics. erefore, we include the diffusion terms modelled by the four transitions Diffc, Diffbm, Diffbs, and DiffE (see Figure 3) in the threedimensional space.
We use a similar approach to the one presented in Figure 1 to model the diffusion dynamics. erefore, the resulting HPN C model is flexible as it permits to adjust the number of grid positions by changing the constants that define the colour set Grid 3D, namely, xdim, ydim, and zdim. In other words, the spatial dimension of the HPN C Ca 2+ model is defined and later adjusted via the constants xdim, ydim, and zdim.    Table 1 for all colourrelated definitions. It is an extension of our model published in [33]. e meaning of arcs and nodes is illustrated in the legend and briefly described in Section 2. e colour code is as follows: dark red: colour set names, red: marking, blue: arc expressions, violet: guards, and black: node names.  //the number of the clusters in x, y and z dir.

Modelling 3D Channel Dynamics.
In this subsection, we explain the modelling of the stochastic part which mainly consists of the dynamics of channel subunits. In our HPN C model, channels can be grouped into clusters. e dynamics of the stochastic part of our HPN C is based on the one given in [36]. In [36], the authors' goal was to construct a simple mathematical model for IP 3 R that can effectively describe the properties of Ca 2+ release events by assuming that every IP 3 R has four equal and independent subunits. Each subunit has three sites of binding: a Ca 2+ activation binding site, an IP 3 binding site, and an inhabitation binding site (section Calcium Release). erefore, every subunit may be found in any of the eight allowed states. ese states control the opening and closing of the corresponding channel, whereby a channel opens if there exist at least three subunits in the active state. us, X ijk denotes the state of a subunit, where i characterises the state of the IP 3 R binding site (i � 0 means not bound), j the state of the Ca 2+ activation binding site, and k the state of inhibiting Ca 2+ binding site. However, according to [36], the states X 010 and X 011 can be neglected (section Calcium Release). erefore the eight states model can be reduced to just six states, and these six are further reduced to four states.
To model the different subunit states, we have created four discrete places, where X 11 denotes a state in which all binding sites are on, while the place X 00 denotes a state in which all binding sites are off. e number of subunits in a given state is described by the number of tokens in the corresponding place. e place X Act determines the number of active subunits for the channel. e transition between different states is naturally of type stochastic. e rate for each transition is defined by the mass action kinetic rate law.
To model a channel in a 3D grid, we use the Ch colour set and the tuple (k, n, j), providing the channel index in 3D. Similarly, for modelling a cluster in a 3D grid, we use the colour set Cl and the tuple (l, m, i) as cluster index. However, channels are scattered in clusters. erefore, we use the compound colour set ChinCl to include channels inside clusters via the variable ((l, m, i), (k, n, j)) that gives the index of a channel inside a cluster as encoded with colours in Figure 3. For instance, the value of the variable of ((2, 2, 2), (2, 1, 3)) refers to a channel located at (2, 1, 3) inside a cluster located at (2, 2, 2). Opening channels depends on their active subunits, where a channel is open if there exist at least three subunits in the active state which is determined by the discrete place X Act .

Modelling the Effect of Channel States on Calcium
Releasing. Releasing of Ca 2+ from ER into cytosol also affects channel opening and closing. erefore, we have constructed a subnet (Figure 3) to account for this behaviour. e place opChans was created with colour set ChinCl to decide the channel state (open/close). e read arc between the place X Act and the immediate transition open with weight three is responsible for switching the channel state from close to open (from 0 to 1), when there are at least three tokens on X Act . e inhibitor arc between the place opChans and the immediate transition open prevents the repeated firing of the transition, as soon as there is a token on the place open. On the contrary, another immediate transition close will switch off the channel state, when the place X Act has less than three tokens. Another inhibitor arc with a weight of three is added for this purpose.
Switching to the deterministic part from the stochastic one is controlled by the continuous transition ER2CT. Firing of the ER2CT transition increases the marking of the continuous place CaC, but in a specific region. e determined grid region, which permits the flow of Ca 2+ from ER into cytosol through the transition ER2CT, is determined in a similar way as in [15,19,20] by the radius of each channel and the number of channels in the open state. is condition is modelled as the rate of the transition ER2CT which is given by [15,19,20].
f(x, y, z, l, m, i) � b(x, y, z, l, m, i) · P ch · (CaER − CaC) , where P ch is a constant, x, y, and z give the grid position in the x, y, and z directions, l, m, and i give the cluster position in the x, y, and z directions, and b(x, y, z, l, m, i) is a Boolean expression having the form: where dx, dy, and dz are the spacing between two grid positions in the x, y, and z directions, CX 0 , CY 0 , and CZ 0 give the position of the first cluster in the x, y, and z directions, Cdx, Cdy, and Cdz give the spacing between two clusters in the x, y, and z directions, and R s is a constant //cluster position in the x, y, and z directions variable k: xdim, n: xdim, j: xdim; //channel position in the x, y, and z directions //functions function IsNeighbour (x, y, z, a, b, c); value representing the radius of a single channel. e variables x, y, z, l, m, and i are specified in the coloured coordinates, while the constants dx, dy, dz, CX 0 , CY 0 , CZ 0 , Cdx, Cdy, and Cdz are specified in the real system coordinates.
e distance between each grid position and the cluster radius is given by Eq. 1. e influx of Ca 2+ from ER into cytosol depends on the closest distance between two grid positions and the cluster radius. Moreover, the cluster radius is not fixed but is determined by the radius of all open channels [19,20].
us, increasing the number of open channels leads to an increasing of the cluster radius, and vice versa. erefore, the coloured stochastic place opChCounter was added with the colour set Cl, for providing the total number of open channels inside every cluster. e transition open can be considered as the controller for the channel opening where tokens are added when the transition fires. In the contrary, the transition close switches off the channel as soon as it fires. Besides, the tokens in the place opChCounter are subsequently used by the deterministic regime to regulate the Ca 2+ flow.

Modelling the Effect of Calcium Releasing on Channel
Transitions. e four transitions of channel subunits are affected by Ca 2+ concentration in the cytosol close to the channel, as discussed in [16,36]. erefore, it is required to sum up all Ca 2+ concentration near the channel and provide the result for the channel transition rates. So, a connection between the deterministic part and the stochastic part is required. In order to fulfil this scenario, we created a coloured continuous place called clCaC to sum up for each channel the nearby Ca 2+ concentrations (Figure 3).
In this paper, we are specifically interested in adjusting the model parameters to determine the number of channels in every cluster as well as the cluster number (demonstrated in the next section). e constants L, M, P determine the number of clusters, while the constants K, N, Q determine the number of channels.

Results and Discussion
is section presents the results of the constructed HPN C model shown in Figure 3 by performing different experiments; it shows how easily we can scale HPN C models by adjusting just a few parameters. Parameters of interest in this sections are L, M, P, K, N, and Q.

Calcium Releasing when a Channel Open.
In this section, we examine Ca 2+ release when a channel is in on state. We run this experiment mainly for two reasons: firstly to check opening and closing of channels based on their number of active subunits and secondly to examine the dependency of calcium releasing when a channel is open. is experiment is configured by setting the constants L � M � P � K � N � Q � 1.
e simulation results are presented in Figures 4 and 5. Figure 4 shows the number of active subunits and the channel states: open/close. Figure 4(a) gives the number of active subunits of the only channel and Figure 4(b) the channel states, both over the whole simulation time. Figure 4(c) gives the number of active subunits in the channel in a small time interval during the simulation time, while Figure 4(d) gives the channel states in the same time interval of the channel subunits. ese results show that the opening of the channel occurs only when the number of active subunits is equal to 3 or more, while closing of the channel occurs when the number of active subunits is less than 3. erefore, these figures confirm that the HPN C model result is valid for the simulation of a single channel in a single cluster. Figure 5 illustrates calcium release from ER to cytosol at different time points. Figure 5(a) presents the calcium concentration inside the cube at the first opening event of the channel, while Figure 5(b) shows the effect of another opening to the channel on the calcium release where Ca 2+ concentration is increased. For further illustration, Figure 5(c) shows the opening of the channel at another time point which clarifies the increasing of the calcium concentration, while Figure 5(d) shows the effect of closing the channel at another time point on Ca 2+ concentration. Please note that although there is no Ca 2+ release due to the channel being closed, there is Ca 2+ transport due to the leak term.

Calcium Diffusion.
e dynamics of intracellular calcium is described by Ca 2+ diffusion. erefore, we have performed two experiments to study the diffusion dynamics of the HPN C model. e first experiment was configured with only one cluster including 24 channels by setting L � M � P � 1, K � 3, N � 4, and Q � 2, while the second experiment was configured with four clusters with every cluster containing 12 channels by setting L � M � 2, P � 1, K � 3, N � 2, and Q � 2. e first experiment shows the effect of the number of channels in the open state on the calcium concentration in the cytosol where increasing the number of open channels also increases the cytosolic calcium concentration. On the contrary, decreasing the number of open channels decreases the cytosolic calcium concentration. e simulation results of this experiment are presented in Figure 6. e number of channels in the open state over the simulation time is presented in Figure 6(a) where this number increases and decreases randomly, while it affects the calcium concentration. For more illustration, consider Figure 6

Validation of the HPN C Model.
In Section 4, three experiments have been presented by adjusting a few parameters only without the need to change the model structure, showing the adaptivity and the flexibility of the HPN C Ca 2+ model. When the model is simulated by using only one cluster, the simulation results agree with the experimental data in [40].
is experiment was implemented to provide evidence for the correctness of the HPN C model by testing the channel opening and closing depending on the number of active subunits. When there are three subunits in the active state, it leads to the opening of the channel. When the number of subunits in the active state is less than three, the channel is switched back to the closed state. e results of active subunits are shown in Figure 4(a) and the corresponding state of the channel is given in Figure 4(b). Besides, the Ca 2+ starts to inflow from the ER membrane into the cytosol after the channel opening. is fact is clarified by considering Figure 5 e second experiment is performed with only one cluster comprising however multiple channels to study the effect of opening multiple channels on the Ca 2+ concentration inside the cytosol as well as comparing the model results with the results of the model presented in [20]. e results of the model in [20] showed that the channel opening and closing occur burst-like and in a stochastic way which is replicated by the HPN C model as shown in Figure 6(a). Moreover, the author in [20] showed the effect of calcium release on the number of open channels. erefore, opening a different number of channels leads to a variation of calcium concentrations in the cytosol at different points of the cube that is showed by the results of HPN C model in Figures 6(b)-6(f ). e results in [22] describe the dynamics of calcium as a wave form, which is also confirmed by the HPN C model. e behaviour observed in [20] is very much similar to the synchronized channel openings in the case of many clusters with respect to oscillation, and the calcium concentration is high only in a narrow range during the release very close to the open channels. Although absolute levels of calcium in the cell stay small, they are enough to devise the synchronized release by cluster-to-cluster coupling. erefore, the behaviour observed in [20] corresponds to the HPN C model results in Figure 7 which is obtained with multiclusters to check the effect of calcium bursting via an open cluster on its neighbours according to the diffusion dynamics. Figure 7 shows the outcome of this experiment. When a cluster opens, other clusters are also highly likely to open. is in turn increases the Ca 2+ concentration inside the cube. A comparison of the behaviour of our HPN C model with results known from similar models in the literature is provided in Table 2.

Simulation Procedure.
e HPN C model in Figure 3 was constructed and simulated by the help of Snoopy [11,41]. Snoopy is a free software that offers a graphical platform to construct and simulate different types of Petri nets including the coloured hybrid Petri nets used in this paper. A complete user manual of Snoopy's HPN C can be found in [42].
Simulation of HPN C in Snoopy is performed by first unfolding the HPN C model into the equivalent HPN one, followed by choosing a suitable hybrid simulation algorithm to execute the unfolded HPN , as presented in [25]. e hybrid simulators of Snoopy simulate the hybrid Petri nets based on the idea introduced in [43] by generating two random numbers. e first random number is used to determine the time of the next stochastic transition to fire, while the second random number is employed to select this transition. Moreover, this procedure is combined with the jump equation that accounts for the evolution of the stochastic transition rates, while the model's deterministic part (represented by a corresponding, automatically generated system of ODEs) is integrated with respect to time.
Unfortunately, the basic idea presented in [43] cannot simulate our HPN C model in reasonable time. erefore, we have added recently more efficient simulation techniques to Snoopy (c.f., [10,30]) which permit to execute larger models in a more efficient way.
In [16,18,19], Ca 2+ hybrid models are simulated by distinguishing between discrete event types according to the different channel transitions. For example, stochastic events that do not impose any effect on the concentration of calcium are not taken into account when calculating the ODE solver step size, while stochastic transitions that affect the  calcium dynamics are taken into consideration when calculating the ODE solver step size. However, the jump equation is not adapted according to this scheme.
Our optimised hybrid algorithms developed in [30,44] can be viewed as a generalisation of the simulation idea of the Ca 2+ model in [16], though they have been independently developed. In other words, we have developed in [30] a general approach to group all reactions of the stochastic regime into independent/dependent reactions based on their connection with the continuous regime. We call the dependent reactions interface reactions. Additionally, we have proposed an approximation of the jump equation where the continuous and stochastic regime become completely independent, but only with regard to the set of interface reactions. In [44], the idea of approximating the jump equation has been optimised by using a similar approach to the one presented in [45]. e runtime evaluation of the newly optimised approach is presented in [44].
In summary, the simulation results presented in this paper are based on the hybrid simulation algorithms introduced in [30,44,46]. e integration time step adaptivity for the underlying ODEs system is achieved by an efficient ODE solver from [47]. e size of the unfolded models as well as the simulation runtime for the three performed experiments are given in Table 3. We adopt a 3D space for     [40], channel opening and closing are based on the activity of subunits where the number of active subunits needs to be at least three for a channel to be opened. In [20], the channel opening and closing occur in a few seconds (e.g., a channel opens approximately every 0.02 s) Figure 4(d) captures a small interval from the simulation time to show the period of channel opening and closing which occurs approximately every 0.02 s.
Opening a channel affects the calcium concentration in cytosol. It increases the value of calcium in cytosol [20]. Figure 5 illustrates the calcium concentration in cytosol at different points of simulation time (e.g., at opening and closing the channel) which shows that the value of calcium concentration increases with channel opening. e results of the model in [20] showed that the opening and closing of the channels occur burst-like and in a stochastic way. In [20] it has been shown that the releasing of calcium affects the number of open channels. us, opening different numbers of channels leads to a variation in the calcium concentrations in the cytosol.
Figures 6(b)-6(f ) describe the dynamics of calcium in cytosol inside a cube, where the calcium diffuses to the adjacent grid points, which increases the chance of channel opening at these points. e results in [22] describe the dynamics of calcium as a wave form. e wave-like shape can be clearly seen in Figure 6.
e behaviour observation described in [20] is very similar to the openings of the synchronized channels in many oscillating clusters, and the calcium concentration is abundant only in a small range during Ca 2+ release due to channel opening. Figure 7 presents the results of the experiment with multiclusters inside the cube which shows opening more than one cluster at the same time where the value of calcium is high at the center of the cluster (narrow range), while low calcium concentrations lead to the opening of other clusters.
this study in a cube form, by uniformly setting D to 5 for all experiments, which corresponds to 9μ × 9μ × 9μ.

Conclusions
is paper has presented the construction and simulation of 3D models using HPN C . To illustrate our idea, we have applied this approach to construct a graphical model for intracellular calcium dynamics in the three-dimensional space exploiting the power of coloured hybrid Petri nets. We explained in more details how the HPN C model for intracellular calcium dynamics is created. e model is divided into two main parts: the continuous and the stochastic part. e continuous part describes the deterministic calcium release between cytosol, ER, and buffers as well as the diffusion process. e stochastic part models the clusters with the channel subunits. e modelling of the continuous and stochastic parts with coloured continuous and stochastic Petri nets classes is thoroughly discussed as well as the switching from the stochastic part to the deterministic part and vice versa. e flexibility of HPN C is showed by changing only a few parameters to run more than one experiment. e results of the developed HPN C have been validated against results known from the literature.
Our HPN C model can be further developed to study the behaviour of calcium under other types of channels such as RyR, where the opening occurs with high probability, even when the calcium concentration is low [22]. Furthermore, we will apply our generic 3D model construction approach using HPN C to other modelling scenarios in future work.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare no conflicts of interest regarding the publication of this paper.