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

Authors: Amr Ismail; Mostafa Herajy; Elsayed Atlam; Monika Heiner
Journal: Modelling and Simulation in Engineering
 (2020)

---

## Abstract

Three-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. The 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 (HPNC) 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, HPNC models are also scalable. This paper shows how HPNC 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 HPNC can yield more flexible models as the structure can be easily scaled up and down by just modifying a few parameters. This advantage of convenient model configuration facilitates the design of different experiments without the need to alter the model structure.


---

## Body

## 1. 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–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 (PNC) [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 onHPN and PNC, coloured Hybrid Petri Nets HPNC have been proposed in [8–10] combining all the merits of HPN and PNC into one net class. The 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, PNC can be converted algorithmically into equivalent low-level PN via an unfolding process [11]. Indeed, HPNC 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 (Ca2+) dynamics [12–15] is an important biological pathway that has been recently widely studied. Intracellular Ca2+ 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 Ca2+ dynamics. So far, many models were introduced to analyse Ca2+ dynamics in two-dimensional space (see, e.g., [14–19]), but only a few models for the three-dimensional space (see, e.g., [20–23]). However, constructing an in silico Ca2+ model cannot be achieved by using a single modelling paradigm (i.e., deterministic or stochastic).On one hand, theCa2+ 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. These 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.Thus, 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.This paper utilisesHPNC [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 Ca2+ behaviour. The 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). The constructed model is validated by conducting a few experiments and comparing their outcome with facts known from the literature.The 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 Ca2+ dynamics as well as the utilisation of HPNC for modelling biological processes. Section 4 shows how HPNC can be used to construct and execute biological models illustrated by the Ca2+ case study. Simulation results of the Ca2+ model are presented in Section 4. We conclude the paper by some closing remarks and an outline of possible future work.


## 2. 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.


### 2.1. Petri Nets

Petri Nets (PN) are an excellent modelling formalism for studying the behaviour of complex systems made of concurrent components. They 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.Many extensions based on the standardPN were proposed, such as Stochastic Petri Nets (SPN) [4], Continuous Petri Nets (CPN) [5], and Hybrid Petri Nets (HPN) [1, 25].SPN extend PN by associating each transition with a random firing rate. The 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]. The 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].The motivation ofGHPN, 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. Their 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. The 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.


### 2.2. Coloured Hybrid Petri Nets

Coloured Petri Nets (PNC) extend low-level PN with a set of finite colour sets [6, 7]. PNC offer a parameterised approach to construct large biological systems; there are many coloured net classes, among them coloured continuous Petri nets (CPNC) [7] and coloured stochastic Petri nets (SPNC) [31].By augmentingGHPN with colours, we obtain a new class of hybrid Petri nets called Coloured Generalised Hybrid Petri Nets (GHPNC) [8]. In GHPNC, 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.GHPNC 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 CPNC to model diffusion in 3D space.Figure 1
An illustration of diffusion in 3D space. (a) Six neighbours for a species at the indexl,m,k in the form: l±1,m,k,l,m±1,k or l,m,z±1. (b) The 26 neighbours due to one of their faces, corners, or edges: 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. (c) HPNC implementation of 3D diffusion. (d) The corresponding colour definitions of the HPNC in (c). IsNeighbourx,y,z,a,b,c is a Boolean function specifying the neighbourhood relation of any two components in 3D space according to (a) or (b).
(a)(b)(c)(d)


### 2.3. Colour Representation of 3D Space

When constructing 3D models, the components of the model can be arranged in a cube form as shown in Figures1(a) and 1(b). Figure 1(a) illustrates a component at index l,m,k and its six neighbours at different indices. These 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. These 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). HPNC 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 theHPNC, but also employs their colour expressions. To present a component we use the graph structure of HPNC 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 Figure1, 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 two-dimensional space [33], we define, in this paper, a three-dimensional compound colour set, Grid3D, 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 IsNeighbourx,y,z,a,b,c:bool IsNeighbour (xdim x, ydim y, zdim z, xdim a, ydim b, zdim c){(a = x + 1 & b = y & c = z)//neighbouring conditions| (a =x − 1 & b = y & c = z)| (a = x & b = y + 1 & c = z)| (a = x & b = y − 1 & c = z)| (a = x & b = y & c = z + 1)| (a = x & b = y & c = z − 1)& (1 ≤a & a ≤ D1) & (1 ≤ b & b ≤ D2) & (1 ≤ c & c ≤ D3) //boundary conditions}where xdim, ydim, and zdim are the same as in Figure 1(d). The function IsNeighbour () returns either true or false depending on the evaluation of the function body. The 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 developedHPNC model, we make use of the component presented in Figure 1 as a building block for modelling calcium diffusion in 3D. This 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 ofGHPNC 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 GHPNC permit to combine different timescales.


### 2.4. Intracellular Calcium Dynamics

This section provides an example of howGHPNC can be employed to construct and simulate 3D models. The case study is about the 3D modelling of intracellular calcium dynamics. This specific example is selected because it is intricate to model it via existing modelling approaches. Moreover, this case study is 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 Ca2+ dynamics, followed by a detailed discussion of the GHPNC model construction.


### 2.5. Calcium Release

Intracellular calcium is organised by the release ofCa2+ ions from the cell membrane or internal stores such as Endoplasmic Reticulum (ER). Releasing of Ca2+ from ER occurs via the opening of inositol-4,5-triphosphate receptor (IP3R) channels (among different types, such as RyR and L-type calcium channels [13]). These channels are randomly grouped into clusters [13]. A stochastic model to study the random opening and closing of IP3R channels is presented in [17]. The 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 Ca2+ release from the ER to cytosol.Figure 2
A schematic diagram describing the different processes and components involved in theCa2+ model according to [19, 20, 38]. When some channels are open, Ca2+ flows from the ER to the cytosol. Ca2+ returns back to the ER via a process called pump. Channel opening is affected by the Ca2+ concentration surrounding the channel. Ca2+ in cytosol can be associated/dissociated with different types of buffer. Ca2+ in cytosol, ER, and Buffers rapidly diffuses to the neighbouring cells.As it is detailed in [17], (IP3R) channels consist of four identical and independent subunits which can be combined with the deterministic part to study the Ca2+ 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 IP3, one for activating Ca2+, and a third for Ca2+ inhibiting. Each of these sites is either busy or idle by their particular regulatory agonist (Ca2+ or IP3) [17, 34]. Building on the experimental studies in [35], Marchant et al. observed that in the absence of bound Ca2+, IP3 does not bind to its activating site. Therefore, the two states of unbound IP3 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 IP3 (the fastest), binding by Ca2+ activation, and binding by Ca2+ inhibition (the slowest). As IP3 dissociation and association is faster than the competing processes, transitions between states of unbound IP3 and bound IP3 are the first to reach the equilibrium state. Therefore, the unbound IP3 state is grouped with the bound IP3 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.The 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 5th state.Finally, the opening of anIP3R channel is influenced by the neighbouring Ca2+ concentration in the cytosol. The 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. The 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 Ca2+ dynamics in two-dimensional space, while in [20–23] Ca2+ diffusion is studied in three-dimensional space. Furthermore, a hybrid model was introduced in [16] to study the Ca2+ 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.This paper constructs anHPNC 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).


## 2.1. Petri Nets

Petri Nets (PN) are an excellent modelling formalism for studying the behaviour of complex systems made of concurrent components. They 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.Many extensions based on the standardPN were proposed, such as Stochastic Petri Nets (SPN) [4], Continuous Petri Nets (CPN) [5], and Hybrid Petri Nets (HPN) [1, 25].SPN extend PN by associating each transition with a random firing rate. The 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]. The 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].The motivation ofGHPN, 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. Their 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. The 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.


## 2.2. Coloured Hybrid Petri Nets

Coloured Petri Nets (PNC) extend low-level PN with a set of finite colour sets [6, 7]. PNC offer a parameterised approach to construct large biological systems; there are many coloured net classes, among them coloured continuous Petri nets (CPNC) [7] and coloured stochastic Petri nets (SPNC) [31].By augmentingGHPN with colours, we obtain a new class of hybrid Petri nets called Coloured Generalised Hybrid Petri Nets (GHPNC) [8]. In GHPNC, 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.GHPNC 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 CPNC to model diffusion in 3D space.Figure 1
An illustration of diffusion in 3D space. (a) Six neighbours for a species at the indexl,m,k in the form: l±1,m,k,l,m±1,k or l,m,z±1. (b) The 26 neighbours due to one of their faces, corners, or edges: 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. (c) HPNC implementation of 3D diffusion. (d) The corresponding colour definitions of the HPNC in (c). IsNeighbourx,y,z,a,b,c is a Boolean function specifying the neighbourhood relation of any two components in 3D space according to (a) or (b).
(a)(b)(c)(d)


## 2.3. Colour Representation of 3D Space

When constructing 3D models, the components of the model can be arranged in a cube form as shown in Figures1(a) and 1(b). Figure 1(a) illustrates a component at index l,m,k and its six neighbours at different indices. These 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. These 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). HPNC 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 theHPNC, but also employs their colour expressions. To present a component we use the graph structure of HPNC 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 Figure1, 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 two-dimensional space [33], we define, in this paper, a three-dimensional compound colour set, Grid3D, 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 IsNeighbourx,y,z,a,b,c:bool IsNeighbour (xdim x, ydim y, zdim z, xdim a, ydim b, zdim c){(a = x + 1 & b = y & c = z)//neighbouring conditions| (a =x − 1 & b = y & c = z)| (a = x & b = y + 1 & c = z)| (a = x & b = y − 1 & c = z)| (a = x & b = y & c = z + 1)| (a = x & b = y & c = z − 1)& (1 ≤a & a ≤ D1) & (1 ≤ b & b ≤ D2) & (1 ≤ c & c ≤ D3) //boundary conditions}where xdim, ydim, and zdim are the same as in Figure 1(d). The function IsNeighbour () returns either true or false depending on the evaluation of the function body. The 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 developedHPNC model, we make use of the component presented in Figure 1 as a building block for modelling calcium diffusion in 3D. This 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 ofGHPNC 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 GHPNC permit to combine different timescales.


## 2.4. Intracellular Calcium Dynamics

This section provides an example of howGHPNC can be employed to construct and simulate 3D models. The case study is about the 3D modelling of intracellular calcium dynamics. This specific example is selected because it is intricate to model it via existing modelling approaches. Moreover, this case study is 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 Ca2+ dynamics, followed by a detailed discussion of the GHPNC model construction.


## 2.5. Calcium Release

Intracellular calcium is organised by the release ofCa2+ ions from the cell membrane or internal stores such as Endoplasmic Reticulum (ER). Releasing of Ca2+ from ER occurs via the opening of inositol-4,5-triphosphate receptor (IP3R) channels (among different types, such as RyR and L-type calcium channels [13]). These channels are randomly grouped into clusters [13]. A stochastic model to study the random opening and closing of IP3R channels is presented in [17]. The 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 Ca2+ release from the ER to cytosol.Figure 2
A schematic diagram describing the different processes and components involved in theCa2+ model according to [19, 20, 38]. When some channels are open, Ca2+ flows from the ER to the cytosol. Ca2+ returns back to the ER via a process called pump. Channel opening is affected by the Ca2+ concentration surrounding the channel. Ca2+ in cytosol can be associated/dissociated with different types of buffer. Ca2+ in cytosol, ER, and Buffers rapidly diffuses to the neighbouring cells.As it is detailed in [17], (IP3R) channels consist of four identical and independent subunits which can be combined with the deterministic part to study the Ca2+ 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 IP3, one for activating Ca2+, and a third for Ca2+ inhibiting. Each of these sites is either busy or idle by their particular regulatory agonist (Ca2+ or IP3) [17, 34]. Building on the experimental studies in [35], Marchant et al. observed that in the absence of bound Ca2+, IP3 does not bind to its activating site. Therefore, the two states of unbound IP3 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 IP3 (the fastest), binding by Ca2+ activation, and binding by Ca2+ inhibition (the slowest). As IP3 dissociation and association is faster than the competing processes, transitions between states of unbound IP3 and bound IP3 are the first to reach the equilibrium state. Therefore, the unbound IP3 state is grouped with the bound IP3 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.The 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 5th state.Finally, the opening of anIP3R channel is influenced by the neighbouring Ca2+ concentration in the cytosol. The 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. The 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 Ca2+ dynamics in two-dimensional space, while in [20–23] Ca2+ diffusion is studied in three-dimensional space. Furthermore, a hybrid model was introduced in [16] to study the Ca2+ 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.This paper constructs anHPNC 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).


## 3. Materials and Methods

This section illustrates in more detail the construction of anHPNC model for intracellular calcium dynamics in three-dimensional space. The model shows how the calcium dynamics between the ER and the cytosol as well as the corresponding diffusion process can be constructed via HPNC. Figure 3 provides the complete HPNC model.Figure 3
TheHPNC model of Intracellular Calcium in 3D. Stochastic channel processes are represented by the part on the right, while the continuous processes of the calcium dynamics are represented by the left part. The upper part at the center of the model represents the connection from the stochastic part to the continuous part, vice versa for the lower part at the center of the model. See Table 1 for all colour-related definitions. It is an extension of our model published in [33]. The meaning of arcs and nodes is illustrated in the legend and briefly described in Section2. The colour code is as follows: dark red: colour set names, red: marking, blue: arc expressions, violet: guards, and black: node names.Table 1
Colour declarations of theHPNC in Figure 3.
Declarations//constants:constantD = 5;//the uniform size of all dimensionsconstantD1 = D;//the size ofx-dimensionconstantD2 = D;//the size ofy-dimensionconstantD3 = D;//the size ofz-dimensionconstantMID = D/2;//the middle point of any dimensionconstantclX0 = MID;//the index of the first cluster inx-dir.constantclY0 = MID;//the index of the first cluster iny-dir.constantclZ0 = MID;//the index of the first cluster in z dir.constantclDx = 1, clDy = 1, clDz = 1;//the distance between two clustersconstantL = 1, M = 1, P = 1;//the number of the clusters inx, y and z dir.constantK = 1, N = 1, Q = 1//the number of the channel inx, y and z dir.//colour sets:coloursetxdim = int with 1 to D1;coloursetydim = int with 1 to D2;coloursetzdim = int with 1 to D3;coloursetGrid3D = product (xdim, ydim, zdim);coloursetrow = int with 1 to L;coloursetcolumn = int with 1 to M;coloursetdim = int with 1 to P;coloursetcrow = int with 1 to K;coloursetccolumn = int with 1 to N;coloursetcdim = int with 1 to Q;coloursetCh = product (crow, ccolumn, cdim);//C incrow, ccolumn, cdim means channelcoloursetCl = product (row, column, dim);coloursetChinCl = product (clus, chans);//variables:variablex:xdim, y:ydim, z:zdim;//grid positions in thex, y, and z directionsvariablea:xdim, b:ydim, c:zdim;//grid positions in thex, y, and z directionsvariablel:xdim, m:ydim, i:zdim;//cluster position in thex, y, and z directionsvariablek:xdim, n:xdim, j:xdim;//channel position in thex, y, and z directions//functionsfunction IsNeighbour (x, y, z, a, b, c);The construction of theHPNC 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 IP3R 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.


### 3.1. Modelling 3D Calcium Release

We first present the deterministic part of theHPNC model in 3D space. This 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 Ca2+ 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 betweenCa2+ 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. The 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 Ca2+ inside a cube, we associate all places with the colour set Grid3D, which is a compound colour set of the dimension xdim×ydim×zdim (see Figure 3 and the declarations given in Table 1). The diffusion process plays a vital role in the calcium dynamics. Therefore, we include the diffusion terms modelled by the four transitions Diffc, Diffbm, Diffbs, and DiffE (see Figure 3) in the three-dimensional space.We use a similar approach to the one presented in Figure1 to model the diffusion dynamics. Therefore, the resulting HPNC model is flexible as it permits to adjust the number of grid positions by changing the constants that define the colour set Grid3D, namely, xdim, ydim, and zdim. In other words, the spatial dimension of the HPNCCa2+ model is defined and later adjusted via the constants xdim, ydim, and zdim.


### 3.2. 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 ourHPNC model, channels can be grouped into clusters.The dynamics of the stochastic part of ourHPNC is based on the one given in [36]. In [36], the authors’ goal was to construct a simple mathematical model for IP3R that can effectively describe the properties of Ca2+ release events by assuming that every IP3R has four equal and independent subunits. Each subunit has three sites of binding: a Ca2+ activation binding site, an IP3 binding site, and an inhabitation binding site (section Calcium Release). Therefore, every subunit may be found in any of the eight allowed states. These 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. Thus, Xijk denotes the state of a subunit, where i characterises the state of the IP3R binding site (i=0 means not bound), j the state of the Ca2+ activation binding site, and k the state of inhibiting Ca2+ binding site. However, according to [36], the states X010 and X011 can be neglected (section Calcium Release). Therefore 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, whereX11 denotes a state in which all binding sites are on, while the place X00 denotes a state in which all binding sites are off. The number of subunits in a given state is described by the number of tokens in the corresponding place. The place XAct determines the number of active subunits for the channel. The transition between different states is naturally of type stochastic. The rate for each transition is defined by the mass action kinetic rate law.To model a channel in a 3D grid, we use theCh 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. Therefore, 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 XAct.


### 3.3. Modelling the Effect of Channel States on Calcium Releasing

Releasing ofCa2+ from ER into cytosol also affects channel opening and closing. Therefore, we have constructed a subnet (Figure 3) to account for this behaviour. The place opChans was created with colour set ChinCl to decide the channel state (open/close). The read arc between the place XAct 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 XAct. The 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 XAct 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 transitionER2CT. Firing of the ER2CT transition increases the marking of the continuous place CaC, but in a specific region. The determined grid region, which permits the flow of Ca2+ 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. This condition is modelled as the rate of the transition ER2CT which is given by [15, 19, 20].(1)fx,y,z,l,m,i=bx,y,z,l,m,i⋅Pch⋅CaER−CaC,where Pch 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 bx,y,z,l,m,i is a Boolean expression having the form:(2)bx,y,z,l,m,i=d12+d22+d32<Rs⋅opChCounter,d1=x⋅dx−CX0+l−1⋅cdx,d2=y⋅dy−CY0+m−1⋅cdy,d3=z⋅dz−CZ0+i−1⋅cdz,where dx, dy, and dz are the spacing between two grid positions in the x, y, and z directions, CX0, CY0, and CZ0 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 Rs is a constant value representing the radius of a single channel. The variables x, y, z, l, m, and i are specified in the coloured coordinates, while the constants dx, dy, dz, CX0, CY0, CZ0, Cdx, Cdy, and Cdz are specified in the real system coordinates.The distance between each grid position and the cluster radius is given by Eq.1. The influx of Ca2+ 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]. Thus, increasing the number of open channels leads to an increasing of the cluster radius, and vice versa.Therefore, the coloured stochastic placeopChCounter was added with the colour set Cl, for providing the total number of open channels inside every cluster. The 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 Ca2+ flow.


### 3.4. Modelling the Effect of Calcium Releasing on Channel Transitions

The four transitions of channel subunits are affected byCa2+ concentration in the cytosol close to the channel, as discussed in [16, 36]. Therefore, it is required to sum up all Ca2+ 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 Ca2+ 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). The constantsL, M, P determine the number of clusters, while the constants K, N, Q determine the number of channels.


## 3.1. Modelling 3D Calcium Release

We first present the deterministic part of theHPNC model in 3D space. This 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 Ca2+ 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 betweenCa2+ 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. The 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 Ca2+ inside a cube, we associate all places with the colour set Grid3D, which is a compound colour set of the dimension xdim×ydim×zdim (see Figure 3 and the declarations given in Table 1). The diffusion process plays a vital role in the calcium dynamics. Therefore, we include the diffusion terms modelled by the four transitions Diffc, Diffbm, Diffbs, and DiffE (see Figure 3) in the three-dimensional space.We use a similar approach to the one presented in Figure1 to model the diffusion dynamics. Therefore, the resulting HPNC model is flexible as it permits to adjust the number of grid positions by changing the constants that define the colour set Grid3D, namely, xdim, ydim, and zdim. In other words, the spatial dimension of the HPNCCa2+ model is defined and later adjusted via the constants xdim, ydim, and zdim.


## 3.2. 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 ourHPNC model, channels can be grouped into clusters.The dynamics of the stochastic part of ourHPNC is based on the one given in [36]. In [36], the authors’ goal was to construct a simple mathematical model for IP3R that can effectively describe the properties of Ca2+ release events by assuming that every IP3R has four equal and independent subunits. Each subunit has three sites of binding: a Ca2+ activation binding site, an IP3 binding site, and an inhabitation binding site (section Calcium Release). Therefore, every subunit may be found in any of the eight allowed states. These 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. Thus, Xijk denotes the state of a subunit, where i characterises the state of the IP3R binding site (i=0 means not bound), j the state of the Ca2+ activation binding site, and k the state of inhibiting Ca2+ binding site. However, according to [36], the states X010 and X011 can be neglected (section Calcium Release). Therefore 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, whereX11 denotes a state in which all binding sites are on, while the place X00 denotes a state in which all binding sites are off. The number of subunits in a given state is described by the number of tokens in the corresponding place. The place XAct determines the number of active subunits for the channel. The transition between different states is naturally of type stochastic. The rate for each transition is defined by the mass action kinetic rate law.To model a channel in a 3D grid, we use theCh 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. Therefore, 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 XAct.


## 3.3. Modelling the Effect of Channel States on Calcium Releasing

Releasing ofCa2+ from ER into cytosol also affects channel opening and closing. Therefore, we have constructed a subnet (Figure 3) to account for this behaviour. The place opChans was created with colour set ChinCl to decide the channel state (open/close). The read arc between the place XAct 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 XAct. The 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 XAct 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 transitionER2CT. Firing of the ER2CT transition increases the marking of the continuous place CaC, but in a specific region. The determined grid region, which permits the flow of Ca2+ 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. This condition is modelled as the rate of the transition ER2CT which is given by [15, 19, 20].(1)fx,y,z,l,m,i=bx,y,z,l,m,i⋅Pch⋅CaER−CaC,where Pch 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 bx,y,z,l,m,i is a Boolean expression having the form:(2)bx,y,z,l,m,i=d12+d22+d32<Rs⋅opChCounter,d1=x⋅dx−CX0+l−1⋅cdx,d2=y⋅dy−CY0+m−1⋅cdy,d3=z⋅dz−CZ0+i−1⋅cdz,where dx, dy, and dz are the spacing between two grid positions in the x, y, and z directions, CX0, CY0, and CZ0 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 Rs is a constant value representing the radius of a single channel. The variables x, y, z, l, m, and i are specified in the coloured coordinates, while the constants dx, dy, dz, CX0, CY0, CZ0, Cdx, Cdy, and Cdz are specified in the real system coordinates.The distance between each grid position and the cluster radius is given by Eq.1. The influx of Ca2+ 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]. Thus, increasing the number of open channels leads to an increasing of the cluster radius, and vice versa.Therefore, the coloured stochastic placeopChCounter was added with the colour set Cl, for providing the total number of open channels inside every cluster. The 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 Ca2+ flow.


## 3.4. Modelling the Effect of Calcium Releasing on Channel Transitions

The four transitions of channel subunits are affected byCa2+ concentration in the cytosol close to the channel, as discussed in [16, 36]. Therefore, it is required to sum up all Ca2+ 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 Ca2+ 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). The constantsL, M, P determine the number of clusters, while the constants K, N, Q determine the number of channels.


## 4. Results and Discussion

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


### 4.1. Calcium Releasing when a Channel Open

In this section, we examineCa2+ 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. This experiment is configured by setting the constants L=M=P=K=N=Q=1.The simulation results are presented in Figures4 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. These 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. Therefore, these figures confirm that the HPNC model result is valid for the simulation of a single channel in a single cluster.Figure 4
Simulation results of theHPNC model given in Figure 3 configured with a single cluster comprising only one channel. (a) Number of subunits in the active state over the entire simulation time. (b) Channel state (open/close) over the entire simulation time. (c) Number of subunits in the active state in a small time interval of the simulation time. (d) Channel state (open/close) of the subunit in the same time interval of the simulation time.
(a)(b)(c)(d)Figure 5
Simulation results of theHPNC model given in Figure 3 with a single cluster comprising one channel. (a) Calcium concentration in cytosol when the only channel opens (for the first time) at time = 7.42 seconds, (b) calcium concentration in cytosol at another opening of the channel at time = 7.45 seconds, (c) calcium concentration in cytosol at another opening of the channel at time = 7.58 seconds, and (d) calcium concentration in cytosol at closing the channel at time = 7.98 seconds.
(a)(b)(c)(d)Figure5 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 Ca2+ 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 Ca2+ concentration. Please note that although there is no Ca2+ release due to the channel being closed, there is Ca2+ transport due to the leak term.


### 4.2. Calcium Diffusion

The dynamics of intracellular calcium is described byCa2+ diffusion. Therefore, we have performed two experiments to study the diffusion dynamics of the HPNC model. The 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.The 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. The simulation results of this experiment are presented in Figure6. The 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(b) which represents the Ca2+ concentration inside the grid (cube form) when six channels are in the open state, while Figure 6(c) shows the same results of the Ca2+ concentration when seven channels are in the open state. Figure 6(d) shows Ca2+ concentration for a few number of open channels. At time = 6.01 seconds, there are no channels in the open state; therefore, the value of Ca2+ concentration dropped down in Figure 6(e), while there are two channels in the open state at time = 6.98 seconds. Therefore, the value of Ca2+ starts increasing again (Figure 6(f)).Figure 6
Simulation results of one cluster containing 24 channels. (a) The number of open channels over the time, (b–f) calcium concentration in cytosol at different time points.
(a)(b)(c)(d)(e)(f)The second experiment shows the diffusion of cytosolic calcium in the cube with three and four clusters, respectively. Clusters are adjusted at different points in the grid where opening of one cluster affects its neighbours. The results of this experiment are presented in Figure7. The effect of opening three clusters is shown in Figures 7(a) and 7(b), where Figure 7(a) presents the density of calcium in the cytosol and Figure 7(b) visualizes the density of calcium in a slice view at the same time point as in Figure 7(a). Figure 7(c) shows the diffusion of calcium inside the grid which leads to the opening of another cluster and the corresponding slice view is provided in Figure 7(d). It is worth mentioning that the slice views of the calcium diffusion in Figures 7(b) and 7(d) provide a clearer visualization of the number of open clusters that what is given in Figures 7(a) and 7(c).Figure 7
Simulation results of theCa2+ model in Figure 3 with more than one cluster, every one contains 12 channels. (a) With three clusters, (b) with three clusters (slice plot), (c) with four clusters, and (d) with four clusters (slice plot).
(a)(b)(c)(d)


### 4.3. Validation of theHPNC Model

In Section4, 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 HPNCCa2+ model. When the model is simulated by using only one cluster, the simulation results agree with the experimental data in [40]. This experiment was implemented to provide evidence for the correctness of the HPNC 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. The 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 Ca2+ starts to inflow from the ER membrane into the cytosol after the channel opening. This fact is clarified by considering Figure 5(c), showing when a channel is in the open state. On the contrary, Ca2+ cannot inflow when the channel is in the closed state; compare Figure 5(d).The second experiment is performed with only one cluster comprising however multiple channels to study the effect of opening multiple channels on theCa2+ concentration inside the cytosol as well as comparing the model results with the results of the model presented in [20]. The 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 HPNC model as shown in Figure 6(a). Moreover, the author in [20] showed the effect of calcium release on the number of open channels. Therefore, 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 HPNC model in Figures 6(b)–6(f). The results in [22] describe the dynamics of calcium as a wave form, which is also confirmed by the HPNC model.The 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. Therefore, the behaviour observed in [20] corresponds to the HPNC 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. This in turn increases the Ca2+ concentration inside the cube. A comparison of the behaviour of our HPNC model with results known from similar models in the literature is provided in Table 2.Table 2
ComparingHPNC model results with results from the literature.
KnownCa2+ behaviour (literature)HPNC model BehaviourIn [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.Figure4(a) shows the number of active subunits and Figure 4(b) the corresponding open channels, when the number of active subunits is more than three.In [20], the channel opening and closing occur in a few seconds (e.g., a channel opens approximately every 0.02 s)Figure4(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].Figure5 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.The results of the model in [20] showed that the opening and closing of the channels occur burst-like and in a stochastic way.Figure6(a) depicts the number of open channels inside a cluster which exhibits a similar burst and happens in a stochastic way.In [20] it has been shown that the releasing of calcium affects the number of open channels. Thus, opening different numbers of channels leads to a variation in the calcium concentrations in the cytosol.Figures6(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.The results in [22] describe the dynamics of calcium as a wave form.The wave-like shape can be clearly seen in Figure6.The 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 Ca2+ release due to channel opening.Figure7 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.


### 4.4. Simulation Procedure

TheHPNC 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 HPNC can be found in [42].Simulation ofHPNC in Snoopy is performed by first unfolding the HPNC model into the equivalent HPN one, followed by choosing a suitable hybrid simulation algorithm to execute the unfolded HPN, as presented in [25]. The hybrid simulators of Snoopy simulate the hybrid Petri nets based on the idea introduced in [43] by generating two random numbers. The 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 HPNC model in reasonable time. Therefore, 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], Ca2+ 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 Ca2+ 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]. The 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]. The integration time step adaptivity for the underlying ODEs system is achieved by an efficient ODE solver from [47]. The 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 this study in a cube form, by uniformly setting D to 5 for all experiments, which corresponds to 9μ × 9μ × 9μ.Table 3
Key measures for the three simulation experiments.
#CustersChannelsUnfolded placesUnfolded transitionsRun time∗ (minutes)1L=1,M=1,P=1K = 1, N = 1, Q = 12,92429,52470.082L=1,M=1,P=1K = 3, N = 4, Q = 23,06229,800260.533L=2,M=2,P=1K = 3, N = 2, Q = 23,21234,462497.73∗Done using OS X, 2.6 GHz Intel Core i7, 16 GB 1600 MHz DDR3.


## 4.1. Calcium Releasing when a Channel Open

In this section, we examineCa2+ 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. This experiment is configured by setting the constants L=M=P=K=N=Q=1.The simulation results are presented in Figures4 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. These 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. Therefore, these figures confirm that the HPNC model result is valid for the simulation of a single channel in a single cluster.Figure 4
Simulation results of theHPNC model given in Figure 3 configured with a single cluster comprising only one channel. (a) Number of subunits in the active state over the entire simulation time. (b) Channel state (open/close) over the entire simulation time. (c) Number of subunits in the active state in a small time interval of the simulation time. (d) Channel state (open/close) of the subunit in the same time interval of the simulation time.
(a)(b)(c)(d)Figure 5
Simulation results of theHPNC model given in Figure 3 with a single cluster comprising one channel. (a) Calcium concentration in cytosol when the only channel opens (for the first time) at time = 7.42 seconds, (b) calcium concentration in cytosol at another opening of the channel at time = 7.45 seconds, (c) calcium concentration in cytosol at another opening of the channel at time = 7.58 seconds, and (d) calcium concentration in cytosol at closing the channel at time = 7.98 seconds.
(a)(b)(c)(d)Figure5 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 Ca2+ 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 Ca2+ concentration. Please note that although there is no Ca2+ release due to the channel being closed, there is Ca2+ transport due to the leak term.


## 4.2. Calcium Diffusion

The dynamics of intracellular calcium is described byCa2+ diffusion. Therefore, we have performed two experiments to study the diffusion dynamics of the HPNC model. The 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.The 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. The simulation results of this experiment are presented in Figure6. The 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(b) which represents the Ca2+ concentration inside the grid (cube form) when six channels are in the open state, while Figure 6(c) shows the same results of the Ca2+ concentration when seven channels are in the open state. Figure 6(d) shows Ca2+ concentration for a few number of open channels. At time = 6.01 seconds, there are no channels in the open state; therefore, the value of Ca2+ concentration dropped down in Figure 6(e), while there are two channels in the open state at time = 6.98 seconds. Therefore, the value of Ca2+ starts increasing again (Figure 6(f)).Figure 6
Simulation results of one cluster containing 24 channels. (a) The number of open channels over the time, (b–f) calcium concentration in cytosol at different time points.
(a)(b)(c)(d)(e)(f)The second experiment shows the diffusion of cytosolic calcium in the cube with three and four clusters, respectively. Clusters are adjusted at different points in the grid where opening of one cluster affects its neighbours. The results of this experiment are presented in Figure7. The effect of opening three clusters is shown in Figures 7(a) and 7(b), where Figure 7(a) presents the density of calcium in the cytosol and Figure 7(b) visualizes the density of calcium in a slice view at the same time point as in Figure 7(a). Figure 7(c) shows the diffusion of calcium inside the grid which leads to the opening of another cluster and the corresponding slice view is provided in Figure 7(d). It is worth mentioning that the slice views of the calcium diffusion in Figures 7(b) and 7(d) provide a clearer visualization of the number of open clusters that what is given in Figures 7(a) and 7(c).Figure 7
Simulation results of theCa2+ model in Figure 3 with more than one cluster, every one contains 12 channels. (a) With three clusters, (b) with three clusters (slice plot), (c) with four clusters, and (d) with four clusters (slice plot).
(a)(b)(c)(d)


## 4.3. Validation of theHPNC Model

In Section4, 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 HPNCCa2+ model. When the model is simulated by using only one cluster, the simulation results agree with the experimental data in [40]. This experiment was implemented to provide evidence for the correctness of the HPNC 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. The 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 Ca2+ starts to inflow from the ER membrane into the cytosol after the channel opening. This fact is clarified by considering Figure 5(c), showing when a channel is in the open state. On the contrary, Ca2+ cannot inflow when the channel is in the closed state; compare Figure 5(d).The second experiment is performed with only one cluster comprising however multiple channels to study the effect of opening multiple channels on theCa2+ concentration inside the cytosol as well as comparing the model results with the results of the model presented in [20]. The 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 HPNC model as shown in Figure 6(a). Moreover, the author in [20] showed the effect of calcium release on the number of open channels. Therefore, 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 HPNC model in Figures 6(b)–6(f). The results in [22] describe the dynamics of calcium as a wave form, which is also confirmed by the HPNC model.The 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. Therefore, the behaviour observed in [20] corresponds to the HPNC 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. This in turn increases the Ca2+ concentration inside the cube. A comparison of the behaviour of our HPNC model with results known from similar models in the literature is provided in Table 2.Table 2
ComparingHPNC model results with results from the literature.
KnownCa2+ behaviour (literature)HPNC model BehaviourIn [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.Figure4(a) shows the number of active subunits and Figure 4(b) the corresponding open channels, when the number of active subunits is more than three.In [20], the channel opening and closing occur in a few seconds (e.g., a channel opens approximately every 0.02 s)Figure4(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].Figure5 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.The results of the model in [20] showed that the opening and closing of the channels occur burst-like and in a stochastic way.Figure6(a) depicts the number of open channels inside a cluster which exhibits a similar burst and happens in a stochastic way.In [20] it has been shown that the releasing of calcium affects the number of open channels. Thus, opening different numbers of channels leads to a variation in the calcium concentrations in the cytosol.Figures6(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.The results in [22] describe the dynamics of calcium as a wave form.The wave-like shape can be clearly seen in Figure6.The 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 Ca2+ release due to channel opening.Figure7 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.


## 4.4. Simulation Procedure

TheHPNC 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 HPNC can be found in [42].Simulation ofHPNC in Snoopy is performed by first unfolding the HPNC model into the equivalent HPN one, followed by choosing a suitable hybrid simulation algorithm to execute the unfolded HPN, as presented in [25]. The hybrid simulators of Snoopy simulate the hybrid Petri nets based on the idea introduced in [43] by generating two random numbers. The 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 HPNC model in reasonable time. Therefore, 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], Ca2+ 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 Ca2+ 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]. The 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]. The integration time step adaptivity for the underlying ODEs system is achieved by an efficient ODE solver from [47]. The 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 this study in a cube form, by uniformly setting D to 5 for all experiments, which corresponds to 9μ × 9μ × 9μ.Table 3
Key measures for the three simulation experiments.
#CustersChannelsUnfolded placesUnfolded transitionsRun time∗ (minutes)1L=1,M=1,P=1K = 1, N = 1, Q = 12,92429,52470.082L=1,M=1,P=1K = 3, N = 4, Q = 23,06229,800260.533L=2,M=2,P=1K = 3, N = 2, Q = 23,21234,462497.73∗Done using OS X, 2.6 GHz Intel Core i7, 16 GB 1600 MHz DDR3.


## 5. Conclusions

This paper has presented the construction and simulation of 3D models usingHPNC. 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 HPNC model for intracellular calcium dynamics is created. The model is divided into two main parts: the continuous and the stochastic part. The continuous part describes the deterministic calcium release between cytosol, ER, and buffers as well as the diffusion process. The stochastic part models the clusters with the channel subunits. The 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. The flexibility of HPNC is showed by changing only a few parameters to run more than one experiment. The results of the developed HPNC have been validated against results known from the literature.OurHPNC 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 HPNC to other modelling scenarios in future work.



---

*Source: 4715172-2020-07-30.xml*